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 : ]))
|
|
|
|