用 SciPy 数值求解 ODE
Posted
技术标签:
【中文标题】用 SciPy 数值求解 ODE【英文标题】:Numerically Solving ODE with SciPy 【发布时间】:2016-01-24 19:38:24 【问题描述】:我坚持将scipy.integrate.odeint
应用于以下非常简单的 ODE:
y(t)/dt = y(t) + t^2 and y(0) = 0
SciPy 计算的解决方案不正确(很可能我在这里混淆了某些东西)——尤其是解决方案不满足初始条件。
import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt
import math
# the definition of the ODE equation
def f(y,t):
return [t**2 + y[0]]
# computing the solution
ts = np.linspace(-3,3,1000)
res = scipy.integrate.odeint(f, [0], ts)
# the solution computed by WolframAlpha [1]
def y(t):
return -t**2 - 2*t + 2*math.exp(t) - 2
fig = plt.figure(1, figsize=(8,8))
ax1 = fig.add_subplot(211)
ax1.plot(ts, res[:,0])
ax1.text(0.5, 0.95,'SciPy solution', ha='center', va='top',
transform = ax1.transAxes)
ax1 = fig.add_subplot(212)
ax1.plot(ts, np.vectorize(y)(ts))
ax1.text(0.5, 0.95,'WolframAlpha solution', ha='center', va='top',
transform = ax1.transAxes)
plt.show()
1 : WolframAlpha: "solve dy(t)/dt = t^2 + y(t), y(0) = 0"
我的错误在哪里?
【问题讨论】:
【参考方案1】:您的 scipy 代码求解了初始条件为 y(-3) = 0
的微分方程,而不是 y(0) = 0
。 odeint
的y0
参数是t
参数第一次给出的值。
在 y(0) = 0 的区间 [-3, 3] 上解决此问题的一种方法是调用 odeint
两次,如下所示:
In [81]: from scipy.integrate import odeint
In [82]: def f(y,t):
....: return [t**2 + y[0]]
....:
In [83]: tneg = np.linspace(0, -3, 500)
In [84]: tpos = np.linspace(0, 3, 500)
In [85]: sol_neg = odeint(f, [0], tneg)
In [86]: sol_pos = odeint(f, [0], tpos)
In [87]: plot(tneg, sol_neg)
Out[87]: [<matplotlib.lines.Line2D at 0x10f890d90>]
In [88]: plot(tpos, sol_pos)
Out[88]: [<matplotlib.lines.Line2D at 0x107a43cd0>]
In [89]: grid(True)
创建
【讨论】:
解决 y(0)=0 但仍计算 [-3,3] 的最简单方法是什么? 可能从 0 计算到 3,然后从 0 计算到 -3以上是关于用 SciPy 数值求解 ODE的主要内容,如果未能解决你的问题,请参考以下文章
如何使用 scipy.integrate.odeint 求解具有时间相关变量的 ODE 系统