Scipy.integrate solve_ivp 的稳定性问题

问题描述

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

当我使用 deSolve 包在 R 中执行此操作时，它工作正常，但我正在尝试使用 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[:,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])

输出：

解决方法

如果您在 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，导致一个非零的固定点{ {1}} 分量，M 为最佳估计值。

前两个分量的流图证实了这一假设，0.53371702 处的鞍点将收敛到稳定点 [0.31024394,0.22563217] 和 [0.07006604,0.53371702] 的区域分开

再次检查 R 和 python 版本中的代码是否确实相同。找出R中的积分方法，大概就是lsoda了。然后找出默认容差，或设置合理的容差（[0.59734069,0.0]），然后在 python 中复制。使用 atol = 1e-11,rtol=1e-13(method) 选项设置 ="BDF" 的积分方法。

numpy ode python scipy scipy