首页 文章

在python中遇到ODEINT问题

提问于
浏览 次
0

我是Python的新手，并试图用它来解决二阶非线性微分方程，特别是电解质中的Poisson-Boltzmann方程 .

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

基本上它描述了静电势（phi）远离电解质中的带电表面的衰减，其中衰减速率由参数k控制 .

  • phi（r） - r的潜力

  • dphi（r） - r的潜在导数

  • r - 距离曲面的距离（在这种情况下我求解r = 1到r = R）

和边界条件

  • phi（1）= 5

  • dphi（R）= 0

问题的代码如下

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)

通常对于低k值，积分工作正常，我可以使用scipy.optimize中的root来根据边界条件求解它，但是当我尝试使用相对较大的k值（0.5和更高）时，积分会遇到问题，输出以下错误

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

运行full_output = 1并查看 tcur 参数后，它似乎顺利计数直到某一点，然后从非常大的值大幅度振荡到非常小的值 . 在同一点 nfe ，功能评估的数量下降到零 . 正确工作时，tcur参数从1到R顺利运行 . 任何人都可以告诉我为什么会发生这种情况？这是我的实施问题还是方程中存在不稳定性？

在此先感谢您的帮助，

戴夫

scipy differential-equations nonlinear-functions poisson odeint

1 回答

  • 0

    ODE可能不稳定 . 相关方程式

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

    有一个 k=1 的解决方案（对于r = 1时的相应初始条件）：

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

    该解决方案在 r=pi/2 处具有奇点 . 数值ODE求解器将无法将该等式积分超出此点 . 有一个类似的方程与一阶导数项（无论如何应该可以忽略不计，接近奇点）的行为似乎是相似的 .

    你遇到的实际问题是使用ODE求解器的射击方法不是解决边界值问题的好方法 - 你应该使用相当稳定的搭配方法 . 参见例如http://www.scholarpedia.org/article/Boundary_value_problem及其中的参考文献 .

    对于Python软件，请参阅https://pypi.python.org/pypi?%3Aaction=search&term=boundary+value+problem&submit=search

    通常也很容易自己编写这样的求解器，因为唯一需要的步骤是将问题离散化为一组（多个）方程，之后 root 可以解决它们 .

    回复于

相关问题