Using NDSolve Iteratively to Solve a Boundary Value Problem

I am trying to solve a simple nonlinear boundary value problem with FEM using the in built routine NDSolveValue. What I am doing is painfully slow, and I am sure you know how to do better. The equations I want to solve are way more complicated, so this one just serves as an illustration of what I want to do.

The equation in question takes the following form

$\partial^2_x\psi+\partial^2_{y}\psi=e^{\psi}$

with the integration domain being $(x,y)\in[0,1]\times[0,1]$. The boundary conditions are pretty simple:

$\partial_x \psi =0$ at $x=1$

and

$\psi=0$ at x =0, y=0 and y =1.

First, I create a grid:

<< NDSolveFEM
Om = Rectangle[];
Bmesh = ToBoundaryMesh[Om,MaxCellMeasure -> {"Length" -> 0.01}];
mesh = ToElementMesh[Om,MaxCellMeasure -> {"Length" -> 0.01}];
pts = mesh["Coordinates"];


Then I solve the linear problem iteratively:

u0[x_, y_] = 0;
u0data = Apply[u0, pts, 2];
Errmax = 10;
Errmin = 10^-6;
dE = (Errmax + Errmin)/2;
Print["Error = ", dE // Dynamic]
While[Errmin < dE < Errmax,
du = NDSolveValue[{
\!$$\*SuperscriptBox[\(du$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "2"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] +
\!$$\*SuperscriptBox[\(du$$,
TagBox[
RowBox[{"(",
RowBox[{"2", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] - E^u0[x, y] du[x, y] == -(\!$$\*SubscriptBox[\(\[PartialD]$$, $$x, x$$]$$u0[x, y]$$\) + \!$$\*SubscriptBox[\(\[PartialD]$$, $$y, y$$]$$u0[x, y]$$\) -
Exp[u0[x, y]]) + NeumannValue[0, x == 1],
DirichletCondition[
du[x, y] == 0, ({x, y} \[Element] Bmesh) && x != 1]},
du, {x, y} \[Element] mesh,
Method -> {"PDEDiscretization" -> {"FiniteElement"}}];
dudata = Apply[du, pts, 2];
u0data = u0data + dudata;
u0 = ElementMeshInterpolation[{mesh}, u0data];
dE = Norm[dudata, Infinity];
Clear[du, dudata];
];
ufinal = u0;
Clear[u0, Errmin, Errmax, dE]


The code is just to slow (I mean, I have implemented this in Fortran spectral collocation methods and using a simple Newton's method and runs in seconds). Is there a way to speed this up in Mathematica?

• I recall a similar question where ElementMeshInterpolation was the bottleneck... – Henrik Schumacher Dec 13 '17 at 0:05
• Is the integration domain still a square in your actual problem? – xzczd Dec 13 '17 at 5:46
• No, my integration domain will be more complicated, that is why I am using finite elements. – user12588 Dec 14 '17 at 8:52

If I understand your question right, you try to find an iterative solution of your pde:

NestList[NDSolveValue[{Laplacian[\[Psi][x, y], {x, y}] ==Exp[#] (1 - # + \[Psi][x, y]) + NeumannValue[0, x == 1], \[Psi][0, y] == 0, \[Psi][x, 0] == 0, \[Psi][x, 1] ==0}, \[Psi][x, y], {x, 0, 1}, {y, 0, 1}] &, 0, 5]


the result (3 iterations) is obtained in 0.3s ( error 10^-9).

• Is there any further interest for an answer of the question ??? – Ulrich Neumann Dec 15 '17 at 17:27
• Your answer looks great, the only catch is that I cannot control the error, and that is not good for my purposes. But thank you for your time! – user12588 Jan 4 '18 at 8:02
• @user12588: Thank you for your feedback. error control would be the same as you used in your approach – Ulrich Neumann Jan 4 '18 at 8:20

In version 12.0 this is straight forward:

sol = NDSolveValue[{Laplacian[u[x, y], {x, y}] - E^u[x, y] == 0,
DirichletCondition[u[x, y] == 0, x != 1]},
u, Element[{x, y}, Rectangle[]]];
Plot3D[sol[x, y], Element[{x, y}, Rectangle[]]]