使用 scipy.integrate.ode 求解 Ode 系统

Question

我正在尝试在 python 中解决 Du/Dt = F(u) 形式的 Odes 系统，我怀疑我可能在某处犯了一个相当愚蠢的错误。

从技术上讲，F(u) 实际上是 u 关于另一个变量 y 的二阶导数，但实际上我们可以将其视为一个系统和某个函数。

#Settings#
minx = -20
h = float(1) 
w = float(10)
U0 = float(10)
Nt = 10
Ny = 10
tmax = 10

v=float(1)

#Arrays#
y = np.linspace(0,h,Ny)
t = np.linspace(0,tmax,Nt)

#Variables from arrays#
dt = t[1]-t[0]
p = [0]*(Nt)
delta = y[1] - y[0]

def zo(y):
    return math.cos(y/(2*math.pi))

z0 = [zo(i) for i in y]

def df(t,v1):

    output = np.zeros(len(y))

    it = 1

    output[0] = math.cos(w*t)

    output[len(y)-1] = math.cos(w*t)

    while it < len(y)-1:

        output[it] = ( v1[it - 1] + v1[it + 1] - 2 * v1[it] ) * ( v / ( ( delta )**2 ))

        it += 1

    return output

r = ode(df).set_integrator('zvode', method='bdf',order =15)

r.set_initial_value(z0, 0)

it=0

while r.successful() and r.t < tmax:

    p[it] = r.integrate(r.t+dt)

    it+=1

print(z0-p[0])
print(p[1])

现在问题是双重的：

-首先，初始的"condition"即p[0]似乎是关闭的。（这可能只是因为 ode 函数的工作方式，所以我不知道这是否正常）

-其次，p[1]和之后的所有p都是0。

因此出于某种原因，ode 函数立即失败...（您可以通过在初始化 p 时将值更改为 1 来检查这一点）

除了我知道这个方法应该有效。毕竟这是 matlab 中 ode45 的 "equivalent" 并且绝对有效。

Answer 1

如果您想使用 Dormand-Price rk45 resp，为什么要选择具有相当高阶的隐式后向微分公式的复杂求解器。 dopri5？

还请更正 df 中的循环缩进。为什么不在 range(1, len(y)-1) 上进行 for 循环？

目前 p[0] 包含第一步之后的解决方案点，位于 t=1*dt。您必须显式分配 p[0]=z0 并启动 it=1 才能在 p 中获得完整的解决方案路径。查看p的长度，可能是你需要Nt+1.

使用 scipy.integrate.ode 求解 Ode 系统

Solving a System of Ode's using scipy.integrate.ode

python

system

scipy

ode