Skip to content

Instantly share code, notes, and snippets.

@PM2Ring
Created January 4, 2023 11:59
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save PM2Ring/6f060da8f837cf124082cc5e0c288a89 to your computer and use it in GitHub Desktop.
Save PM2Ring/6f060da8f837cf124082cc5e0c288a89 to your computer and use it in GitHub Desktop.
Lorenz attractor
""" Lorenz attractor
From https://en.wikipedia.org/wiki/Lorenz_system
Modified for Sagemath by PM 2Ring 2023.01.04
"""
import numpy as np
from scipy.integrate import odeint
# Build Bezier curves
def bez(pos, tgt):
pos_tgt = zip(pos, tgt)
p0, t0 = next(pos_tgt)
row = [p0]
curves = []
for p, t in pos_tgt:
s = abs(p - p0) / 3
row.extend([p0 + t0*s, p - t*s, p])
curves.append(row)
row = []
p0, t0 = p, t
return curves
@interact
def _(rho=28, sigma=10, beta=8/3, hi=40, step=1/10, frame=True,
dark=True, color='skyblue', thickness=1, auto_update=False):
theme = 'dark' if dark else 'light'
def f(state, t):
x, y, z = state
# Derivatives
return vector([sigma * (y - x), x * (rho - z) - y, x * y - beta * z])
state0 = [1.0, 1.0, 1.0]
times = np.arange(0.0, hi, step)
states = odeint(f, state0, times)
states = list(map(vector, states))
grads = [f(state, t).normalized()
for state, t in zip(states, times)]
P = bezier3d(bez(states, grads), color=color, thickness=thickness)
P.show(frame=frame, theme=theme, projection='orthographic')
@PM2Ring
Copy link
Author

PM2Ring commented Jan 4, 2023

Here's a live demo running in SageMathCell.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment