Lorenz Attractor
Overview
The Lorenz Equations are a mathematical model that represents atmospheric convection through a system of ordinary differential equations.
System
$$ \begin{align*} {{dx} \over {dt}} =& - \sigma x + \sigma y \\ {{dy} \over {dt}} =& - xz + \rho x - y \\ {{dz} \over {dt}} =& xy - \beta z \end{align*} $$
Variables
- $x(t)$: Represents the $x$ coordinate of a particle at time $t$.
- $y(t)$: Represents the $y$ coordinate of a particle at time $t$.
- $z(t)$: Represents the $z$ coordinate of a particle at time $t$.
Parameters
- $\sigma$: The Prandtl Number, a parameter related to viscosity and thermal conductivity.
- $\rho$: The Rayleigh Number, a parameter related to the method of heat transfer in fluids.
- $\beta$: According to Wikipedia1, it seems to be a parameter related to dimensionality. Textbooks often do not explain what it exactly is.
Initially, these parameters are often set as $\sigma = 10$, $\rho = 28$, $\displaystyle \beta = {{8} \over {3}}$, following Lorenz’s convention. The particle movement described by the Lorenz attractor becomes chaotic when the Rayleigh Number exceeds $\rho \approx 24.74$. It responds sensitively to even minor changes in initial conditions, making prediction very difficult. It embodies the idea that the flap of a butterfly’s wings could potentially lead to a storm.
Bifurcation
The diagram above shows the bifurcation diagram of the Lorenz attractor. The horizontal axis gradually changes the value of $\rho$ from $25$ to $325$, while the vertical axis plots the peaks in the $z$ direction from the trajectory after a sufficiently long time. The movement appears very unstable when $r= \rho$ is small, but becomes stable and possesses a fixed trajectory as it increases beyond a certain level. The following explanation is recommended:
The graph above illustrates the trajectory when $\gamma=100$. Usually, if $\gamma$ is large enough, the movement may seem chaotic initially but becomes increasingly periodic over time.
The graph above shows only the particle trajectories after a long time, removing the initial ones. Although not completely periodic, the behavior becomes almost predictable.
When graphing $xz$, the trajectory looks as shown above. Initially chaotic, over time, the orbit stabilizes, and the peaks become distinctly around three. Even in the bifurcation diagram, the same interpretation can be offered.
When $\gamma = 28$, however, the motion remains chaotic over a sufficiently long time, resulting in a scattered, thickened appearance in the bifurcation diagram.
Conversely, when $\gamma = 300$, it shows two distinct peaks as below. In the bifurcation diagram, this results in two lines appearing for a single $\gamma$ value.
Code
MATLAB
Below is the MATLAB code to draw the trajectory of the Lorenz attractor.
sigma = 10;
beta = 8/3;
rho = 28;
f = @(t,v) [-sigma*v(1) + sigma*v(2); rho*v(1) - v(2) - v(1)*v(3); -beta*v(3) + v(1)*v(2)];
[t,v] = ode45(f,[0 10000],[1 1 1]); % Runge-Kutta 4th/5th order ODE solver
x=v(:,1); x=x(50000:end);
y=v(:,2); y=y(50000:end);
z=v(:,3); z=z(50000:end);
plot3(x,y,z)
%plot(x,z)
Fortran
Below is Fortran code that creates the source for the bifurcation diagram and returns it as an OUTPUT.csv
file. The method used to solve the ODE is RK4.
program LorenzBifurcation
implicit none
double precision :: sigma=10.d0, rho=28.d0, beta=8.d0/3.d0
Integer, parameter :: D = 3
double precision :: v(D)=1.d0, h=0.001d0, z1, z2, z3
integer :: i, j, k, preitr=10**4-1000, endstep=10**4
character :: answer
abstract interface
function func (v)
Integer, parameter :: D = 3
double precision :: v(D), func(D)
end function func
end interface
open(UNIT=34,FILE="OUTPUT.csv",ACTION="WRITE")
rho=5.d0
h=0.01d0
do while(rho<325)
print '(A,f6.4,A)', "Now", (rho-5.d0)/320.d0, "%"
do i=1, preitr
v = RK4(v,f)
end do
z1 = preitr
z2 = preitr
z3 = preitr
do i=preitr+1, endstep
v = RK4(v,f)
z1 = z2
z2 = z3
z3 = v(3)
if(z1<z2 .and. z2>z3) then
write (34,*) rho, ",", z2
end if
end do
v=1.d0
rho=rho+0.1d0
end do
contains
function RK4(v,f)
procedure (func) :: f
double precision :: v(D), V1(D), V2(D), V3(D), V4(D), RK4(D)
V1 = f(v)
V2 = f(v + (h/2)*V1)
V3 = f(v + (h/2)*V2)
V4 = f(v + h*V3)
RK4 = v + (h/6)*(V1 + 2*V2 + 2*V3 + V4)
end function RK4
function f(v)
double precision :: v(D), f(D)
f(1) = sigma*(v(2)-v(1))
f(2) = v(1)*(rho-v(3))-v(2)
f(3) = v(1)*v(2) - beta*v(3)
end function f
end program LorenzBifurcation
Below is MATLAB code that uses the OUTPUT.csv
file obtained above to draw the bifurcation diagram.
[file,path] = uigetfile('*.csv');
data = csvread([path,file]);
sz = 0.6;
scatter(data(:,1),data(:,2),sz,'.');
Julia
Below is Julia code to find the Lorenz attractor and draw it.
using DifferentialEquations
using Plots
function parameterized_lorenz!(du,u,p,t)
x,y,z = u
σ,ρ,β = p
du[1] = dx = σ*(y-x)
du[2] = dy = x*(ρ-z) - y
du[3] = dz = x*y - β*z
end
u0 = [1.0,0.0,0.0]
tspan = (0.0,100.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(parameterized_lorenz!,u0,tspan,p)
solution = solve(prob)
plot(solution, vars = (1,2,3))