# NDSolve Stiffness and Strange Behaviour with Dirichlet Conditions

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

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.

<< NDSolveFEM
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.

• What is the purpose of "InterpolationOrder"->2? Dec 18 '19 at 8:55
• Just for the memory of what is used here. You can remove this option. Dec 18 '19 at 10:01