求解带迭代参数的微分方程

Question

我想使用（scipy.integrate.odeint）和（scipy.integrate.odeint）解微分方程。我有一个简单的例子：

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.我是否使用ode错误？如何减少错误？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.

我想使用（scipy.integrate.odeint）和（scipy.integrate.odeint）解微分方程。我有一个简单的示例：dy / dt = f [i] * t，f是与t [i]对应的参数，就像...

Answer 1

您正在做的事情更接近于积分y'（t）= t ^ 2，y（0）= 0，导致y（t）= t ^ 3/3。将t ^ 2视为f * t并将f分解为t的阶跃函数版本只会给它带来一点微扰。