Блог

Я использую SciPy для решения системы дифференциальных уравнений, и у меня возникли некоторые проблемы со стабильностью решения. По сути, существуют две переменные состояния, и они должны привести к двум относительно стабильным решениям, одно из которых асимптотически приближается к нулю (по сути, нулю), а другое-к 0,59.

Когда я делаю это в R с пакетом deSolve, это работает нормально, но я пытаюсь использовать scipy, и я получаю странную вещь, когда окончательное решение, переменные состояния нестабильны, даже если они должны быть? Значения переменной состояния, которые должны быть равны нулю, начинаются с ряда значений порядка 7,41 e-323, а затем на один шаг повышаются до 5,3 e-001, а затем возвращаются вниз. Я понятия не имею, почему это так, но мне интересно, может ли кто-нибудь дать какие-либо советы о том, как а) исправить это или б) использовать другой вариант, кроме 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])
 

Выход:

Если вы напечатаете первый десятичный знак после точки последнего M компонента в сетке C-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 .