The Lorenz Attractor: Exploring Chaos with Numerical Simulations

Publish on: Classify at:
Words: 3125 Read:≈ 15min
Summary
This blog explores the Lorenz attractor, a fascinating example of chaos theory, through numerical simulations. It delves into the mathematical foundations, explains key parameters, and demonstrates visualization techniques using various programming tools like Julia, Python, MATLAB, and COMSOL. The post compares methods like Euler and Dormand-Prince, highlighting their strengths and weaknesses, while showcasing the beauty and complexity of chaotic systems.

Introduction

Lorenz attractor is a 3D structure that visualizes chaotic behavior, originally discovered by Edward Lorenz in 1963 while he was studying atmospheric convection (weather patterns). It is famous for being a “strange attractor.” This means that while the system is deterministic (it follows strict rules), its long-term behavior is unpredictable and never repeats itself. The trajectory loops around two points, creating a shape that famously resembles a butterfly.

From a technical standpoint, the Lorenz system is nonlinear, aperiodic, three-dimensional, and deterministic. While originally for weather, the equations have since been found to model behavior in a wide variety of systems, including lasers, dynamos, electric circuits, and even some chemical reactions

Here is the famous system:

dxdt=σ(y−x),dydt=x(ρ−z)−y,dzdt=xy−βz. \begin{aligned} \frac{dx}{dt} &= \sigma (y - x), \\ \frac{dy}{dt} &= x (\rho - z) - y, \\ \frac{dz}{dt} &= xy - \beta z. \end{aligned}

These equations describe how the state of the atmosphere changes over time (tt).

  • xx: The rate of convection (how fast the fluid is moving, proportional to the intensity of the convection).
  • yy: The horizontal temperature difference (proportional to the temperature difference between the rising and falling air currents).
  • zz: The vertical temperature gradient (proportional to the distortion of the vertical temperature profile).

The constants σ,ρ,β\sigma, \rho, \beta are parameters representing physical properties of the system.

  • σ\sigma is the Prandtl number . This is essentially the ratio of the fluid’s viscosity to its thermal conductivity. A high number means the fluid is thick and holds onto momentum (like oil), while a low number means heat spreads through it faster than it can physically move (like mercury).
  • ρ\rho is the  Rayleigh number . It represents the temperature difference between the hot bottom and the cold top.
  • β\beta relates to the physical dimensions of the fluid layer itself. It describes the aspect ratio (width vs. height) of the circular currents in the fluid.

Visualization

When ρ=28\rho=28, σ=10\sigma=10, and β=83\beta=\dfrac{8}{3}, the Lorenz system has chaotic solutions (but not all solutions are chaotic). Almost all initial points will tend to an invariant set — the Lorenz attractor — which is a strange attractor, a fractal, and a self-excited attractor with respect to all three equilibria. To simulate the Lorenz attractor, we use numerical methods to solve the equations. Two common methods are:

  1. Euler Method: A simple and fast method, but less accurate for chaotic systems.
  2. Dormand-Prince Method: A more advanced method (used in ode45 and RK45) that adapts the step size for better accuracy.

Below, we demonstrate how to visualize the Lorenz attractor using these methods in different programming languages and simulation tools.

Code Examples

COMSOL-Dormand Prince 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

% Lorenz Attractor using COMSOL with Dormand-Prince Method

import com.comsol.model.*
import com.comsol.model.util.*

model = ModelUtil.create('Model');

% ==================== 1. Global Parameters ====================
params = model.param;
params.set('sigma', '10');
params.set('rho', '28');
params.set('beta', '8/3');
params.set('dt', '0.02');
params.set('Tfinal', '30');

% ==================== 2. Component and Physics ====================
model.component.create('comp1', true);
comp = model.component('comp1');
comp.geom.create('geom1', 3);

comp.physics.create('ge', 'GlobalEquations', 'geom1');
ge = comp.physics('ge').feature('ge1');
ge.set('name', {'u'; 'v'; 'w'});
ge.set('equation', {'ut-sigma*(v-u)'; 'vt-rho*u+v+u*w'; 'wt-u*v+beta*w'});
ge.set('initialValueU', [2; 10; 10]);

% ==================== 3. Study and Solver ====================
model.study.create('std1');
study = model.study('std1');
study.create('time', 'Transient');
study.feature('time').set('tlist', 'range(0,dt,Tfinal)');

model.sol.create('sol1');
sol = model.sol('sol1');
sol.attach('std1');
sol.create('st1', 'StudyStep');
sol.create('v1', 'Variables');
sol.create('t1', 'Time');

timeSolver = sol.feature('t1');
timeSolver.set('tlist', 'range(0,dt,Tfinal)');
timeSolver.set('odesolvertype', 'explicit');
timeSolver.set('rkmethod', 'dopri5');

study.runNoGen;

% ==================== 4. Post-processing ====================
model.result.create('pg1', 'PlotGroup3D');
pg = model.result('pg1');
pg.set('looplevel', {'last'});
pg.create('traj1', 'PointTrajectories');

traj = pg.feature('traj1');
traj.set('expr', {'u', 'v', 'w'});
traj.set('linetype', 'tube');
traj.set('tuberadiusscaleactive', true);
traj.set('tuberadiusscale', 0.1);
traj.create('col1', 'Color');
traj.feature('col1').set('expr', 'sqrt(ut^2+vt^2+wt^2)');

% ==================== 5. Save Model ====================
mphsave(model, fullfile(pwd, 'Lorenz_attractor_DormandPrince.mph'));

Julia - Euler

Julia Euler
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

using Plots

# Define the Lorenz attractor
@kwdef mutable struct Lorenz
    dt::Float64 = 0.02
    σ::Float64 = 10
    ρ::Float64 = 28
    β::Float64 = 8 / 3
    x::Float64 = 2
    y::Float64 = 1
    z::Float64 = 1
end

# Update the state of the Lorenz attractor
function step!(l::Lorenz)
    dx = l.σ * (l.y - l.x)
    dy = l.x * (l.ρ - l.z) - l.y
    dz = l.x * l.y - l.β * l.z
    l.x += l.dt * dx
    l.y += l.dt * dy
    l.z += l.dt * dz
end

# Create an instance of the Lorenz attractor
attractor = Lorenz()

# Initialize a 3D plot
plt = plot3d(
    1,
    xlim=(-30, 30),
    ylim=(-30, 30),
    zlim=(0, 60),
    title="Lorenz Attractor",
    marker=2,
    legend=false,
    xlabel="X Axis",
    ylabel="Y Axis",
    zlabel="Z Axis",
    size=(1200, 800),
)

# Build an animated GIF
anim = @animate for _ = 1:1500
    step!(attractor)
    push!(plt, attractor.x, attractor.y, attractor.z)
end every 10

gif(anim, joinpath(pwd(), "lorenz_attractor_euler_jl.gif"), fps=30)

Julia - Dormand Prince 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

using Plots
using DifferentialEquations

# Define the Lorenz attractor
@kwdef mutable struct Lorenz
    dt::Float64 = 0.02
    σ::Float64 = 10
    ρ::Float64 = 28
    β::Float64 = 8 / 3
    x::Float64 = 2
    y::Float64 = 1
    z::Float64 = 1
end

# Define the Lorenz system for DifferentialEquations.jl
function lorenz!(du, u, p, t)
    σ, ρ, β = p
    du[1] = σ * (u[2] - u[1])          # dx/dt
    du[2] = u[1] * (ρ - u[3]) - u[2]   # dy/dt
    du[3] = u[1] * u[2] - β * u[3]     # dz/dt
end

# Create an instance of the Lorenz attractor
attractor = Lorenz()

# Define initial conditions and solve the Lorenz system
u0 = [attractor.x, attractor.y, attractor.z]
p = (attractor.σ, attractor.ρ, attractor.β)
tspan = (0.0, 1500 * attractor.dt)
prob = ODEProblem(lorenz!, u0, tspan, p)
sol = solve(prob, DP5(), saveat=attractor.dt)

# Initialize a 3D plot
plt = plot3d(
    1,
    xlim=(-30, 30),
    ylim=(-30, 30),
    zlim=(0, 60),
    title="Lorenz Attractor",
    marker=2,
    legend=false,
    xlabel="X Axis",
    ylabel="Y Axis",
    zlabel="Z Axis",
    size=(1200, 800),
)


# Build an animated GIF
anim = @animate for _ = 1:1500
    step!(attractor)
    push!(plt, attractor.x, attractor.y, attractor.z)
end every 10

gif(anim, joinpath(pwd(), "lorenz_attractor_dormandprince_jl.gif"), fps=30)

Python - Euler

Python Euler
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation


# Define the Lorenz attractor with default parameters
class Lorenz:
    def __init__(self, dt=0.02, sigma=10, rho=28, beta=8 / 3, x=2, y=1, z=1):
        self.dt = dt
        self.sigma = sigma
        self.rho = rho
        self.beta = beta
        self.x = x
        self.y = y
        self.z = z

    # Update the state of the Lorenz attractor using the Euler method
    def step(self):
        dx = self.sigma * (self.y - self.x)
        dy = self.x * (self.rho - self.z) - self.y
        dz = self.x * self.y - self.beta * self.z
        self.x += self.dt * dx
        self.y += self.dt * dy
        self.z += self.dt * dz
        return self.x, self.y, self.z


# Create an instance of the Lorenz attractor
attractor = Lorenz()

# Initialize a 3D plot
fig = plt.figure(figsize=(12, 8), dpi=100)
ax = fig.add_subplot(111, projection="3d")
ax.set_xlim(-30, 30)
ax.set_ylim(-30, 30)
ax.set_zlim(0, 60)
ax.set_title("Lorenz Attractor")
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")

# Initialize the plot with an empty series
(line,) = ax.plot([], [], [], lw=0.8, color="#007acc", alpha=0.9)
(points,) = ax.plot(
    [],
    [],
    [],
    marker="o",
    markersize=3,
    markerfacecolor="#00d2ff",
    markeredgecolor="black",
    markeredgewidth=0.6,
    alpha=0.9,
)

# Data storage for the animation
x_data, y_data, z_data = [], [], []


# Update the plot for each frame
def update(frame):
    global x_data, y_data, z_data
    x, y, z = attractor.step()
    x_data.append(x)
    y_data.append(y)
    z_data.append(z)
    line.set_data(x_data, y_data)
    line.set_3d_properties(z_data)
    points.set_data(x_data, y_data)
    points.set_3d_properties(z_data)
    return line, points


# Create an animated GIF
ani = FuncAnimation(fig, update, frames=1500, interval=10, blit=True)

# Save the animation as a GIF
ani.save("lorenz_attractor_euler_py.gif", writer="pillow", fps=30)

Python - Dormand Prince 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from matplotlib.animation import FuncAnimation


# Define the Lorenz attractor with default parameters
class Lorenz:
    def __init__(self, dt=0.02, sigma=10, rho=28, beta=8 / 3, x=2, y=1, z=1):
        self.dt = dt
        self.sigma = sigma
        self.rho = rho
        self.beta = beta
        self.x = x
        self.y = y
        self.z = z

    # Function to define the Lorenz system for solve_ivp
    def equations(self, t, state):
        x, y, z = state
        dx = self.sigma * (y - x)
        dy = x * (self.rho - z) - y
        dz = x * y - self.beta * z
        return [dx, dy, dz]


# Create an instance of the Lorenz attractor
attractor = Lorenz()

# Define initial conditions and time span for the solver
u0 = [attractor.x, attractor.y, attractor.z]
t_span = (0, 1500 * attractor.dt)
t_eval = np.arange(0, 1500 * attractor.dt, attractor.dt)

# Solve the Lorenz system using Dormand-Prince method
sol = solve_ivp(attractor.equations, t_span, u0, t_eval=t_eval, method="RK45")

# Extract solution points
x_data, y_data, z_data = sol.y

# Initialize a 3D plot
fig = plt.figure(figsize=(12, 8), dpi=100)
ax = fig.add_subplot(111, projection="3d")
ax.set_xlim(-30, 30)
ax.set_ylim(-30, 30)
ax.set_zlim(0, 60)
ax.set_title("Lorenz Attractor")
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")

# Initialize the plot with an empty series
(line,) = ax.plot([], [], [], lw=0.8, color="#007acc", alpha=0.9)
(points,) = ax.plot(
    [],
    [],
    [],
    marker="o",
    markersize=3,
    markerfacecolor="#00d2ff",
    markeredgecolor="black",
    markeredgewidth=0.6,
    alpha=0.9,
)


# Function to update the plot for each frame
def update(frame):
    line.set_data(x_data[:frame], y_data[:frame])
    line.set_3d_properties(z_data[:frame])
    points.set_data(x_data[:frame], y_data[:frame])
    points.set_3d_properties(z_data[:frame])
    return line, points


# Create an animated GIF
ani = FuncAnimation(fig, update, frames=len(t_eval), interval=10, blit=True)
ani.save("lorenz_attractor_dormandprince_py.gif", writer="pillow", fps=30)

MATLAB - Euler

MATLAB Euler
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

% Lorenz Attractor using Euler Method

% Define parameters for the Lorenz attractor
dt = 0.02;
sigma = 10;
rho = 28;
beta = 8/3;

% Initial conditions
x = 2;
y = 1;
z = 1;

% Number of steps for the simulation
num_steps = 1500;

% Initialize arrays to store the trajectory
x_vals = zeros(1, num_steps);
y_vals = zeros(1, num_steps);
z_vals = zeros(1, num_steps);

% Store initial conditions
x_vals(1) = x;
y_vals(1) = y;
z_vals(1) = z;

% Euler method to update the Lorenz attractor
for i = 2:num_steps
    dx = sigma * (y - x);
    dy = x * (rho - z) - y;
    dz = x * y - beta * z;

    x = x + dt * dx;
    y = y + dt * dy;
    z = z + dt * dz;

    x_vals(i) = x;
    y_vals(i) = y;
    z_vals(i) = z;
end

% Plot the Lorenz attractor
figure;
plot3(x_vals, y_vals, z_vals, 'LineWidth', 1.5);
axis equal;
view(30, 30);
grid on;
xlim([-30, 30]);
ylim([-30, 30]);
zlim([0, 60]);
title('Lorenz Attractor');
xlabel('X Axis');
ylabel('Y Axis');
zlabel('Z Axis');

set(gcf, 'Position', [100, 100, 1200, 800]);
saveas(gcf, fullfile(pwd, 'lorenz_attractor_euler_matlab.png'));

MATLAB - Dormand Prince

MATLAB DP
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

%% Lorenz Attractor using Dormand-Prince Method (ode45)

% Define parameters for the Lorenz attractor
dt = 0.02;
sigma = 10;
rho = 28;
beta = 8/3;

% Initial conditions
x0 = 2;
y0 = 1;
z0 = 1;

% Define the Lorenz system
lorenz = @(t, u) [
                  sigma * (u(2) - u(1)); % dx/dt
                  u(1) * (rho - u(3)) - u(2); % dy/dt
                  u(1) * u(2) - beta * u(3) % dz/dt
                  ];

% Define initial conditions and solve the Lorenz system
u0 = [x0; y0; z0];
tspan = [0, 1500 * dt];
opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-9);
[t, sol] = ode45(lorenz, tspan, u0, opts);

% Extract the trajectory
x_vals = sol(:, 1);
y_vals = sol(:, 2);
z_vals = sol(:, 3);

% Plot the Lorenz attractor
figure;
plot3(x_vals, y_vals, z_vals, 'LineWidth', 1.5);
axis equal;
grid on;
xlim([-30, 30]);
ylim([-30, 30]);
zlim([0, 60]);
daspect([1 1 1]);
view(30, 30);
title('Lorenz Attractor');
xlabel('X Axis');
ylabel('Y Axis');
zlabel('Z Axis');
set(gcf, 'Position', [100, 100, 1200, 800]);
saveas(gcf, fullfile(pwd, 'lorenz_attractor_dormandprince_matlab.png'));

Maple

Maple
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

# Lorenz Attractor: using odeplot

restart;
with(plots):

# Parameters
dt := 0.02:
sigma := 10:
rho := 28:
beta := 8/3:

# Initial conditions
x0 := 2:
y0 := 1:
z0 := 1:

# Define and solve the system
lorenz_sys := {
    diff(x(t), t) = sigma * (y(t) - x(t)),
    diff(y(t), t) = x(t) * (rho - z(t)) - y(t),
    diff(z(t), t) = x(t) * y(t) - beta * z(t),
    x(0) = x0, y(0) = y0, z(0) = z0
}:

sol := dsolve(lorenz_sys, numeric, method = rkf45, maxfun = 0):

# Create static 3D plot
odeplot(sol, [x(t), y(t), z(t)], 0 .. 30,
    axes = boxed,
    labels = ["X Axis", "Y Axis", "Z Axis"],
    title = "Lorenz Attractor",
    color = blue,
    thickness = 1,
    numpoints = 1500,
    orientation = [45, 70],
    view = [-30 .. 30, -30 .. 30, 0 .. 60]
);

Mathematica

Mathematica
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

(* Lorenz Attractor using Dormand-Prince method in Mathematica *)

(* Define the Lorenz attractor parameters *)
dt = 0.02;
σ = 10;
ρ = 28;
β = 8/3;

(* Initial conditions *)
x0 = 2;
y0 = 1;
z0 = 1;

(* Define the Lorenz system *)
lorenzSystem = {
    x'[t] == σ (y[t] - x[t]),           (* dx/dt *)
    y'[t] == x[t] (ρ - z[t]) - y[t],    (* dy/dt *)
    z'[t] == x[t] y[t] - β z[t]         (* dz/dt *)
};

(* Initial conditions *)
initialConditions = {x[0] == x0, y[0] == y0, z[0] == z0};

(* Time span *)
tspan = 1500 dt;

(* Solve the Lorenz system using NDSolve with ExplicitRungeKutta (Dormand-Prince) *)
sol = NDSolve[
    Join[lorenzSystem, initialConditions],
    {x, y, z},
    {t, 0, tspan},
    Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5},
    MaxSteps -> Infinity
];

(* Create the 3D plot *)
ParametricPlot3D[
    Evaluate[{x[t], y[t], z[t]} /. sol],
    {t, 0, tspan},
    PlotRange -> {{-30, 30}, {-30, 30}, {0, 60}},
    PlotLabel -> "Lorenz Attractor",
    AxesLabel -> {"X Axis", "Y Axis", "Z Axis"},
    PlotStyle -> Blue,
    ImageSize -> {1200, 800},
    BoxRatios -> {1, 1, 1},
    PlotPoints -> 1500
]

R

R
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

# Load required packages
library(deSolve)
library(plotly)

# Define the Lorenz attractor parameters
dt <- 0.02
sigma <- 10
rho <- 28
beta <- 8 / 3

# Initial conditions
x0 <- 2
y0 <- 1
z0 <- 1

# Define the Lorenz system
lorenz <- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
        dx <- sigma * (y - x) # dx/dt
        dy <- x * (rho - z) - y # dy/dt
        dz <- x * y - beta * z # dz/dt
        list(c(dx, dy, dz))
    })
}

# Set up parameters and initial state
parameters <- c(sigma = sigma, rho = rho, beta = beta)
state <- c(x = x0, y = y0, z = z0)

# Time span
times <- seq(0, 1500 * dt, by = dt)

# Solve using Dormand-Prince method
sol <- ode(
    y = state,
    times = times,
    func = lorenz,
    parms = parameters,
    method = "ode45"
)

# Convert to data frame
sol_df <- as.data.frame(sol)

# Create interactive 3D plot using plotly
fig <- plot_ly(
    sol_df,
    x = ~x,
    y = ~y,
    z = ~z,
    type = "scatter3d",
    mode = "lines",
    line = list(color = "blue", width = 2)
) %>%
    layout(
        title = "Lorenz Attractor",
        scene = list(
            xaxis = list(title = "X Axis", range = c(-30, 30)),
            yaxis = list(title = "Y Axis", range = c(-30, 30)),
            zaxis = list(title = "Z Axis", range = c(0, 60))
        )
    )

# Display the plot
fig

Comparison of Tools and Methods

MethodAdvantagesDisadvantages
EulerSimple and fastLess accurate for chaotic systems
Dormand-PrinceAdaptive step size, more accurateSlower, more complex implementation
COMSOLBuilt-in solver, easy visualizationRequires proprietary software
JuliaHigh performance, concise syntaxRequires familiarity with Julia
PythonWidely used, rich librariesSlower than compiled languages
MATLABExcellent for numerical computationExpensive, less flexible for animation
RGood for statistical analysisLimited support for 3D visualization

The choice of method depends on the trade-off between speed, accuracy, and the tools available.

Conclusion

The Lorenz attractor is a beautiful example of chaos theory, demonstrating how simple equations can produce complex, unpredictable behavior. In this blog, we explored various ways to visualize the attractor using different numerical methods and programming tools.

Each method has its strengths and weaknesses, and the choice depends on your goals and resources. Whether you’re a beginner or an expert, experimenting with the Lorenz attractor is a great way to learn about dynamical systems and numerical simulation.

Reference