CODEWITHSUNDEEP
scipydifferential-equationsnonlinear-functionspoissonodeint
David Gillespie
David Gillespie

Reputation: 1

Having problems with ODEINT in python

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.

and the boundary conditions

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

Upvotes: 0

Views: 1710

Answers (2)

Tarik
Tarik

Reputation: 11209

I modified the code to use solve_ivp instead of odeint. I selected Radau as the solving method and generated the solutions for various values of k.

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

ks = [0.02, 0.05, 0.1, 0.5, 0.75, 1, 2, 3, 4]
phi = 5 
dphi = -10
R = 21

r0 = 1 
rf = R
r_span = np.array([r0, rf])
r_eval = t_eval=np.linspace(1,R,1017)

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


results = [solve_ivp(deriv, r_span, np.array([phi,dphi]), t_eval=r_eval, method='Radau', vectorized=True, args=(k,)) for k in ks]

for k, result in zip(ks, results):
    plot(result.t, result.y[0], label=f'k = {k}')
    
plt.legend()
plt.grid()

Plot of solutions: Plot of Solutions

As seen in the plots, the solution either shoot up or down to + or - infinity

Upvotes: 0

pv.
pv.

Reputation: 35115

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.

Upvotes: 0

Related Questions