%pylab inline
import odespy
Populating the interactive namespace from numpy and matplotlib
print odespy.list_available_solvers()
['AdamsBashMoulton2', 'AdamsBashMoulton3', 'AdamsBashforth2', 'AdamsBashforth3', 'AdamsBashforth4', 'Backward2Step', 'BackwardEuler', 'BogackiShampine', 'CashKarp', 'CrankNicolson', 'Dop853', 'Dopri5', 'DormandPrince', 'Euler', 'EulerCromer', 'Fehlberg', 'ForwardEuler', 'Heun', 'Leapfrog', 'LeapfrogFiltered', 'Lsoda', 'Lsodar', 'Lsode', 'Lsodes', 'Lsodi', 'Lsodis', 'Lsoibt', 'MidpointImplicit', 'MidpointIter', 'RK2', 'RK3', 'RK4', 'RKC', 'RKF45', 'RKFehlberg', 'Radau5', 'Radau5Explicit', 'Radau5Implicit', 'RungeKutta1', 'RungeKutta2', 'RungeKutta3', 'RungeKutta4', 'ThetaRule', 'Trapezoidal', 'Vode', 'lsoda_scipy', 'odefun_sympy']
def f(y,t,r,K):
return r*y*(1-y/K)
r,K=0.1,1
solver = odespy.Vode(f, f_args=(r,K))
solver.set_initial_condition(0.1)
t = linspace(0,100)
y,t = solver.solve(t)
plot(t,y);
t
array([ 0. , 2.04081633, 4.08163265, 6.12244898, 8.16326531, 10.20408163, 12.24489796, 14.28571429, 16.32653061, 18.36734694, 20.40816327, 22.44897959, 24.48979592, 26.53061224, 28.57142857, 30.6122449 , 32.65306122, 34.69387755, 36.73469388, 38.7755102 , 40.81632653, 42.85714286, 44.89795918, 46.93877551, 48.97959184, 51.02040816, 53.06122449, 55.10204082, 57.14285714, 59.18367347, 61.2244898 , 63.26530612, 65.30612245, 67.34693878, 69.3877551 , 71.42857143, 73.46938776, 75.51020408, 77.55102041, 79.59183673, 81.63265306, 83.67346939, 85.71428571, 87.75510204, 89.79591837, 91.83673469, 93.87755102, 95.91836735, 97.95918367, 100. ])
def f(y, t, r, K):
return r[0]*y[0]*(1-y.sum()/K[0]), r[1]*y[1]*(1-y.sum()/K[1])
r,K = (0.1,0.11),(1,1.1)
solver = odespy.Vode(f, f_args=(r,K))
solver.set_initial_condition((0.1,0.1))
t = linspace(0,1000,10000)
y,t = solver.solve(t)#, lambda u, t, step_no: any(u[step_no]<1e-6))
plot(t,y)
plot(t,y.sum(axis=1), label='sum')
legend(labels=['y0','y1'], loc='center right')
axhline(y=1, color='k', ls='--')
axhline(y=0.95, color='k', ls='--')
ylim(0,1.1);
plot(t,y)
xlim(0,50)
(0, 50)
print odespy.list_available_solvers()
['AdamsBashMoulton2', 'AdamsBashMoulton3', 'AdamsBashforth2', 'AdamsBashforth3', 'AdamsBashforth4', 'Backward2Step', 'BackwardEuler', 'BogackiShampine', 'CashKarp', 'CrankNicolson', 'Dop853', 'Dopri5', 'DormandPrince', 'Euler', 'EulerCromer', 'Fehlberg', 'ForwardEuler', 'Heun', 'Leapfrog', 'LeapfrogFiltered', 'Lsoda', 'Lsodar', 'Lsode', 'Lsodes', 'Lsodi', 'Lsodis', 'Lsoibt', 'MidpointImplicit', 'MidpointIter', 'RK2', 'RK3', 'RK4', 'RKC', 'RKF45', 'RKFehlberg', 'Radau5', 'Radau5Explicit', 'Radau5Implicit', 'RungeKutta1', 'RungeKutta2', 'RungeKutta3', 'RungeKutta4', 'ThetaRule', 'Trapezoidal', 'Vode', 'lsoda_scipy', 'odefun_sympy']
def f(y, OD, r, K, P, m, Dmax, Km):
N,F = y
Dn = Dmax * N/(N + F + Km * OD)
Df = Dmax * F/(N + F + Km * OD)
dOD = r * (1 - OD/K)
dN = P - m*N/OD - Dn
dF = m*N/OD - Df
return dN/dOD, dF/dOD
r, K, P, m, Dmax, Km = 0.1, 1, 250, 100, 230, 35
solver = odespy.Dopri5(f, f_args=(r, K, P, m, Dmax, Km))
solver.set_initial_condition((0.,0.))
OD = linspace(0.1,K*0.999,100)
y,OD = solver.solve(OD)#, lambda u, OD, step_no: any(u[step_no]<1e-6))
plot(OD,y)
legend(labels=['N','F'], loc='upper left')
xlabel('OD')
ylabel('F');
def FvsOD(OD, r, K, P, m, Dmax, Km):
solver = odespy.Dopri5(f, f_args=(r, K, P, m, Dmax, Km))
solver.set_initial_condition((0.1,0.))
y,OD = solver.solve(OD)
N,F = y
return F