0

I am relatively new to Python and trying to use it to solve a second order nonlinear differential equation, specifically the Poisson-Boltzmann equation in an electrolyte.

phi''(r) + (2/r)*phi'(r) = (k^2)*sinh(phi(r))

Essentially it describes the decay of electrostatic potential (phi) away from a charged surface in an electrolyte with the rate of decay governed by a paramter k.

  • phi(r) - the potential at r
  • dphi(r) - derivative of potential at r
  • r - distance from the surface (I am solving for r = 1 to r = R in this case)

and the boundary conditions

  • phi(1) = 5
  • dphi(R) = 0

The problem bit of code is as follows

from scipy.integrate import odeint
from scipy.optimize import root
from pylab import * # for plotting commands

k = 0.5 
phi = 5 
dphi = -10
R = 21

def deriv(z,r): # return derivatives of the array z   (where z = [phi, phi'])
    return np.array([
    (z[1]),
    ((k**2)*sinh(z[0]))-((2/r)*z[1])
    ])


result = odeint(deriv,np.array([phi,dphi]),np.linspace(1,R,1017), full_output = 1)

Generally for low values of k the integration works fine and I can use root from scipy.optimize to solve it subject to the boundary conditions, however when I try to use relatively large values of k (0.5 and higher) the integration runs into problems and outputs the following error

Excess work done on this call (perhaps wrong Dfun type).

Having run it with full_output = 1 and had a look at the tcur parameter it seems to count up smoothly until a certain point and then oscillate wildly from very large to very small values. At the same point nfe the number of function evaluations drops to zero. When working correctly the tcur parameter runs smoothly from 1 to R. Could anyone enlighten me to why this is happening? Is it a problem with my implementation or is there an instability in the equation?

Thanks in advance for any help,

Dave

David Gillespie
  • 1
  • 1

1 Answers1

0

The ODE is probably unstable. The related equation

phi''(r) = (k^2)*sinh(phi(r))

has a solution for k=1 (for corresponding initial conditions at r=1):

phi(r) = 2 arctanh(sin(r))

The solution has a singularity at r=pi/2. A numerical ODE solver will not be able to integrate the equation beyond this point. It's plausible that a similar equation with the first-derivative term (which should be negligible close to singularities anyway) behaves similarly.

The actual problem that you have is that a shooting method using an ODE solver is not a good way to solve boundary value problems --- you should use collocation methods, which are fairly stable. See e.g. http://www.scholarpedia.org/article/Boundary_value_problem and references therein.

For Python software, see https://pypi.python.org/pypi?%3Aaction=search&term=boundary+value+problem&submit=search

It's usually also very easy to write such solvers yourself, as the only needed step is discretization of the problem to a set of (many) equations, after which root can solve them.

pv.
  • 33,875
  • 8
  • 55
  • 49