python - 在 python 中使用 scipy.integrate.odeint 时遇到问题

Question

我正在尝试使用 odeint 来解决问题。我的代码如下:

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

eta=1.24e-9/2
def fun(x):
    f=1.05e-8*eta*x**(1.5)*np.exp(13.6/x)
    return (np.sqrt(1.+4*f)-1)/2./f
x=np.arange(0,1,0.001)
y=odeint(fun,x,0)[0]
plt.plot(x,y)
plt.plot(x,x)
plt.show()

如果两条曲线相同，显然是错误的。如果我绘制该函数，它将看起来像一个阶梯函数，它在大约 0.3 之前非常非常小，然后呈指数上升到 1。你能帮我找出它有什么问题吗？谢谢!

Answer 1

您的代码有几个问题，如果您阅读 docstring for odeint，您可能可以自己解决其中的大部分问题更仔细。

为了帮助您入门，以下是使用 odeint 求解标量微分方程的简单示例。我不会尝试理解(并可能调试)您的函数，而是使用一个非常简单的等式。我将求解方程 dy/dt = a * y，初始条件为 y(0) = 100。一旦您使这个示例正常运行，您可以修改 fun 来解决您的问题。

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


def fun(y, t, a):
    """Define the right-hand side of equation dy/dt = a*y""" 
    f = a * y
    return f


# Initial condition
y0 = 100.0

# Times at which the solution is to be computed.
t = np.linspace(0, 1, 51)

# Parameter value to use in `fun`.
a = -2.5

# Solve the equation.
y = odeint(fun, y0, t, args=(a,))

# Plot the solution.  `odeint` is generally used to solve a system
# of equations, so it returns an array with shape (len(t), len(y0)).
# In this case, len(y0) is 1, so y[:,0] gives us the solution.
plt.plot(t, y[:,0])
plt.xlabel('t')
plt.ylabel('y')
plt.show()

这是情节:

可以在 SciPy Cookbook 中找到使用 odeint 的更复杂的示例(向下滚动到标有“常微分方程”的项目符号)。