用迭代参数求解微分方程

Question

我正在学习使用 (scipy.integrate.odeint) 和 (scipy.integrate.ode) 求解微分方程。我有一个简单的例子：

dy/dt=f[i]*t

f是t[i]对应的参数，同代码中的例子；即

t[0]=0.0, f[0]=0.0

t[1]=0.1, f[1]=0.1

...

t[10]=1.0, f[1]=1.0

手动结果应该是：

y=1/2*f[i]*t**2，因为y的初始值为零

那么y的数值结果应该是 [0.0, 0.0005, 0.004, 0.0135, 0.032, 0.0625, 0.108, 0.1715, 0.256, 0.3645, 0.5]。但是当我使用 scipy.integrate.ode 时，我得到了不同的结果。我的问题是： 1.我用错了吗？我怎样才能减少错误？ 2.我可以用odeint或者其他方法解决这个问题吗？

代码如下：

import matplotlib.pyplot as pl
import numpy as np
import sympy as sp
from scipy.integrate import odeint
from scipy.integrate import ode
import numpy as np

def func(t, y, f):
    return f*t

t=np.linspace(0.0, 1.0, 11)
f=np.linspace(0.0, 1.0, 11)
dt = t[1]-t[0]

sol= np.empty_like(t)
r = ode(func).set_integrator("dopri5")
r.set_initial_value(0, 0).set_f_params(f[0])

# result of ode
for i in xrange(len(t)):
    r.set_f_params(f[i])
    r.integrate(r.t+dt)
    sol[i] = r.y

res=[]
# result of t**3/3
for a in np.linspace(0.0, 1, 11):
    f=(a**3)/3
    print f
    res.append(f)

# result3
res2=[]
for n in range(0, 11):
    dt=0.1
    y= t[n]**3/3 - dt*t[n]**2/4 - dt**2*t[n]/12
    res2.append(y)

pl.plot(sol)
pl.plot(res)  
pl.plot(res2)
pl.show()

我已将此示例扩展到二维微分方程：

du/dt=-u(v-f[i])

dv/dt=v(f[i]-u)

初始值：u(0)=v(0)=1。下面是代码：

import matplotlib.pyplot as pl
import numpy as np
import sympy as sp
from scipy.integrate import odeint
from scipy.integrate import ode
from numpy import array

def func(t, y, f):
    u,v=y
    dotu=-u*(v-f)
    dotv=v*(f-u)
    return array([dotu, dotv])

t=np.linspace(0.0, 10, 11)
f=np.linspace(0.0, 20, 11)

dt = t[1]-t[0]

# result correct
y0=array([1.0, 1.0])
sol= np.empty([11, 2])

sol[0] = array([1.0, 1.0])
r = ode(func).set_integrator("dopri5")
r.set_initial_value(t[0], sol[0]).set_f_params(f[0])

for i in range(len(t)-1):
    r.set_f_params(f[i])
    r.integrate(r.t+dt)
    sol[i+1] = r.y

pl.plot(sol[:,0])

但我收到一条错误消息：

Traceback (most recent call last): File "C:\Users\odeint test.py", line 26, in <module> sol[0] = array([1.0, 1.0]) ValueError: setting an array element with a sequence.

Answer 1

你所做的更接近于积分y'(t)=t^2, y(0)=0, 导致y(t)=t^3/3。将 t^2 分解为 f*t 并将 f 分解为 t 的阶跃函数版本只会增加一个小的扰动。

t[i]*t对t[i]..t[i+1]的积分是

y[i+1]-y[i] = t[i]/2*(t[i+1]^2-t[i]^2) 
= (t[i+1]^3-t[i]^3)/3 - (t[i+1]-t[i])^2*(t[i]+2t[i+1])/6
= (t[i+1]^3-t[i]^3)/3 - dt*(t[i+1]^2-t[i]^2)/4 - dt^2*(t[i+1]-t[i])/12

总计约为

y[n] = t[n]^3/3 - dt*t[n]^2/4 - dt^2*t[n]/12

如何获得正确的解决方案

sol= np.empty_like(t)

设置初始值

sol[0] = 0
r = ode(func).set_integrator("dopri5")

使用初始点作为初始点，既明确索引0处的点是固定的，又"used up"

r.set_initial_value(sol[0],t[0]).set_f_params(f[0])

# result of ode

从 t[i] 点到 t[i+1] 点。以 i+1=len(t) 或 i=len(t)-1

结尾

for i in xrange(len(t)-1):
    r.set_f_params(f[i])
    r.integrate(r.t+dt)

t[i]+dt 的值是 t[i+1]

的值

    sol[i+1] = r.y

通过这些更改，数值解与手动计算的解一致。

用迭代参数求解微分方程

solve differential equation with iterated parameter

python

ode

如何获得正确的解决方案