# Solve PDE system for Helmholtz decomposition

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.

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.

• Thanks! I tested code for several linear vector fields. In this case, it is not even necessary for them to specify the boundary conditions. But let's take a vector field as an example from here. In this case, the solution has not been found. Solution of the Helmholtz decomposition must be this potential, and the divergence-free term $r$ is $0$. I assume that the lack of a solution is due either to the unsolvable structure of the PDE, or to the DSolve algorithm, or boundary conditions can help. Are my thoughts correct?
– ayr
Commented Nov 20, 2022 at 4:07
• @dtn When b.c. is not imposed, Helmholtz decomposition is not unique: physics.stackexchange.com/q/10522/18778 The solution on that page essentially impose the constraint $\mathbf{r}=\mathbf{0}$. If you do the same on the solution in my answer, you'll obtain the same result. Commented Nov 20, 2022 at 4:42
• Yes, Helmholtz decomposition is not unique, all depends on the choice of boundary/initial conditions. I have one more question. When we compose a system of high-dimensional equations, the curl turns into a tensor and the system of equations becomes overdetermined. There can be two options. The first is to define more boundary/initial conditions in order to achieve a "balance". The second could be the "anti-curl" operator (math.stackexchange.com/questions/81405/anti-curl-operator). I didn't quite understand how to apply it in high dimensional. Let's try?
– ayr
Commented Nov 21, 2022 at 7:05
• @dtn Er… what do you mean by "When we compose a system of high-dimensional equations, the curl turns into a tensor and the system of equations becomes overdetermined."? High dimensional cases are just similar to 2D cases e.g. $\mathbf{f}=(f_1,f_2,f_3)$ results in 4 equations for $G$, $r_1$, $r_2$, $r_3$, aren't they? Commented Nov 21, 2022 at 8:41
• @dtn You're solving for the $\mathbf{R}$ where $\mathbf{r}=\nabla \times \mathbf{R}$? If so, then this $\mathbf{R}$ is not a vector (in narrow sense) but a 4×4 matrix. Try field = Table[f[i, j][x, y, z, u], {i, 4}, {j, 4}]; Curl[field, {x, y, z, u}]. This topic is a bit involved, I suggest reading this post: mathematica.stackexchange.com/a/191445/1871 Commented Nov 21, 2022 at 10:54