CODEWITHSUNDEEP
pythonscipyodedifferential-equations
ValientProcess
ValientProcess

Reputation: 1801

Solving ODE set with a step-function parameter using odeint

I'm trying to solve the following system numerically with odeint in Python:enter image description here:

My problem is that k1 depends on the value of p. When p<=0.01 k1=0 and otherwise k1=5.

I'm not sure how to write the system with the above constraint. I can start with:

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

Where flip can be the 0.01 value, but I'm not sure how to implement the if statement for this case, because I need the integration process to check the value of p after each iteration.

Upvotes: 1

Views: 1968

Answers (2)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

To reconstruct the values of k1 and p after the integration, put this computation into an extra function

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

and call this then in the ODE derivatives function

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

and later after the integration to plot along the solution components

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

Using the per comment corrected threshold value flip=-0.01 one obtains indeed oscillating behavior.

plot of solution and decision variables

Upvotes: 3

xdze2
xdze2

Reputation: 4151

Here is a solution using 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');

result graph

It seems to be a difficult problem to solve numerically. Other tested solvers don't converge, and the use of events is not helpful as they are more and more frequent near the steady state solution

Edit: indeed changing the flipping value to -0.01 instead of +0.01 leads to a limit cycle solution, for which the events method of solve_ivp is working:

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)

now the solution is:

solution with p=-0.01

Upvotes: 1

Related Questions