# Numerically solving 2 nonlinear PDEs of 2nd and 1st order

I want to compute flux coordinates $$\{\psi,\theta,\chi\}$$ as functions of cylindrical coordinates $$\{r,\theta,z\}$$ in the problem of ballooning mode instability in mirror traps (also called open traps) for plasma confinement and heating. Magnetic field in a mirror trap is expressed through the flux potentials as follows. In the cylindrical coordinates

$$\vec{B}(r,z)= \left\{\frac{1}{r}\frac{\partial \psi (r,z)}{\partial z},0,-\frac{1}{r}\frac{\partial \psi (r,z)}{\partial r}\right\}$$

or

$$\vec{B}(r,z)= \left\{G(r,z)\frac{\partial \chi (r,z)}{\partial r},0,G(r,z)\frac{\partial \chi (r,z)}{\partial z}\right\}$$

Function $$\psi(r,z)$$ can be found as a solution of the Grad-Shafranov equation, however method of searching of functions $$G$$ and $$\chi$$ is not widely known. Moreover, for investigating plasma stability it is better to express $$r$$, $$z$$ and $$G$$ as functions of $$\psi$$ and $$\chi$$. Appropriate equations were derived in 2011 by Arsenin and Terekhin:

(Eq.9)

$${\partial_\psi}\left(G(\psi ,\chi ) {\partial_\psi r^2(\psi ,\chi )}\right) +{\partial_\chi}\left(\frac{{\partial_\chi r^2(\psi,\chi)}}{G(\psi ,\chi ) r^2(\psi ,\chi )}\right)=0$$,

(Eq.15)

$$2 \pi \frac{\partial p(\psi )}{\partial \psi } \left(\left(G(\psi ,\chi ) {\partial_\psi {r^2}(\psi ,\chi )}\right)^2 +\frac{\left({\partial_\chi {r^2}(\psi ,\chi )}\right)^2}{{r^2}(\psi ,\chi )}\right)+\frac{\partial G(\psi ,\chi )^2}{\partial \psi }=0$$

Here $$r^2(\psi,\chi)$$ is the square of $$r(\psi,\chi)$$ and $$p(\psi)$$ is radial profile of the plasma pressure. These equations should be solved on a rectangle area

$$\Omega=\text{Rectangle}[\{0,0\},\{\psi_m,\chi_m\}]$$

with the boundary conditions

$$r^2(\psi,\chi)==0$$,

$$G(\psi,\chi)==1$$,

$$G(\psi,\chi)\ \partial_\psi r^2(\psi, \chi)=2/B_0(\chi)$$

at the axis $$\psi=0$$ of the axial symmetry, where$$B_0(\chi)$$ is the magnetic field on the axis. Instead of the last (Neuman) BC we can try posing Dirichlet BC at the opposite edge of the rectangle:

$$r^2(\psi_m, \chi)=2\psi_m/B_0(\chi)$$

At the boundaries $$\chi=0$$ and $$\chi=\chi_m$$ we have the following conditions:

$$\partial_\chi r^2(\psi,\chi)=0$$.

Formal solution of Eq.15 can be found with DSolve. It expresses $$G(\psi,\chi)$$ through a double integral of $$r^2(\psi,\chi)$$. Calculating the double integral might be time-consuming so finding both $$r2=r^2$$ and $$G$$ with the use of NDSolve seems preferable. Unfortunately, straightforward attaсk on the problem failed with the error message:

"The PDE coefficient ... does not evaluate to a numeric scalar
at the coordinate {0.0025000000000000005,2}; it evaluated to


EDIT #1

As it is said in FEMDocumentation/tutorial/FiniteElementProgramming, to solve nonlinear PDE it should be put in an inactive form. Here is my second attempt to solve the above formulated problem.

Here is inactive version of my code. Note that $$G2$$ and $$r2$$ are squares of $$G$$ and $$r$$ respectively.

Introduce some functions and parameters:

ClearAll[\[Psi], \[Chi], \[Theta], B, B0, G, r, z, r2]
(* set functions p[\[Psi]] and B0[\[Chi]]*)
funs = {B0 ->
Function[\[Chi],
B00 (1 + (K - 1) Sin[(\[Pi] \[Chi])/(2 \[Chi]m)]^2)]
, p -> Function[\[Psi], p0 (1 - \[Psi]/\[Psi]m)]}

(* set input parameters *)params = {\[Chi]m ->
B00 Sqrt[K] L, \[Psi]m -> B00 a^2/2, B00 -> 1, L -> 1, K -> 16,
a -> 0.1, p0 -> A/(4 \[Pi]), A -> 1}

(* define region where we look for a solution *)
\[CapitalOmega] = Rectangle[{0, 0}, {\[Psi]m, \[Chi]m}] //. params


Prepare functions r2 and G2 for setting Initial Seeding an BCsNote that we put into initial seeding and BCs a so called paraxial (i.e. near-axis solution of our equations):

(* Define functions for Initial Seeding and BCs *)
r2[\[Psi]_, \[Chi]_] = (2 \[Psi])/
B0[\[Chi]] - (
B0[\[Chi]] Derivative[B0][\[Chi]]^2 -
8 \[Pi] Derivative[p] -
B0[\[Chi]]^2 (B0^\[Prime]\[Prime])[\[Chi]])/(
2 B0[\[Chi]]^3) \[Psi]^2 /. funs //. params
G2[\[Psi]_, \[Chi]_] =
1 - (8 \[Pi] \[Psi] Derivative[p])/B0[\[Chi]]^2 /. funs //.
params


Set Eq. 9 in active and inactive forms and check that they are equivalent:

eq09 = \!$$\*SubscriptBox[\(\[PartialD]$$, $$\[Psi]$$]$$( \*SqrtBox[\(G2[\[Psi], \[Chi]]$$]
\*SubscriptBox[$$\[PartialD]$$, $$\[Psi]$$]r2[\[Psi], \[Chi]])\)\) + \
\!$$\*SubscriptBox[\(\[PartialD]$$, $$\[Chi]$$]$$( \*FractionBox[\( \*SubscriptBox[\(\[PartialD]$$, $$\[Chi]$$]r2[\[Psi], \[Chi]]\), $$\*SqrtBox[\(G2[\[Psi], \[Chi]]$$] r2[\[Psi], \[Chi]]\)])\)\) - 0;
\[DoubleStruckCapitalC]09 = {{Sqrt[G2[\[Psi], \[Chi]]], 0}, {0, 1/(
Sqrt[G2[\[Psi], \[Chi]]] r2[\[Psi], \[Chi]])}};
\[DoubleStruckCapitalC]09 // MatrixForm
\[DoubleStruckCapitalC]09 .
Eq09 = Inactive[Div][\[DoubleStruckCapitalC]09 .
r2[\[Psi], \[Chi]], {\[Psi], \[Chi]}], {\[Psi], \[Chi]}]
Activate[Eq09]
% == eq09


Set Eq. 15 in active and inactive forms and check that they are equivalent:

(* Eq. 15 from the above cited paper *)
eq15 = \!$$\*SubscriptBox[\(\[PartialD]$$, $$\[Psi]$$]$$G2[\[Psi], \[Chi]]$$\) +
2 \[Pi] (\!$$\*SubscriptBox[\(\[PartialD]$$, $$\[Psi]$$]$$p[\[Psi]]$$\)) \
(G2[\[Psi], \[Chi]] (\!$$\*SubscriptBox[\(\[PartialD]$$, $$\[Psi]$$]$$r2[\[Psi], \ \[Chi]]$$\))^2 + (\!$$\*SubscriptBox[\(\[PartialD]$$, $$\[Chi]$$]$$r2[\[Psi], \ \[Chi]]$$\))^2/r2[\[Psi], \[Chi]]) - 0
\[DoubleStruckCapitalC]15 = {{G2[\[Psi], \[Chi]], 0}, {0, 1/
r2[\[Psi], \[Chi]]}};
\[DoubleStruckCapitalC]15 // MatrixForm;

(* This version of Eq. 15 begin to work somehow *)
Eq15 = D[G2[\[Psi], \[Chi]], \[Psi]] + 2 \[Pi] \!$$\*SubscriptBox[\(\[PartialD]$$, $$\[Psi]$$]$$p[\[Psi]]$$\) Inactive[
Grad][r2[\[Psi], \[Chi]], {\[Psi], \[Chi]}] . \
\[DoubleStruckCapitalC]15 .

Eq15 = D[G2[\[Psi], \[Chi]], \[Psi]] + 2 \[Pi] \!$$\*SubscriptBox[\(\[PartialD]$$, $$\[Psi]$$]$$p[\[Psi]]$$\) Inactive[
r2[\[Psi], \[Chi]], {\[Psi], \[Chi]}] . \
\[DoubleStruckCapitalC]15 .

Activate[Eq15]
% == eq15 // Simplify


Set boundary conditions:

BCD17 = DirichletCondition[r2[\[Psi], \[Chi]] == 0, \[Psi] == 0];
BCD20 = DirichletCondition[
r2[\[Psi], \[Chi]] == r2[\[Psi], \[Chi]], \[Psi] == 0];
BCD21 = DirichletCondition[
G2[\[Psi], \[Chi]] == G2[\[Psi], \[Chi]], \[Psi] == 0];
BCN22 = NeumannValue[(* G[\[Psi],\[Chi]]\!$$\*SubscriptBox[\(\[PartialD]$$, $$\[Psi]$$]$$r2[\[Psi], \[Chi]]$$\)\
\[Equal] *)2/B0[\[Chi]], \[Psi] == 0];


Combine a system of PDEs:

(* Сombine system of PDEs *)
(* Version with Neuman and Dirichlet BCs on lower edge of
\[CapitalOmega] *)
PDEs = {Eq09 == BCN22, Eq15 == 0, BCD17, BCD21} /. funs //. params;

(* Version with Dirichlet BCs on lower and uppers edges of
\[CapitalOmega] *)
PDEs = {Eq09 == 0, Eq15 == 0, BCD17, BCD20, BCD21} /. funs //. params


Solve it:

"Solution"
NDSolveValue[PDEs, {r2,
G2}, {\[Psi], \[Chi]} \[Element] \[CapitalOmega],
InitialSeeding -> {r2[\[Psi], \[Chi]] == r2[\[Psi], \[Chi]],
G2[\[Psi], \[Chi]] == G2[\[Psi], \[Chi]]}]


Unfortunately this story ends with the following error:

FindRoot::dfmin: The minimal damping factor of 1/10000 has been reached.

NDSolveValue::fempsf: PDESolve could not find a solution.


P.S.

In my other Question How to pose Dirichlet and Neumann BCs on same boundary? I noted that NDSolve cannot solve Laplace equation if Dirichlet and Neuman BCs are posed on same edge of the rectangle region $$\Omega$$. This is also the case with the above problem. Therefor instead of Neuman BC

$$G(\psi,\chi)\partial_\psi r^2(\psi, \chi)=2/B_0(\chi)$$

at $$\psi=0$$ it is better to pose Dirichlet BC

$$r^2(\psi,\chi)=2\psi/B_0(\chi)$$

at $$\psi=\psi_m$$.

Any help would be greatly appreciated.

• I get a different error: NDSolveValue::femnlmdor::"The maximum derivative order of the nonlinear PDE coefficients for \ the Finite Element Method is larger than 1. It may help to rewrite \ the PDE in inactive form." Apr 13, 2021 at 5:56
• You are probably better off using the TensprProductGrid method. See the section "What Triggers the Use of the Finite Element Method" here Apr 13, 2021 at 5:59
• @user21: I have rewritten the PDEs in inactive form and got another error: FindRoot::dfmin: The minimal damping factor of 1/10000 has been reached. I don't understand how to bypass it. Apr 13, 2021 at 13:51
• That typically happens if the nonlinear solver can not solve the PDE. Try to solve a linear problem (eliminate the nonlinear terms) and use the result of that as a initial guess. Probably if would be useful if you could edit your post with an update of the new code that you have. Apr 13, 2021 at 13:57
• @user21: Please look at 2nd version of my code. I have rewrittem PDEs in inactive form. Apr 14, 2021 at 7:32

To solve a system of nonlinear equations with FEM we use the method of the false transient. First we make replacement $$x=\psi, y=\chi$$, and remove all parameters with

funs = {B0 -> Function[y, B00 (1 + (k - 1) Sin[(\[Pi] y)/(2 ym)]^2)],
p -> Function[x, p0 (1 - x/xm)]}; params = {ym -> B00 Sqrt[k] L,
xm -> B00 a^2/2, p0 -> A/(4 \[Pi])}; num = {B00 -> 1, L -> 1,
k -> 16, a -> 0.1, A -> 1};


For example Eq15/. funs /. params /. num. Second, we remove singularity at x=0 with small x0=10^-6. The solution of the problem converges quickly and finally we have

Clear["Global*"]
Needs["NDSolveFEM"]; x0 = 10^-6;
reg = Rectangle[{x0, 0}, {0.005, 4}]; dt = 1/10; R2[x_, y_] := 1.;
g2[x_, y_] := 1.;
eq15 = Derivative[1, 0][G2][x, y] -
100*(D[R2[i - 1][x, y], y]^2/R2[i - 1][x, y] +
G2[x, y]*D[R2[i - 1][x, y], x]^2);
bc = {DirichletCondition[{r2[x,
y] == (2*x)/(1 + 15*Sin[(Pi*y)/8]^2)}, x == x0],
DirichletCondition[r2[x, y] == (2*x)/(1 + 15*Sin[(Pi*y)/8]^2),
x == 0.005]};
eq9 = -((D[g2[i - 1][x, y], y]*Derivative[0, 1][r2][x, y])/(2*
g2[i - 1][x, y]^(3/2)*R2[i - 1][x, y])) -
(Derivative[0, 1][r2][x, y]*
D[R2[i - 1][x, y], y])/(Sqrt[g2[i - 1][x, y]]*
R2[i - 1][x, y]^2) +
Derivative[0, 2][r2][x,
y]/(Sqrt[g2[i - 1][x, y]]*
R2[i - 1][x, y]) + (D[g2[i - 1][x, y], x]*
Derivative[1, 0][r2][x, y])/
(2*Sqrt[g2[i - 1][x, y]]) +
Sqrt[g2[i - 1][x, y]]*
Derivative[2, 0][r2][x, y] + (r2[x, y] - R2[i - 1][x, y])/dt;

Do[R2[i] =
NDSolveValue[{-((D[g2[i - 1][x, y], y]*
Derivative[0, 1][r2][x, y])/(2*g2[i - 1][x, y]^(3/2)*
R2[i - 1][x, y])) - (Derivative[0, 1][r2][x, y]*
D[R2[i - 1][x, y], y])/(Sqrt[g2[i - 1][x, y]]*
R2[i - 1][x, y]^2) +
Derivative[0, 2][r2][x,
y]/(Sqrt[g2[i - 1][x, y]]*
R2[i - 1][x, y]) + (D[g2[i - 1][x, y], x]*
Derivative[1, 0][r2][x, y])/(2*Sqrt[g2[i - 1][x, y]]) +
Sqrt[g2[i - 1][x, y]]*
Derivative[2, 0][r2][x, y] + (r2[x, y] - R2[i - 1][x, y])/dt ==
0, bc}, r2, Element[{x, y}, reg], Method -> {"FiniteElement"}];
g2[i] = NDSolveValue[{Derivative[1, 0][G2][x, y] -
100*(D[R2[i - 1][x, y], y]^2/R2[i - 1][x, y] +
G2[x, y]*D[R2[i - 1][x, y], x]^2) == 0,
DirichletCondition[G2[x, y] == 1, x == x0]}, G2,
Element[{x, y}, reg], Method -> {"FiniteElement"}] // Quiet;, {i,
10}]


Visualization of two last iterations shows that numerical solution converges

Table[{ContourPlot[R2[i][x, y], Element[{x, y}, reg],
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotRange -> All, Contours -> 10],
ContourPlot[g2[i][x, y], Element[{x, y}, reg],
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
Contours -> 10]}, {i, 9, 10}] The implementation of two boundary conditions at $$\psi=0$$ is given by

bc = {NeumannValue[- Sqrt[g2[i - 1][x, y]] 2/(1 + 15*Sin[(Pi*y)/8]^2),
x == x0],
DirichletCondition[r2[x, y] == (2*x)/(1 + 15*Sin[(Pi*y)/8]^2),
x == 0.005 || x == x0]};
Do[R2[i] =
NDSolveValue[{-((D[g2[i - 1][x, y], y]*
Derivative[0, 1][r2][x, y])/(2*g2[i - 1][x, y]^(3/2)*
R2[i - 1][x, y])) - (Derivative[0, 1][r2][x, y]*
D[R2[i - 1][x, y], y])/(Sqrt[g2[i - 1][x, y]]*
R2[i - 1][x, y]^2) +
Derivative[0, 2][r2][x,
y]/(Sqrt[g2[i - 1][x, y]]*
R2[i - 1][x, y]) + (D[g2[i - 1][x, y], x]*
Derivative[1, 0][r2][x, y])/(2*Sqrt[g2[i - 1][x, y]]) +
Sqrt[g2[i - 1][x, y]]*
Derivative[2, 0][r2][x, y] + (r2[x, y] - R2[i - 1][x, y])/dt ==
bc[], bc[]}, r2, Element[{x, y}, reg],
Method -> {"FiniteElement"}];
g2[i] = NDSolveValue[{Derivative[1, 0][G2][x, y] -
100*(D[R2[i - 1][x, y], y]^2/R2[i - 1][x, y] +
G2[x, y]*D[R2[i - 1][x, y], x]^2) == 0,
DirichletCondition[G2[x, y] == 1, x == x0]}, G2,
Element[{x, y}, reg], Method -> {"FiniteElement"}] // Quiet;, {i,
10}]


The numerical solution is not much differ from above.

• Thank you for your solution! Unfortunatly G2 has unphysical huge peak near x=x0. For the monent I don'r know why it is so. Also, I wonder Is there a secret meaning to write Derivative [1, 0] [G2] [x, y] instead of D[G2[x, y],x] in eq15? Apr 17, 2021 at 10:43
• @IgorKotelnikov All "unphysical" results come from boundary conditions. Therefore, we can play with bc to get physical results. The secret meaning is that we convert all expressions in Raw Input Form as it required for this Forum :) Apr 17, 2021 at 10:57
• Actually, I developed a similar code but solved 2 equations at once by single call to NDSolve. Why did you make 2 calls to NDSolve? In my method, I see unphysical negative r2 near axis x=x0 with same BCs. Another difference is that I did not load Needs["NDSolveFEM"];. I don't see where you have used this package/ Apr 17, 2021 at 11:41
• @IgorKotelnikov We integrate equations with Method -> {"FiniteElement"}, but second equation for G2 is not defined well for this method, it is why I put // Quiet in the last line in Do loop. Apr 17, 2021 at 12:51
• Thank you so much again. I think that "unphysical" solution is the result of ill BCs. BCs for $r2$ should be $r2=0$ and $\partial_\psi r2=2/B_0(\chi)$ at $\psi=0$. Is it possible to pose 2 BCs on same edge of reg? Apr 18, 2021 at 10:43