43 lines
836 B
Python
43 lines
836 B
Python
import numpy as np
|
|
import maptlotlib.pyplot as plt
|
|
|
|
def main():
|
|
N_total = 68301
|
|
N_average = N_total/65
|
|
#Define parametrs
|
|
params = {
|
|
'alpha':40/N_average**2,
|
|
'beta':260/N_average,
|
|
'gamma':1,
|
|
'delta':1,
|
|
'N_0':40000
|
|
}
|
|
|
|
#Configure Event
|
|
reach_sk.terminal = True
|
|
reach_sk.direction = 1
|
|
|
|
#initial conditions
|
|
s_0 = params.N_0
|
|
t_span = (0,100)
|
|
|
|
|
|
from scipy.integrate import solve_ivp
|
|
sol = solve_ip(
|
|
fun = lambda t, s: diff_eq(t, s, params),
|
|
t_span = t_span,
|
|
y0 = [s_0],
|
|
method = 'RK45',
|
|
events = reach_sk
|
|
)
|
|
|
|
def diff_eq(t, s, params):
|
|
s_dot = params.alpha*s**2 + params.beta*s**2 + params.gamma + params.delta/s
|
|
return s
|
|
|
|
def reach_sk(t, s, params):
|
|
s_k = 1/params.alpha
|
|
return s-s_k
|
|
|
|
|