3
$\begingroup$

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 enter image description here

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}]

enter image description here

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]

enter image description here

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

enter image description here

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

$\endgroup$
1
$\begingroup$

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]

Figure 1

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ What is the purpose of "InterpolationOrder"->2? $\endgroup$ – user21 Dec 18 '19 at 8:55
  • $\begingroup$ Just for the memory of what is used here. You can remove this option. $\endgroup$ – Alex Trounev Dec 18 '19 at 10:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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