📂Dynamics

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 Wikipedia¹, 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))