Scipy.integrate solve_ivp 的稳定性问题

Question

我正在使用 SciPy 来求解微分方程组，但我在解的稳定性方面遇到了一些问题。本质上，有两个状态变量，它们应该走向两个相对稳定的解，一个渐进趋近于零（本质上为零），一个趋近于0.59。

当我在 R 中使用 deSolve 包执行此操作时，它工作正常，但我正在尝试使用 scipy 并且我得到了一个奇怪的东西，最终的解决方案，状态变量不稳定即使他们应该是？应该为零的状态变量的值从 7.41e-323 顺序的一系列值开始，然后向上跳到 5.3e-001 一步，然后返回。我不知道这是为什么，但我想知道是否有人可以就如何 a) 解决此问题或 b) 使用 scipy?

以外的其他选项提供任何建议

在 R（lsoda 包）和 Maple 中尝试这个已经产生了稳定的解决方案。

谢谢！

代码：

# %%
import numpy as np
import scipy.integrate as spi
from scipy.integrate import solve_ivp

# %%
x_coords = np.arange(0.01,1,0.05)
#this makes a mesh grid at those points (the whole square)
M_temp, C_temp = np.meshgrid(x_coords, x_coords)

#only includes the points I want (the lower triangular in M x C space)
C_points = C_temp[M_temp + C_temp <=1]
M_points = M_temp[M_temp + C_temp <=1]
T_points = 1 - C_points - M_points
# initialize array
M_array = np.zeros((50000,len(C_points)))
C_array = np.zeros((50000,len(C_points)))

# %%
for i in range(0,len(C_points)):
    def rhs(s,v):
        a = 0.25
        g = 0.3
        z = 0.05
        r = 0.55 
        d = 0.24 
        y = 0.77 
        return [r*v[0]*v[2] + z*r*v[2] - d*v[0] - a*v[0]*v[1], a*v[0]*v[1] -
        (g*v[1])/(v[1]+v[2]) + y*v[1]*v[2],-r*v[0]*v[2] - z*r*v[2] + d*v[0]+
        (g*v[1])/(v[1]+v[2]) - y*v[1]*v[2]]

    res = solve_ivp(rhs, (0, 50000),
        [C_points[i], M_points[i], T_points[i]],
        t_eval =np.arange(0,50000,1)) #solves from t =0 -> t = 50000
    M_array[:,i] = res.y.T[:,1]
    C_array[:,i] = res.y.T[:,0] #[0:50,2]

trajectories = np.ones(len(C_points)*len(M_array[:,i]),\
    dtype = {'names': ['M', 'C'], 'formats':[np.float64, np.float64]})
trajectories['M'] = M_array.flatten()
trajectories['C'] = C_array.flatten()

print(trajectories['M'][1049900:1049999])

输出：

Answer 1

如果您打印 C-M 网格中最后一个 M 分量的点后的第一个小数点，您会得到

0.01: 0 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5  
0.06: 0 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5  
0.11: 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5  
0.16: 0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5 5  
0.21: 0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5  
0.26: 0 0 0 0 5 5 5 5 5 5 5 5 5 5 5  
0.31: 0 0 0 0 0 5 5 5 5 5 5 5 5 5  
0.36: 0 0 0 0 0 0 5 5 5 5 5 5 5  
0.41: 0 0 0 0 0 0 0 5 5 5 5 5  
0.46: 0 0 0 0 0 0 0 0 5 5 5  
0.51: 0 0 0 0 0 0 0 0 0 5  
0.56: 0 0 0 0 0 0 0 0 0  
0.61: 0 0 0 0 0 0 0 0  
0.66: 0 0 0 0 0 0 0  
0.71: 0 0 0 0 0 0  
0.76: 0 0 0 0 0  
0.81: 0 0 0 0  
0.86: 0 0 0  
0.91: 0 0  
0.96: 0

这个模式强烈表明观察到的模式不是随机误差，而是系统的一个属性，一组初始值，大致为C0 < M0，导致一个不动点具有非零 M 分量，最佳估计值 0.53371702。

前两个分量的流图证实了这一假设，[0.31024394, 0.22563217] 处的鞍点将汇聚到稳定点 [0.07006604, 0.53371702] 和 [0.59734069, 0.0]

的区域分开

再次检查 R 和 python 版本中的代码确实是一样的。找出R中的积分方式，大概是lsoda。然后找出默认公差，或设置合理的公差 (atol = 1e-11, rtol=1e-13)，然后在 python 中复制这些公差。使用method(="BDF")选项设置solve_ivp.

的积分方式

Scipy.integrate solve_ivp 的稳定性问题

Problems with stability of Scipy.integrate solve_ivp

python

numpy

scipy

ode

代码：

输出：