I'm trying to solve a tricky 2nd order ODE (the 1D Poisson-Boltzmann equation) with NDSolve for the electric potential at a distance (r) from the center of a charged particle of radius (a) at a density (d) with charge (q); it has two Neumann boundary conditions (BCs). It's an intricate equation and it isn't surprising it's hard [for me] to get a numerical solution. I feel I am close, but as yet no cigar...

soln =
 With[{a = 56/2*10^-9, d = 2.22814*10^19, q = 900},
{f''[r] + (2 f'[r])/r == (1.03975*10^8)^2 Sinh[f[r]],
 f'[a] == -((q (.714295*10^-9))/a^2),
 f'[a*((4. π d a^3)/3.)^(-1/3)] == 0}, f,
{r, 2*10^-8, 2.5*10^-7}]] // Flatten

FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations.>>

NDSolve::berr: The scaled boundary value residual error of 8.240798196770146`*^9 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.>>

I feel the two are likely related - if it had more iterations it may converge to a solution that satisfies the boundary values. This explains the first of the following attempts at resolving the issue: (listed in the order I suspect might be most valuable)

  • Instructing it to run more iterations, but MaxSteps doesn't do the job and MaxIterations is not supported as an option of NDSolve; I feel inept not finding how to specify this.

  • I've read that using Method -> {"Shooting","StartingInitialConditions" -> {f'[a] ==1.25((q (.714295*10^-9))/a^2)}} i.e. simply adding a coefficient to start slightly (and I've tried not so slightly) away from the initial condition might work. This returns singularity / stiffness errors. This question seems to indicate that if I choose the starting IC more wisely, it may eliminate the stiffness issues, but I'm at a loss as to what else to try.

  • Changing the coefficient of the RHS of the equation (the value 1.03975*10^8) corresponding to a parameter that will eventually be fed into ParametricNDSolve. This leads to stiffness issues that can be resolved with Method -> "StiffnessSwitching, only to return the same errors quoted above.

  • Distances (a), (r), that in the density (d), etc. are in meters; I've similar problems in the more natural length scale of nanometers.

  • Lowering the AccuracyGoal and PrecisionGoal. Simply increases the discrepancy with the BC.

  • Adding understood BC at infinity (f[inf] = 0). Mathematica crashes, understandably.

Does anyone have any idea how to resolve the issue? For convenience, the function to plot the result of the integration over the concerning range (with a gridlines at each BC point, where the green ones are the value of the function f[r] itself that I really aim to acquire):

Module[{a = 56/2*10^-9, d = 2.22814*10^19},
 Plot[{f[r]} /. soln, {r, 2*10^-8, 2.5*10^-7}, PlotStyle -> Thick,
  GridLines -> {{{a, Red}, {2.20457*10^-7, Darker@Green}},  
   {{f[a*((4 π d a^3)/3.)^(-1/3)] /. soln, Darker@Green}}}]]

1 Answer 1


A step backward i.e. turn to finite difference method (FDM) seems not to be a bad choice here. I'll use pdetoae for the generation of difference equations:

With[{a = 56/(2 10^9), d = 2.22814 10^19, q = 900}, 
 domain = Rationalize[#, 0] &@{a, a ((4. π d a^3)/3.)^(-1/3)};
 eq = f''[r] + (2 Derivative[1][f][r])/r == (1.03975 10^8)^2 Sinh[f[r]];
 bc = Rationalize[#, 0] &@{f'[a] == -((q .714295)/(10^9 a^2)), 
    f'[a ((4. π d a^3)/3.)^(-1/3)] == 0};]

points = 100;
grid = Array[# &, points, domain];
difforder = 2;
(* Definition of pdetoae isn't included in this code piece,
   please find it in the link above. *)
ptoafunc = pdetoae[f[r], grid, difforder];
ae = ptoafunc[eq][[2 ;; -2]];
aebc = ptoafunc[bc];
soldata = FindRoot[Rationalize[{ae, aebc} // Flatten, 0], {f@#, 1/2} & /@ grid, 
    WorkingPrecision -> 16][[All, -1]];
ListLinePlot[soldata, DataRange -> domain, PlotRange -> All]

Mathematica graphics

  • $\begingroup$ Thank you, xzczd. This seems to work and is very impressive, but in my inexperience with FDM and your code, I can't fully tell why. I am going to leave the question open for awhile to see if someone can figure out a solution with traditional NDSolve, but if no one does I will certainly accept your answer. While I suppose I still won't understand the internals if they do solve it with NDSolve, it would be good for future users' and my peace of mind. That's no slight to the internals of your code; I should learn FDM instead of relying on the almighty WRI. Thank you again. :-) $\endgroup$
    – Ghersic
    Oct 9, 2017 at 12:11
  • $\begingroup$ @Ghersic Feel free to wait for a better answer. :) $\endgroup$
    – xzczd
    Oct 9, 2017 at 12:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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