4
$\begingroup$

I am trying to solve a second-order elliptic non-linear PDE on a disk, and it seems that the "solution" given by NDSolve (more precisely, NDSolveValue) is either incorrect or is only correct at C^0 level (so its derivatives are inaccurate). In particular, it seems that if I plug the "solution" in the PDE I do not get zero -- in fact, I get O(1) errors!

Please consider the following example. The PDE is

usol = NDSolveValue[{\!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y}\), \(2\)]\(u[x, y]\)\) + 
    0.5 Exp[-u[x, y]] - 2. Exp[u[x, y]] == 0., 
  DirichletCondition[u[x, y] == 2. Log[1. - 0.5 Sin[2. ArcTan[x, y]]],
    True]}, u, {x, y} \[Element] Disk[]]

The plot of the solution looks fine

Plot3D[usol[x, y], {x, y} \[Element] Disk[], PlotRange -> Full, 
 AxesLabel -> Automatic]

Now I apply the PDE on the solution

Eqsol[x_, y_] := 
 Derivative[2, 0][usol][x, y] + Derivative[0, 2][usol][x, y] + 
  0.5 Exp[-usol[x, y]] - 2. Exp[usol[x, y]]

And if I plot it I get

Plot3D[Eqsol[x, y], {x, y} \[Element] Disk[], PlotRange -> Full, 
 AxesLabel -> Automatic]

Since this should vanish for an actual solution, we have an O(1) discrepancy. In an average sense things are just slightly better, but still not good, as we can see by integrating the square of Eqsol,

NIntegrate[Eqsol[x, y]^2, {x, y} \[Element] Disk[]]
0.128622

So what could be happening? Am I doing something wrong? How can I obtain accurate derivatives of the solution of the PDE (ideally up to second order, but I only really need first derivatives for my main purposes)? Thank you.

Note: I tried playing with AccuracyGoal, PrecisionGoal and InterpolationOrder without success. I also tried other equations, boundary conditions and domain shapes. The only case that I got really nice results (Eqsol ~ 10^-12) was for a linear equation, with sufficiently simple boundary conditions (such as x^2- y^2) in a rectangular domain.

Improve this question
$\endgroup$
0

1 Answer 1

Reset to default
3
$\begingroup$

You can try with finer discretization by making MaxCellMeasure smaller:

usol = NDSolveValue[{\!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y}\), \(2\)]\(u[x, y]\)\) + 
          0.5 Exp[-u[x, y]] - 2. Exp[u[x, y]] == 0., 
      DirichletCondition[
     u[x, y] == 2. Log[1. - 0.5 Sin[2. ArcTan[x, y]]],
         True]}, u, {x, y} \[Element] Disk[], 
   Method -> {FiniteElement, 
     MeshOptions -> {MaxCellMeasure -> 0.0001, 
       MeshElementType -> TriangleElement}}];

Eqsol[x_, y_] := Derivative[2, 0][usol][x, y] + Derivative[0, 2][usol][x, y] +
   0.5 Exp[-usol[x, y]] - 2. Exp[usol[x, y]];

Plot3D[Evaluate@Eqsol[x, y], {x, y} \[Element] Disk[], PlotRange -> Full, AxesLabel -> Automatic]

Finer mesh Error

However, please read this Mathematica FEM Tutorial about why your way of verifying the result of NDSolveValue gives you falsely exaggerated error. It has to do with calculating the derivative of interpolating function, which is not numerically stable.

Improve this answer
$\endgroup$
4
  • 2
    $\begingroup$ Nice answer. One point: It’s not “stability” per se, although magnification of error is indeed an issue with numerical differentiation. In this case, the derivative should be the exact derivative of an order-2 interpolation, which is not even continuous but piecewise constant. The order of the solution must be higher than the order of the derivative to get a passable approximation to the derivative. That’s not possible with FEM (currently). (+1) $\endgroup$
    – Michael E2
    Commented Jul 10, 2021 at 17:10
  • $\begingroup$ Thanks for the answer. Using a smaller "MaxCellMeasure" seems to help. I was indeed thinking that the problem was about taking derivatives of the interpolated function generated by NDSolve. Somehow there must be some very fine "wigglyness" in usol such that |usol - uexact| ~ 0, but |usol' - uexact'| is inaccurate (where uexact is the exact solution of the PDE), is that right? $\endgroup$
    – RAS
    Commented Jul 10, 2021 at 19:42
  • $\begingroup$ The thing is that I really need a way to compute derivatives of the solution up to first order (and maybe second order). I suppose that one possibility would be to define the derivatives numerically, like (f[x+a]-f[x])/a, taking a to be small but still larger than the mesh size in NDSolve. But is there a more efficient/elegant way to do it? $\endgroup$
    – RAS
    Commented Jul 10, 2021 at 19:42
  • $\begingroup$ I think the first derivatives of the solution look just fine: Plot3D[Evaluate@{Derivative[1, 0][usol][x, y], Derivative[0, 1][usol][x, y]}, {x, y} \[Element] Disk[]] $\endgroup$
    – Domen
    Commented Jul 10, 2021 at 19:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.