# 1D ODE: singularity or stiff system suspected

I'm trying to solve the nonlinear Poisson-Boltzmann equation in a cylindrical geometry. In this example I am surrounding a central charged rod with a neutralizing shell of charge, although I get the same errors if I don't have the neutralizing shell (the shell adds a discontinuity). I have also included code for the linearized version, so I have a good idea of what the solution looks like, and there really isn't anything weird going on like singularities or even oscillations. The code works fine for central charge densities up to -1.9 (as below) but fails as soon as it bump it to -2. I'd like to get it up to the -6 or -7 range. The linearized version uses a boundary conditions of phi=0 at infinity, but the nonlinear version uses a Neumann condition at 3. If I could push that out to 5 or so, that would help, but the code also fails if the outer boundary is more distant than 3 or so. I tried changing the WorkingPrecision to 30, and the code took much longer to run (10-15 seconds) but failed at the same charge density. Are there any other integration methods I should try? Or is there a way to feed it the linearized version as an initial guess?

Here is the code. This version works, but bumping sD to -2.0 causes a failure.

Module[{sD = -1.9, shellsize = 1.784}, Show[
(*this first module is the nonlinear version*)
Module[{sigmainner = sD, csaltM = 0.16, Rinner = 1., Router = 3.,
Rshell = shellsize, ShellQ = -1, e0, func, Lbjerrum = 0.701,
kappa, csalt, Yinner, Youter, Yshell, rhoshell},
csalt =
csaltM (6.022*10^23)/
0.001 (1/10^9)^3(*convert Molar units to nm^-3*);
kappa = Sqrt[
8 Pi csalt Lbjerrum];(*Debye screening parameter in inverse nm*)
Yinner = kappa Rinner; Youter = kappa Router;
Yshell = kappa Rshell;(*rescaled dimensionless variables*)
e0 = (2 *sigmainner  Lbjerrum)/(
Rinner kappa);(*inner boundary condition*)
rhoshell = (ShellQ sigmainner)/(kappa^2 Pi (Rshell^2 - Rinner^2))
4 Pi Lbjerrum;(*charge density in surrounding shell. ShellQ=
1 means the shell charge is equal and opposite to the inner \
cylinder*)
func =
NDSolve[{\[Psi]''[r] + 1/r \[Psi]'[r] ==
Sinh[\[Psi][r]] - rhoshell HeavisideTheta[Yshell - r], \[Psi]'[
Youter] == 0, -\[Psi]'[Yinner] == e0}, \[Psi], {r, Yinner,
Youter}, Method -> {"StiffnessSwitching"}];
Plot[Evaluate[{\[Psi][kappa r]} /. func], {r, Rinner, Router},
PlotRange -> All]
],
(*this second module is the linearized version*)
Module[{a, b, c, e0, kappa, Rdna = 1.0, Rh1, rhoH, rhoD, Phi0,
baa = 0.3, nu = 0.6, FperL, sigmaD = sD, csalt, Qh = 43, Lh = 110,
Lbjerrum = 0.7},
kappa = Sqrt[8 Pi csalt Lbjerrum];
e0 = (2 *sigmaD  Lbjerrum)/( Rdna kappa);
Rh1 = shellsize;
rhoD = Rdna*kappa; rhoH = Rh1*kappa;
csalt = 0.16 (6.022*10^23)/0.001 (1/10^9)^3(*convert to nm^-3*);
Phi0 = 1/(2 csalt) Abs[sigmaD]/(Pi (Rh1^2 - Rdna^2));
a = -((-e0 BesselI[1, rhoH] BesselK[0, rhoH] -
e0 BesselI[0, rhoH] BesselK[1, rhoH] +
Phi0 BesselI[1, rhoD] BesselK[1, rhoH])/(BesselK[1,
rhoD] (BesselI[1, rhoH] BesselK[0, rhoH] +
BesselI[0, rhoH] BesselK[1, rhoH])));
b = -((Phi0 BesselK[1, rhoH])/(
BesselI[1, rhoH] BesselK[0, rhoH] +
BesselI[0, rhoH] BesselK[1, rhoH]));
c = -((-e0 BesselI[1, rhoH] BesselK[0, rhoH] -
Phi0 BesselI[1, rhoH] BesselK[1, rhoD] -
e0 BesselI[0, rhoH] BesselK[1, rhoH] +
Phi0 BesselI[1, rhoD] BesselK[1, rhoH])/(BesselK[1,
rhoD] (BesselI[1, rhoH] BesselK[0, rhoH] +
BesselI[0, rhoH] BesselK[1, rhoH])));
Show[Plot[(a BesselK[0, kappa r] + b BesselI[0, kappa r] +
Phi0), {r, Rdna, Rh1}, PlotStyle -> {Red, Dashing[0.01]},
PlotRange -> {{Rdna, Rh1 + 2}, Automatic}],
Plot[(c BesselK[0, kappa r]), {r, Rh1, Rh1 + 2},
PlotStyle -> {Red, Dashing[0.01]}]]
]]]


The fundamental problem here is with Shooting, not Stiffness. A bit of experimentation shows that

func = NDSolve[{\[Psi]''[r] + 1/r \[Psi]'[r] ==
Sinh[\[Psi][r]] - rhoshell HeavisideTheta[Yshell - r], \[Psi]'[
Youter] == 0, -\[Psi]'[Yinner] == e0}, \[Psi], {r, Yinner, Youter},
Method -> {"Shooting", "StartingInitialConditions" -> \[Psi][Yinner] == -2}]


provides the solution for sD = - 7. For values of sD < 3.8, use \[Psi][Yinner] == -1 instead. Note that I have verified this approach only for the nonlinear case, although I expect that it would work well for the linearized case as well.

• Thanks! The linear code is an analytic solution so that will always work, even if it is less accurate. That seems to be important because the Shooting method seems to be very sensitive to the starting initial condition when I push the outer boundary out to 4 or 5. With an outer boundary of 5 a starting condition of 2.5 doesn't work but 2.4 does work. But between the analytic approximation and the shorter boundary I can iterate until I get the right answer. Thanks again! Dec 5, 2022 at 14:07
• @JeremyDS You are correct that Shooting often can be quite sensitive to the initial guess. However, your particular problem is surprising insensitive to the initial guess, although the sensitivity will grow with the domain of integration. There are, of course, alternatives to Shooting, such as Finite Element. By the way, you might consider using the asymptotic symbolic solution as the outer boundary condition, instead of 0. Dec 5, 2022 at 16:06
• It is a nice solution (+1). It is good to know that NDSolve also working for some cases. Dec 11, 2022 at 5:21
• Thanks, @AlexTrounev. I agree that NDSolve could use some improvements. Dec 11, 2022 at 5:27
• @bbgodfrey You are right, they could improve NDSolve. Why not to detect automatically initial condition \[Psi][Yinner] == -2? Dec 11, 2022 at 5:49

We can solve this problem with using LDG method as follows. First, we get

Clear["Global*"]; Get["NumericalDifferentialEquationAnalysis"];


Second step, we run LDG code with parameters

sD = -1.9; shellsize = 1.784; Youter =
kappa Router; par = {sigmainner = sD, csaltM = 0.16, Rinner = 1.,
Router = 3., Rshell = shellsize, ShellQ = -1,
Lbjerrum = 0.701}; csalt =
csaltM (6.022*10^23)/
0.001 (1/10^9)^3(*convert Molar units to nm^-3*);
kappa = Sqrt[
8 Pi csalt Lbjerrum];(*Debye screening parameter in inverse \
nm*)Yinner = kappa Rinner; Youter = kappa Router;
Yshell =
kappa Rshell;(*rescaled dimensionless variables*)e0 = (2*
sigmainner Lbjerrum)/(Rinner kappa);(*inner boundary \
condition*)rhoshell = (ShellQ sigmainner)/(kappa^2 Pi (Rshell^2 -
Rinner^2)) 4 Pi Lbjerrum;(*charge density in surrounding \
shell.ShellQ=1 means the shell charge is equal and opposite to the \
inner cylinder*)

f = {-1/t v[t] + Sinh[p[t]] - rhoshell HeavisideTheta[Yshell - t],
v[t]}; ini = {-e0, b1}; bc1 = 0;
LDGODEs[M0_, nn_, ns_, np_, f_, ini_, tmax_, tmin_] :=
Module[{dx = (tmax - tmin)/nn, A = Array[a, {M0 + 1, nn, 2}]},
xl = Table[tmin + l*dx, {l, 0, nn}]; UT[m_, t_] := EulerE[m, t];
psi[m_, k_, t_] :=
Piecewise[{{UT[
m, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])],
xl[[k]] <= t <= xl[[k + 1]]}, {0, True}}];
g = Table[
GaussianQuadratureWeights[np, xl[[i]], xl[[i + 1]]], {i, nn}];
dp = Table[
D[UT[m, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])],
t], {k, nn}];
rul = {p[t] ->
Sum[a[m, k, 2] psi[m - 1, k, t], {m, 1, M0 + 1}, {k, 1, nn}],
v[t] -> Sum[
a[m, k, 1] psi[m - 1, k, t], {m, 1, M0 + 1}, {k, 1, nn}]};
eq = Flatten[
Table[-Sum[
a[i + 1, n,
ks] (Table[(psi[i, n, t] If[j == 0, 0,
dp[[n]] /. m -> j]), {t, g[[n]][[All, 1]]}] .
g[[n]][[All, 2]]), {i, 0,
M0}] - (Table[(f[[ks]] /. rul) psi[j, n, t], {t,
g[[n]][[All, 1]]}] . g[[n]][[All, 2]]), {j, 0, M0}, {n, 1,
nn}, {ks, 1, 2}] +
Table[Sum[
a[i + 1, n,
ks] (psi[i, n, xl[[n + 1]]] psi[j, n, xl[[n + 1]]]), {i, 0,
M0}] - Sum[
If[n < 2, ini[[ks]]/(M0 + 1),
a[i + 1, n - 1, ks] psi[i, n - 1, xl[[n]]]] psi[j, n,
xl[[n]]], {i, 0, M0}], {j, 0, M0}, {n, 1, nn}, {ks, 2 ns}]];
eq1 = (v[t] - bc1) /. rul /. t -> tmax;
eqn = Join[eq, {eq1}];
var = Join[Flatten[A], {b1}];
sol1 =
FindRoot[Table[eqn[[i]] == 0, {i, Length[eqn]}],
Table[{var[[i]], 1/10}, {i, Length[var]}]]; sol1];
M0 = 4; nn = 10; ns = 1; np = 5; tmin = Yinner; tmax = Youter;
ldgsol = LDGODEs[M0, nn, ns, np, f, ini, tmax, tmin]; // AbsoluteTiming
PSI = Interpolation[
Table[{t,
Evaluate[
Sum[a[i + 1, k, 2] psi[i, k, t], {i, 0, M0}, {k, nn}] /.
sol1]}, {t, tmin, tmax, .01}]];
pn = Plot[PSI[kappa r], {r, Rinner, Router}, PlotRange -> All,
Frame -> True, PlotStyle -> {Green, Dashed}] Plot pn computed with LDG we can compare with NDSolve and linear solution using your module

p1 = Module[{sD = -1.9, shellsize = 1.784},
Show[(*this first module is the nonlinear version*)
Module[{sigmainner = sD, csaltM = 0.16, Rinner = 1., Router = 3.,
Rshell = shellsize, ShellQ = -1, e0, func, Lbjerrum = 0.701,
kappa, csalt, Yinner, Youter, Yshell, rhoshell},
csalt = csaltM (6.022*10^23)/
0.001 (1/10^9)^3(*convert Molar units to nm^-3*);
kappa =
Sqrt[8 Pi csalt Lbjerrum];(*Debye screening parameter in inverse \
nm*)Yinner = kappa Rinner; Youter = kappa Router;
Yshell = kappa Rshell;(*rescaled dimensionless variables*)
e0 = (2*sigmainner Lbjerrum)/(Rinner kappa);(*inner boundary \
condition*)
rhoshell = (ShellQ sigmainner)/(kappa^2 Pi (Rshell^2 -
Rinner^2)) 4 Pi Lbjerrum;(*charge density in surrounding \
shell.ShellQ=
1 means the shell charge is equal and opposite to the inner \
cylinder*)
func = NDSolve[{\[Psi]''[r] + 1/r \[Psi]'[r] ==
Sinh[\[Psi][r]] -
rhoshell HeavisideTheta[Yshell - r], \[Psi]'[Youter] ==
0, -\[Psi]'[Yinner] == e0}, \[Psi], {r, Yinner, Youter},
Method -> {"StiffnessSwitching"}];
Plot[Evaluate[{\[Psi][kappa r]} /. func], {r, Rinner, Router},
PlotRange -> All, PlotStyle -> Blue,
Frame -> True]],(*this second module is the linearized version*)
Module[{a, b, c, e0, kappa, Rdna = 1.0, Rh1, rhoH, rhoD, Phi0,
baa = 0.3, nu = 0.6, FperL, sigmaD = sD, csalt, Qh = 43,
Lh = 110, Lbjerrum = 0.7}, kappa = Sqrt[8 Pi csalt Lbjerrum];
Rh1 = shellsize;
rhoD = Rdna*kappa; rhoH = Rh1*kappa;
csalt = 0.16 (6.022*10^23)/0.001 (1/10^9)^3(*convert to nm^-3*);
Phi0 = 1/(2 csalt) Abs[sigmaD]/(Pi (Rh1^2 - Rdna^2));
a = -((-e0 BesselI[1, rhoH] BesselK[0, rhoH] -
e0 BesselI[0, rhoH] BesselK[1, rhoH] +
Phi0 BesselI[1, rhoD] BesselK[1, rhoH])/(BesselK[1,
rhoD] (BesselI[1, rhoH] BesselK[0, rhoH] +
BesselI[0, rhoH] BesselK[1, rhoH])));
b = -((Phi0 BesselK[1, rhoH])/(BesselI[1, rhoH] BesselK[0, rhoH] +
BesselI[0, rhoH] BesselK[1, rhoH]));
c = -((-e0 BesselI[1, rhoH] BesselK[0, rhoH] -
Phi0 BesselI[1, rhoH] BesselK[1, rhoD] -
e0 BesselI[0, rhoH] BesselK[1, rhoH] +
Phi0 BesselI[1, rhoD] BesselK[1, rhoH])/(BesselK[1,
rhoD] (BesselI[1, rhoH] BesselK[0, rhoH] +
BesselI[0, rhoH] BesselK[1, rhoH])));
Show[
Plot[(a BesselK[0, kappa r] + b BesselI[0, kappa r] + Phi0), {r,
Rdna, Rh1}, PlotStyle -> {Red, Dashing[0.01]},
PlotRange -> {{Rdna, Rh1 + 2}, Automatic}],
Plot[(c BesselK[0, kappa r]), {r, Rh1, Rh1 + 2},
PlotStyle -> {Red, Dashing[0.01]}]]]]];


Show pn and p1 in one plot

Show[p1, pn] We see that in this case LDG solution not so differ from NDSolve. Finally, we put sD = -6, compute LDG solution and compare with linear solution All solutions in one plot look like • Wow. This is great. Thanks! Dec 4, 2022 at 3:19
• @JeremyDS You are welcome! Dec 4, 2022 at 4:33