在 python 中使用 ODEINT 时遇到问题

Question

我对 Python 比较陌生，并尝试使用它来求解二阶非线性微分方程，特别是电解质中的泊松-玻尔兹曼方程。

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。有人能告诉我为什么会发生这种情况吗？是我的实现有问题还是方程不稳定？

预先感谢您的帮助，

戴夫

Answer 1

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

可以求解它们。

Answer 2

我修改了代码以使用

solve_ivp

而不是

odeint

。我选择

Radau

作为求解方法并生成不同 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 = 1017
r_span = np.array([r0, R])
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()

解图：