使用 dopri5 的 Python 绘图答案

【问题标题】：Python plotting using dopri5使用 dopri5 的 Python 绘图
【发布时间】：2016-02-13 05:06:05
【问题描述】：

我的目标是绘制以下一组耦合 ODE：

odeint 方法给我带来了问题（例如某些 x_i 采用分析上不会发生的负值），因此我决定使用四阶 Runge-Kutta 求解器。我是按照这里的例子使用dopri5方法（Using adaptive step sizes with scipy.integrate.ode）：

import scipy as sp
import pylab as plt
import numpy as np
import scipy.integrate as spi

#Constants
c13 = 6.2
c14 = 1.0
c21 = 7.3
c32 = 2.4
c34 = 12.7
c42 = 5.7
c43 = 5.0

e12 = 1.5
e23 = 2.5
e24 = 2.0
e31 = 3.0
e41 = 4.8

#Time
t_end = 700
t_start = 0
t_step = 1
t_interval = sp.arange(t_start, t_end, t_step)

#Initial Condition
ic = [0.2,0.3,0.3,0.5]

def model(t,ic):
    Eqs= np.zeros((4))
    Eqs[0] = (ic[0]*(1-ic[0]*ic[0]-ic[1]*ic[1]-ic[2]*ic[2]-ic[3]*ic[3])-c21*((ic[1]*ic[1])*ic[0])+e31*((ic[2]*ic[2])*ic[0])+e41*((ic[3]*ic[3])*ic[0]))
    Eqs[1] = (ic[1]*(1-ic[0]*ic[0]-ic[1]*ic[1]-ic[2]*ic[2]-ic[3]*ic[3])+e12*((ic[0]*ic[0])*ic[1])-c32*((ic[2]*ic[2])*ic[1])-c42*((ic[3]*ic[3])*ic[1]))
    Eqs[2] = (ic[2]*(1-ic[0]*ic[0]-ic[1]*ic[1]-ic[2]*ic[2]-ic[3]*ic[3])-c13*((ic[0]*ic[0])*ic[2])+e23*((ic[1]*ic[1])*ic[2])-c43*((ic[3]*ic[3])*ic[2]))
    Eqs[3] = (ic[3]*(1-ic[0]*ic[0]-ic[1]*ic[1]-ic[2]*ic[2]-ic[3]*ic[3])-c14*((ic[0]*ic[0])*ic[3])+e24*((ic[1]*ic[1])*ic[3])-c34*((ic[2]*ic[2])*ic[3]))
    return Eqs

ode =  spi.ode(model)

ode.set_integrator('dopri5')
ode.set_initial_value(ic,t_start)
ts = []
ys = []

while ode.successful() and ode.t < t_end:
    ode.integrate(ode.t + t_step)
    ts.append(ode.t)
    ys.append(ode.y)

t = np.vstack(ts)
x1,x2,x3,x4 = np.vstack(ys).T

plt.subplot(1, 1, 1)
plt.plot(t, x1, 'r', label = 'x1')
plt.plot(t, x2, 'b', label = 'x2')
plt.plot(t, x3, 'g', label = 'x3')
plt.plot(t, x4, 'purple', label = 'x4')
plt.xlim([0,t_end])
plt.legend()
plt.ylim([-0.2,1.3])

plt.show()

产生以下情节：

但是，我的图有些奇怪——x1 值似乎随机飙升到 1 以上。因为 (1,0,0,0) 是动态系统的一个不动点，所以它没有多大意义对我来说它高于一个（如果尖峰显示出某种重复模式可能是有原因的，但它们在图中看起来是随机的，所以我想知道这是否与数值积分有更多关系而不是实际的动态）。改变参数也会改变那个尖峰（我试过的一些参数值仍然有那个尖峰，但它要小得多）。因此，我有两个问题：

1) 是什么导致这些尖峰出现？我最初的想法是“t_step”很大，所以我试着玩弄它（这引出了我的第二个问题）

2) 我尝试将步长从 1 更改为 0.1。然而，它产生了这个图，实际上看起来与步长为 1 时有很大不同（我认为较小的步长会产生更准确的图？）在这个图中，似乎每个 x_i 在 1 附近花费的时间比如果步长是 1 而不是 0.1，并且尖峰的高度都是相等的——考虑到这些图有一些实质性差异，我怎么知道哪个图更准确？

我正在研究的问题的真正症结是查找 x3 和 x4 之间的交互，但我想确保模拟设置正确，以便我可以获得关于 x3 和 x4 的准确结果！

【问题讨论】：

标签： python plot ode

【解决方案1】：

根据 Gronwall 引理，“nice”ODE 的误差传播由以 Lipschitz 常数 L 为因子的指数控制。这意味着在时间 t 发生的任何本地错误贡献在时间 T 时（可能）被放大了 exp(L(T-t)) 倍。

慷慨的人可以为您的系统使用L=10。在10 的时间间隔内，这给出了exp(100)=2.688e+43 的误差放大率，或数值与初始值和精确解的完全解耦。

因此，它更令人惊讶更令人惊讶，直到t=80，保留了两种解决方案的相似性。

关于理论的注释，尖峰：没有什么限制变量保持在区间 [0,1] 内或具有一些具有相同效果的守恒量。您可以做的最好的事情是使用像V(x)=sum(x*x) 这样的 Ljapunov 函数，然后给出将解限制在两个球体之间的边界。

固定点是所有双曲线，雅各比亚都有级别缺乏，无论不包括稳定的固定点。然后，尖峰只是这些双曲点的稳定和不稳定的歧管/轴的方向的结果。

【讨论】：