使用 SciPy 求解 ODE 会在 double_scalars 错误中遇到无效值 [英] Solving ODE with SciPy gives invalid value encountered in double_scalars error

问题描述

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

将 matplotlib.pyplot 导入为 plt将 numpy 导入为 np从 scipy.integrate 导入 o​​deintdef poly(y,xi,n):θ, phi = ydydt = [phi/(xi**2), -(xi**2)*theta**n]返回 dydt

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

y0 = [1.,0.]xi = np.linspace(0., 16., 201)对于范围内的 n(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: 在 double_scalars 中遇到无效值
app.launch_new_instance()

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

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

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

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

解决方案

以下对我有用(看起来类似于

<小时>

最初的问题使用了多方指数的分数值，因此可以使用诸如 follow 之类的东西来调查这些值......

for n in range(12):sol = odeint(poly, y0, xi, args=(n/2.,))如果 np.isnan(sol[:,1]).any():塞子 = np.where(np.isnan(sol[:,1]))[0][0]ax.plot(xi[:stopper], sol[:stopper, 0])别的:ax.plot(xi, sol[:, 0])

I'm trying to solve the Lane-Emden equation for arbitrary values of n (polytropic index). In order to use SciPy, I have represented the second-order ODE as a set of two coupled first-order ODEs. I have the following code:

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

Here, I defined 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()

However, running this code results in the following error:

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

At first, I thought this was a result of the singularity of the equation at xi = 0, so I changed my integration interval:

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

This fixes the problem at n = 0., but not at other values of n. The plots I get don't make any sense and are just plain wrong.

Why do I get this error message, and how can I fix my code?

解决方案

The following works for me (which looks similar to the Wikipedia entry). The derivative was written incorrectly (negative on the wrong element) and input to the derivative is changed simply to 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()

---

The original question used fractional values of the polytropic index, so something like follow could be used to investigate those values ...

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])

这篇关于使用 SciPy 求解 ODE 会在 double_scalars 错误中遇到无效值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！