Solve PDE system for Helmholtz decomposition

Question

In the article Analytical Helmholtz Decomposition and Potential Functions for many n-dimensional unbounded vector fields on page 3 presents the formulas (6-8).

If I understand these formulas correctly, then for example, for a two-dimensional vector field $f(x,y)=\langle f_1(x,y),f_2(x,y) \rangle$, this system of equations will look like this:

$\begin{cases} f_1(x,y) = g_1(x,y)+r_1(x,y) \\ f_2(x,y) = g_2(x,y)+r_2(x,y) \\ g_1(x,y)=\frac{\partial}{\partial x} G(x,y) \\ g_2(x,y)=\frac{\partial}{\partial y} G(x,y) \\ \frac{\partial}{\partial x}r_1(x,y)+\frac{\partial}{\partial y}r_2(x,y)=0 \end{cases}$

We can simplify this expression as follows:

\begin{cases} f_1(x,y) = \frac{\partial}{\partial x} G(x,y)+r_1(x,y) \\ f_2(x,y) = \frac{\partial}{\partial y} G(x,y)+r_2(x,y) \\ \frac{\partial}{\partial x}r_1(x,y)+\frac{\partial}{\partial y}r_2(x,y)=0 \end{cases}

This is a system of partial differential equations, the unknowns are: $[G(x,y),r_1(x,y),r_2(x,y)]$.

I will try to solve this system of equations for a 2D vector field $f=\langle x,y \rangle$:

f = {x, y};

eqns = {f[[1]] == D[G[x, y], x] + Subscript[r, 1][x, y], 
   f[[2]] == D[G[x, y], y] + Subscript[r, 2][x, y], 
   D[Subscript[r, 1][x, y], x] + D[Subscript[r, 2][x, y], y] == 0};

DSolve[eqns, {G[x, y], Subscript[r, 1][x, y], 
   Subscript[r, 2][x, y]}, {x, y}];

But Mathematica could not solve such a system:

I'm trying to figure this out, but my experience with PDE isn't good enough. What is the problem of the impossibility of obtaining solutions to such a PDE system connected with: incorrect notation of differential operators, system is redefined, has an unsolvable structure, needs additional transformation?

I would be glad for a hint.

EDIT №2: Helmholtz decomposition of a 4-dimensional vector field

(***For coordinates***)

coords = {x, y, z, u};

(***Given a 4-D vector field***)

f = {y + x, -x - u^2, 1, z^3 + 3 y};

(***We make a replacement and find the quasipotential by the method \
proposed in https://math.stackexchange.com/a/2056657/656085 ***)

F = ReplaceAll[f, {x -> t x, y -> t y, z -> t z, u -> t u}];

(***Then we find the quasipotential.For a conservative vector \
field,it coincides with the true potential function***)

G = Integrate[F.coords, {t, 0, 1}] // Expand

(***Form the vector R - divergence-free component of the Helmholtz \
decomposition***)

R = Table[r[i][x, y, z, u], {i, 1, 4}];

Next, we use the Helmholtz expansion in a new form:

$\textbf{f}=\textbf{Grad} f(x,y,z,u)+\textbf{R}(x,y,z,u)+\delta r(x,y,z,u)$

here $\delta r(x,y,z,u)$ is a function that is added to each of the equations of the system of partial differential equations for the Helmholtz expansion. This addition allows solving this system of equations.

eqns = Thread[Grad[G, coords] + \[Delta]r[x, y, z, u] + R == f] //Expand

(***Begin to solve the system of equations.We express R***)

solR = Solve[eqns, R] // Expand // First

(***Substitute into the divergence-free condition***)

Div[solR[[All, 2]], coords] == 0

(***And solve the equation for \[Delta]r[x,y,z,u]***)

aux = DSolve[Div[solR[[All, 2]], coords] == 0, {\[Delta]r[x, y, z, u]}, 
coords] // Expand // First

(***We substitute and check if the divergence-free component meets \
the condition***)

Div[ReplaceAll[solR[[All, 2]], aux[[1]]], coords] == 0

(***Checking Curl of the gradient component of the vector field.It \
is conservative.***)

Curl[Grad[ReplaceAll[G + \[Delta]r[x, y, z, u], aux[[1]]], coords], coords] // Simplify // Normal

It seems to me that this component $\delta r(x,y,z,u)$ somehow modifies both the $\textbf{R}$ and the gradient component $\textbf{Grad} f(x,y,z,u)$. To what it would be more correct to attribute it, I do not know yet.

Answer 1

It's well known that, though DSolve has been improving in recent versions, it's still far beyond perfect, so let's help it a bit. Obviouly, $r_1$ and $r_2$ can be easily eliminated:

f = {x, y};
eqns = {f[[1]] == D[G[x, y], x] + Subscript[r, 1][x, y], 
        f[[2]] == D[G[x, y], y] + Subscript[r, 2][x, y], 
        D[Subscript[r, 1][x, y], x] + D[Subscript[r, 2][x, y], y] == 0};

rulerxy = Solve[eqns // Most, {Subscript[r, 1][x, y], Subscript[r, 2][x, y]}][[1]]
(* {Subscript[r, 1][x, y] -> x - Derivative[1, 0][G][x, y], 
    Subscript[r, 2][x, y] -> y - Derivative[0, 1][G][x, y]} *)

{funcr1[x_, y_], funcr2[x_, y_]} = 
  {Subscript[r, 1][x, y], Subscript[r, 2][x, y]} /. rulerxy;

neweqn = Last@eqns /. {Subscript[r, 1] -> funcr1, Subscript[r, 2] -> funcr2}
(* 2 - Derivative[0, 2][G][x, y] - Derivative[2, 0][G][x, y] == 0 *)

This is Poisson's equation, and DSolve can find a general solution for it:

ruleG = DSolve[neweqn, G, {x, y}][[1]]
(* {G -> Function[{x, y}, x^2 + C[1][I x + y] + C[2][-I x + y]]} *)

C[1] and C[2] are arbitrary functions, we can substitute in any function to obtain a specific Helmholtz decomposition e.g.:

test = {G[x, y], {Subscript[r, 1][x, y], Subscript[r, 2][x, y]}} /. 
    rulerxy /. ruleG /. C[_] -> (0 &)
(* {x^2, {-x, y}} *)

(* Check: *)
Grad[#, {x, y}] + #2 & @@ test
(* {x, y} *)

Remark

Though it happens to work in this case, usually the "first general solution, then particular solution" strategy won't work well in PDE analysis AFAIK. The more common approach is to directly solve for particular solution satisfying certain initial/boundary conditions.