2
$\begingroup$

In attempting to solve a fourth order system I have encountered an issue with the way NDSolve uses boundary conditions. Consider the following three different attempts:

Nep = 1/2;
Nhp = 2;
Nen = 2;
Nhn = 1/2;
xmin = -2;
xmax = 2;

sol1 = NDSolve[{2 ϕ''[x] == 
     Ne[x] - Nh[x] - UnitStep[x] (Nen - Nhn) - 
      UnitStep[-x] (Nep - Nhp), -Nh[x] ϕ'[x] - Nh'[x] == 
     0, -Ne[x] ϕ'[x] + Ne'[x] == 0, ϕ[xmin] == 
     0, ϕ'[xmin] == 0, ϕ'[xmax] == 0, 
    Ne[xmin] Nh[xmin] == 1}, {ϕ, Ne, Nh}, {x, xmin, xmax}];
Plot[Evaluate[ϕ[x] /. sol1], {x, xmin, xmax}, PlotRange -> All, 
 AxesLabel -> {"\!\(\*OverscriptBox[\(x\), \(~\)]\)", 
   "\!\(\*OverscriptBox[\(ϕ\), \(~\)]\)"}]

sol2 = NDSolve[{2 ϕ''[x] == 
     Ne[x] - Nh[x] - UnitStep[x] (Nen - Nhn) - 
      UnitStep[-x] (Nep - Nhp), -Nh[x] ϕ'[x] - Nh'[x] == 
     0, -Ne[x] ϕ'[x] + Ne'[x] == 0, ϕ[0] == 
     0, ϕ'[xmin] == 0, ϕ'[xmax] == 0, 
    Ne[xmin] Nh[xmin] == 1}, {ϕ, Ne, Nh}, {x, xmin, xmax}];
Plot[Evaluate[ϕ[x] /. sol2], {x, xmin, xmax}, PlotRange -> All, 
 AxesLabel -> {"\!\(\*OverscriptBox[\(x\), \(~\)]\)", 
   "\!\(\*OverscriptBox[\(ϕ\), \(~\)]\)"}]

sol3 = NDSolve[{2 ϕ''[x] == 
     Ne[x] - Nh[x] - UnitStep[x] (Nen - Nhn) - 
      UnitStep[-x] (Nep - Nhp), -Nh[x] ϕ'[x] - Nh'[x] == 
     0, -Ne[x] ϕ'[x] + Ne'[x] == 0, ϕ[xmax] == 
     0, ϕ'[xmin] == 0, ϕ'[xmax] == 0, 
    Ne[xmin] Nh[xmin] == 1}, {ϕ, Ne, Nh}, {x, xmin, xmax}];
Plot[Evaluate[ϕ[x] /. sol3], {x, xmin, xmax}, PlotRange -> All, 
 AxesLabel -> {"\!\(\*OverscriptBox[\(x\), \(~\)]\)", 
   "\!\(\*OverscriptBox[\(ϕ\), \(~\)]\)"}]

The difference between the three is just in the position at which a value of $\phi$ is specified. Given that there is no explicit dependence on $\phi$ in the system of equations, this shouldn't affect the solution up to the addition of a constant. However, only the second, where $\phi$ is specified at $x=0$ runs successfully. Why is this happening when all three problems are well-posed? Is there a fix?

Thanks in advance for any help.

$\endgroup$
2
  • $\begingroup$ In my run, your sol1 also successed by with a warning throwed out, sol3 failed... $\endgroup$
    – Harry
    Nov 3 '20 at 8:29
  • $\begingroup$ I'm not really counting that as a success when it isn't giving the correct answer. $\endgroup$
    – Chris
    Nov 3 '20 at 8:43
3
$\begingroup$

UnitStep causes the problems. Substitute it by a very narrow ArcTan and it works fine. Like

sol1 = First@NDSolve[{2 \[Phi]''[x] == Ne[x] - Nh[x] - 3 (1/Pi ArcTan[x/10^-7]), -Nh[x] \[Phi]'[x] - Nh'[x] == 0, -Ne[x] \[Phi]'[x] + Ne'[x] == 0, \[Phi][xmin] == 0, \[Phi]'[xmin] == 0, \[Phi]'[xmax] == 0,  Ne[xmin] Nh[xmin] == 1}, {\[Phi], Ne, Nh}, 
                        {x, xmin, xmax}] 

May be you need higher WorkingPrecision.

$\endgroup$
2
$\begingroup$

I guess it's again the "Shooting" method that's not robust enough. The problem can be circumvented with the help of new-in-12 nonlinear FiniteElement method.

Since FiniteElement is chosen, the pre-processor of NDSolve can no longer handle the equation system (at least now), so let's help it a bit:

Nep = 1/2;
Nhp = 2;
Nen = 2;
Nhn = 1/2;
xmin = -2;
xmax = 2;

midsol=DSolve[{-Ne[x] ϕ'[x] + Ne'[x] == 0, -Nh[x] ϕ'[x] - Nh'[x] == 0}, {Nh, Ne},x][[1]]

csol = Solve[1 == Ne@xmin Nh@xmin /. midsol, C[1]][[1]]

neweq = 2 ϕ''[x] == 
     Ne[x] - Nh[x] - UnitStep[x] (Nen - Nhn) - UnitStep[-x] (Nep - Nhp) /. midsol /. 
   csol /. C[2] -> c
(*
2 ϕ''[x] == -c E^-ϕ[x] + E^ϕ[x]/c + (3 UnitStep[-x])/2 - (3 UnitStep[x])/2
 *)

psol = ParametricNDSolveValue[neweq, ϕ, {x, xmin, xmax}, c, 
  Method -> "FiniteElement"]

ϕ'[xmin] == 0, ϕ'[xmax] == 0 is omitted, because they're the default setting of FiniteElement in this case, you may check the document of NeumannValue for more info.

para1 = FindRoot[psol[c][xmin] == 0, {c, 1}]
(* {c -> 1.70937} *)

para2 = FindRoot[psol[c][0] == 0, {c, 1}]
(* {c -> 1.} *)

para3 = FindRoot[psol[c][xmax] == 0, {c, 1}]
(* {c -> 0.585012} *)

Plot[psol[c][x] /. {para1, para2, para3} // Evaluate, {x, xmin, xmax}]

enter image description here

$\endgroup$

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.