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）

和边界条件

  • φ(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 平稳运行。有人能告诉我为什么会这样吗？是我的实现有问题还是等式不稳定？

提前感谢您的帮助，

戴夫

4

1 回答 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可以求解它们。

于 2014-12-14T13:58:13.407 回答