BZ Reaction

4271 days ago by medlock

def BZ(t, X, params): (x, y, z) = X (epsilon, delta, q, f) = params dx = (q * y - x * y + x * (1 - x)) / epsilon dy = (- q * y - x * y + 2 * f * z) / delta dz = x - z return [dx, dy, dz] def BZJacobian(t, X, params): (x, y, z) = X (epsilon, delta, q, f) = params dFxdx = (1 - 2 * x - y) / epsilon dFxdy = (q - x) / epsilon dFxdz = 0 dFydx = - y / delta dFydy = (- q - x) / delta dFydz = 2 * f / delta dFzdx = 1 dFzdy = 0 dFzdz = -1 dFxdt = 0 dFydt = 0 dFzdt = 0 return [[dFxdx, dFxdy, dFxdz], [dFydx, dFydy, dFydz], [dFzdx, dFzdy, dFzdz], [dFxdt, dFydt, dFzdt]] 
       
epsilon = 5e-5 delta = 2e-4 q = 8e-4 f = 0.5 
       
T = ode_solver() T.algorithm = 'bsimp' T.function = BZ T.jacobian = BZJacobian T.y_0 = [1, 1, 1] T.ode_solve(t_span = [0, 30], num_points = 2001, params = [epsilon, delta, q, f]) 
       
t = [S[0] for S in T.solution] x = [S[1][0] for S in T.solution] y = [S[1][1] for S in T.solution] z = [S[1][2] for S in T.solution] logx = map(log, x) logy = map(log, y) logz = map(log, z) 
       
line(zip(t, logx), color = 'blue') + line(zip(t, logy), color = 'red') + line(zip(t, logz), color = 'green') 
       
line3d(zip(logx, logy, logz), color = 'red') + line3d(zip(logx[-500 : ], logy[-500 : ], logz[-500 : ]), color = 'blue', thickness = 2) + point3d([logx[0], logy[0], logz[0]], color = 'red') 
       
line(zip(logx[-500 : ], logy[-500 : ])) 
       
line(zip(logx[-500 : ], logz[-500 : ])) 
       
line(zip(logy[-500 : ], logz[-500 : ]))