Решение ODE, заданного параметром пошаговой функции, с использованием odeint

Question

Я пытаюсь численно решить следующую систему с помощью odeint в Python: :

Моя проблема в том, что k1 зависит от значения p. Когда p<=0.01 k1=0 и в противном случае k1=5.

Я не уверен, как написать систему с вышеуказанным ограничением. Я могу начать с:

def sys(x,t,flip):
S,T,p = x

dx = [1-T-k1*T,
     0.1*(1-S)-k1*S,
     -0.2*(1-T-k1*T)+0.1*(1-S)-k1*S]
return dx

Где flip может быть значением 0.01, но я не уверен, как реализовать оператор if для этого случая, потому что мне нужен процесс интеграции для проверкизначение p после каждой итерации.

Answer 1

Чтобы восстановить значения k1 и p после интегрирования, поместите это вычисление в дополнительную функцию

def kp(T,S):
    p = -0.2*T + S
    k1lo, k1hi = 0, 5
    flip = -0.01
    k1 = 0.5*(k1lo+k1hi + (k1hi-k1lo)*np.tanh((p-flip)/1e-8))
    return k1,p

и вызовите ее затем в производной функции ODE

def ode(t, TS):
    T, S = TS
    k1,_ = kp(T,S)
    dTdt = 1 - T - k1*T
    dSdt = 0.1*(1 - S) - k1*S
    return dTdt, dSdt

и позже после интегрирования для построения графика вдоль компонентов решения

T,S = sol.y;
k,p = k1(T,S);

Используя скорректированное пороговое значение для каждого комментария flip=-0.01, можно получить действительно колебательное поведение.

Answer 2

Вот решение, использующее solve_ivp :

import numpy as np
import matplotlib.pylab as plt

from scipy.integrate import solve_ivp

def ode(t, TS):
    T, S = TS

    p = -0.2*T + S
    k1 = 0 if p <= 0.01 else 5

    dTdt = 1 - T - k1*T
    dSdt = 0.1*(1 - S) - k1*S

    return dTdt, dSdt

# Solve
TS_start = (0.7, 0.4)
t_span = [0, 3]
sol = solve_ivp(ode, t_span, TS_start, method='RK45',
                rtol=1e-5)
print(sol.message)
print('nbr eval', sol.nfev)

# Graph
plt.plot(sol.t, sol.y[0, :], label='T')
plt.plot(sol.t, sol.y[1, :], label='S'); plt.legend();
plt.xlabel('time');

Кажется, трудно решить проблему численно,Другие проверенные решатели не сходятся, и использование событий бесполезно, так как они все чаще встречаются вблизи стационарного решения

Редактировать: действительно, изменение значения переворачивания на -0,01 вместо +0,01 приводит крешение с предельным циклом, для которого работает events метод solve_ivp:

import numpy as np
import matplotlib.pylab as plt

from scipy.integrate import solve_ivp

def ode(t, TS):
    T, S = TS

    p = -0.2*T + S
    k1 = 0 if sign_change(t, TS) <= 0 else 5

    dTdt = 1 - T - k1*T
    dSdt = 0.1*(1 - S) - k1*S

    return dTdt, dSdt

def sign_change(t, TS):
    T, S = TS
    p = -0.2*T + S
    flip = -0.01
    return p - flip

# Solve
TS_start = [0.3, 0.1]
t_span = [0, 5]
sol = solve_ivp(ode, t_span, TS_start,
                method='RK45', events=sign_change)
print(sol.message)
print('nbr eval', sol.nfev)

# Graph
plt.plot(sol.t, sol.y[0, :], '-x', label='T')
plt.plot(sol.t, sol.y[1, :], label='S'); plt.legend();
plt.xlabel('time');

for t in sol.t_events[0]:
    plt.axvline(x=t, color='gray', linewidth=1)

Теперь решение: