python - 使用 SciPy 求解 ODE 给出在 double_scalars 错误中遇到的无效值

Question

我正在尝试针对 n(多变指数)的任意值求解 ​​Lane-Emden 方程。为了使用 SciPy，我将二阶 ODE 表示为一组两个耦合的一阶 ODE。我有以下代码:

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

def poly(y,xi,n):  
    theta, phi = y  
    dydt = [phi/(xi**2), -(xi**2)*theta**n]  
    return dydt

在这里，我定义了 phi = theta'

y0 = [1.,0.]  
xi = np.linspace(0., 16., 201)    
for n in range(0,11):  
    sol = odeint(poly, y0, xi, args=(n/2.,))  
    plt.plot(xi, sol[:, 0], label=str(n/2.))  
plt.legend(loc='best')  
plt.xlabel('xi')  
plt.grid()  
plt.show()

但是，运行此代码会导致以下错误:

RuntimeWarning: invalid value encountered in double_scalars
app.launch_new_instance()

起初，我认为这是方程在 xi = 0 处的奇异性的结果，所以我更改了积分区间:

xi = np.linspace(1e-10, 16., 201)

这解决了 n = 0 时的问题，但不是 n 的其他值。我得到的情节没有任何意义，而且完全错误。

为什么我会收到此错误消息，我该如何修复我的代码？

Answer 1

以下对我有用(看起来类似于 Wikipedia entry )。导数写错了(负数在错误的元素上)，导数的输入被简单地更改为 n。

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

def poly(y,xi,n):  
    theta, phi = y  
    dydxi = [-phi/(xi**2), (xi**2)*theta**n]  
    return dydxi

fig, ax = plt.subplots()

y0 = [1.,0.]  
xi = np.linspace(10e-4, 16., 201)

for n in range(6):
    sol = odeint(poly, y0, xi, args=(n,))
    ax.plot(xi, sol[:, 0])

ax.axhline(y = 0.0, linestyle = '--', color = 'k')
ax.set_xlim([0, 10])
ax.set_ylim([-0.5, 1.0])
ax.set_xlabel(r"$\xi$")
ax.set_ylabel(r"$\theta$")
plt.show()

---

最初的问题使用了多变指数的分数值，所以可以使用类似下面的内容来研究这些值......

for n in range(12):
    sol = odeint(poly, y0, xi, args=(n/2.,))

    if np.isnan(sol[:,1]).any():
        stopper = np.where(np.isnan(sol[:,1]))[0][0]
        ax.plot(xi[:stopper], sol[:stopper, 0])
    else:
        ax.plot(xi, sol[:, 0])