NDSolve Stiffness and Strange Behaviour with Dirichlet Conditions

Question

I working with the following semiconductor potential equilibrium equation for a diode

$$\psi(x)''=e^{\psi(x)}-e^{-\psi(x)}+10^6\text{sgn}(x-(l+r)/2)$$

where $l$ and $r$ are left and right boundaries respectively. The boundary values are chosen exactly so that $\psi''(l)=\psi''(r)=0$. Or

$$e^{\psi(l)}-e^{-\psi(l)}=10^6\implies \psi(l)=\sinh ^{-1}(5\cdot 10^5)\text{ also }\psi(r)=-\psi(l)$$

This is essential from preventing the equation from blowing up

Here's the NDSolve workflow

{l, r} = {0, 1};
bval = ArcSinh[500000];
ode = f''[x] == 2 Sinh[f[x]] + 10^6 Sign[x - (l + r)/2];
F = NDSolveValue[{ode, f[l] == bval, f[r] == -bval}, f, {x, l, r}]

Which gives a stiffness warning

but when I specify the boundary conditions like this

F = NDSolveValue[{ode, DirichletCondition[f[x] == bval, x == l], 
   DirichletCondition[f[x] == -bval, x == r]}, 
  f, {x, l, r}]

I get a solution with no warnings. And it seems not very accurate i.e the PDE residual is far from zero

Plot[{F''[x] - (E^(F[x]) - E^(-F[x]) + 
     10^6 Sign[x - (left + right)/2])}, {x, l, r}, 
 PlotRange -> Full]

If I change the scale to say {l, r} = {0, .02/0.0033837837}; (a number from the actual physics)

I get another error and a non-sense solution from the second method

Can anyone help clear up whats going on? I'm not even looking for a miraculous solution I just want to understand this behaviour.

Thanks

Answer 1

If the solution is received with options Automatic , then checking the number of elements shows

F["ElementMesh"]
(*Out[]:ElementMesh[{{0., 1.}}, {LineElement["<" 20 ">"]}]*)

Use the options here.

<< NDSolve`FEM`
l = 0; r = 1;
bval = ArcSinh[500000]; ode = 
 f''[x] == 2 Sinh[f[x]] + 10^6 Sign[x - (l + r)/2]; F = 
 NDSolveValue[{ode, DirichletCondition[f[x] == bval, x == l], 
   DirichletCondition[f[x] == -bval, x == r]}, f, {x, l, r}, 
  Method -> {"FiniteElement", "InterpolationOrder" -> {f -> 2}, 
    "MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}]

Plot[F[x], {x, 0, 1}]

Plot[{F''[x] - (E^(F[x]) - E^(-F[x]) + 10^6 Sign[x - (l + r)/2])}, {x,
   l, r}, PlotRange -> Full]

Now check

F["ElementMesh"]

(*Out[]= ElementMesh[{{0., 1.}}, {LineElement["<" 10000 ">"]}]*)

Thus, a poor solution to the equation was obtained on a grid of 20 elements.