The Lorenz system of equations came about as a result of the study of atmospheric physics.
It is eponymously named after physicist Edward Lorenz, who in his quest to model trends in climatology and meteorology, stumbled on what is the field of chaotic dynamical systems.
In essence, the study of chaotic systems seeks to discover the underlying trends and correlations in otherwise inexplicable phenomena.
The field is vast and worth dedicating an entire lifetime of research. And there are those in academia who do just that.
However, what is arguably more valuable, is using these lessons gleamed in chaos theory, to uncover insights in other fields.
We shall be doing that today via the differential equation toolkit in Julia.
The following code samples are modified from the documentation below.
I particularly like Julia for tasks like this, for the insane amount of ease and convenience it provides when it comes to scientific computation.
It is only a matter of time before it becomes the king of scientific computing languages, overtaking both R and Python.
I shall be examining and explaining how the Julia code, line by line, how it works, and how we work through the following ordinary differential equations (ODEs).
I have a particular interest in differential equations for they arise practically in all fields of physics, and considerable amount of research in the realm of computer science is dedicated to developing algorithms to solve them.
The Lorenz system
At first glance, they appear simple enough.
Just 3 linear equations right? Right?
As we shall see later, they describe a rich system of correlations that are invisible to the human eye.
Julia Code
#In the Julia command line#
using Pkg
Pkg.add("DifferentialEquations")
Pkg.add("Plots")using DifferentialEquations
using Plots
Building the function.
Note the scientific notation. We can just type in the symbols directly as variables.
#Writing a parameterized functionfunction parameterized_lorenz!(du,u,p,t)
x, y, z = u #variables are part of vector array u σ, ρ, β, = p #coefficients are part of vector array p
du[1] = dx = σ*(y-x)
du[2] = dy = x*(ρ-z) - y
du[3] = dz = x*y - β*z
end.
.
.< parameterized_lorenz! (generic function with 1 method) >
Feeding the function the appropriate inputs.
#Initial conditions. x=1.0, y=0.0, z=0.0
u0 = [1.0, 0.0, 0.0]#Timespan of the simulation. 100 in this case.
tspan = (0.0, 100.0)#Coefficients of the functions.
p = [10.0, 28.0, 8/3]#Feeding the inputs to the solver prob = ODEProblem(parameterized_lorenz!, u0, tspan, p)
sol = solve(prob)
The numerical output of the solver:
The 3 elements in the numerical solution are as follows:
- “automatic order switching interpolation”
Different differential equations require different trade offs in terms of granularity and speed. The algorithm utilized can be further configured for different combinations. In this case we are using the default algorithm.
2. “t: 1242 element vector”.
The solution vector array returns an array with a granularity of 1242 time steps.
3. “u: 1242-element Vector”
Similarly the solution vector array returns a solutions for each variable at each moment in time.
retcode: Success
Interpolation: automatic order switching interpolation
t: 1242-element Vector{Float64}:
0.0
0.00020190266322792406
0.0022209292955071643
0.008092024514959373
0.017621941854931468
0.029910970730013992
0.046119886045152166
0.06633761003486688
0.09145306990996667
0.12186146157177435
0.15775265249378573
0.19917912194536747
0.23735551696369733
⋮
99.14124015678328
99.23005160643454
99.3178426998792
99.39089107520101
99.46975368622516
99.54941086677553
99.63372391029684
99.70543691911548
99.79383354268911
99.90535455139218
99.99916037878022
100.0
u: 1242-element Vector{Vector{Float64}}:
[1.0, 7.0, 0.0]
[1.0157322915791744, 7.0034699810052485, 0.001424457841035686]
[1.1713386826104133, 7.042273078544965, 0.016891096442499136]
[1.60812010145366, 7.196207031607893, 0.07444157353741501]
[2.2799815664984635, 7.5686982713074045, 0.2083531444865257]
[3.110193937388093, 8.256085485026166, 0.4611897477586723]
[4.20373407541577, 9.492552123156415, 0.9606655860252622]
[5.661209508205023, 11.525389422787784, 1.9400058922165015]
[7.749623385496089, 14.731104165543684, 3.9773408171241535]
[10.806830843846502, 19.246723699443496, 8.351046992399898]
[14.925214594925762, 23.671518939598823, 17.361338489508825]
[18.585923260717127, 22.195831734843853, 31.716044613009842]
[18.071904208857134, 12.17331991676706, 41.081896247131]
⋮
[-9.183457977837184, -12.682011018250817, 17.64472825931256]
[-12.321421331276596, -13.355987197450428, 26.506979959195224]
[-10.49509552307489, -6.614740940422851, 31.26114857062963]
[-6.412316590547746, -2.347421430959372, 28.77038858949141]
[-3.2434635259271523, -1.1729141982820335, 24.469982740650746]
[-1.942226612879031, -1.347551644896042, 20.53334362551526]
[-1.7799839534618225, -1.9964184236190687, 17.090940494285892]
[-2.2255447014614846, -2.969365928483296, 14.77883507206152]
[-3.5474340484603615, -5.177797852868663, 12.944117866721324]
[-7.152802707828483, -10.573895913693827, 14.093251538190133]
[-11.671892462551858, -14.878852040827208, 22.03887865658922]
[-11.706741167586157, -14.885100008376064, 22.141663040463467]If we were to plot a graph of our solution, this is what we will see
plot1 = plot(sol,
vars = (1,2,3),
xlabel="x",
ylabel="y",
zlabel="z",
label="x against y against z")plot2 = plot(sol, label = ["x" "y" "z"])plot(plot1, plot2, layout = (1, 2))
The aesthetically interesting blue plot is the famed Lorenz attractor, otherwise known as the butterfly plot.
Here we see the first glimpses of the nature of chaotic systems.
If we look at how they vary against time, it’s very difficult to spot the nature of the correlations with each other, and all we can see is that they vary erratically with time.
However, if we were to convert the plot into a phase space plot, we start to get a clearer picture of the nature of the correlations, if they even exist at all.
Extreme sensitivity to initial conditions
In the next segment, we shall demonstrate the key characteristic of chaotic systems. How the change of one variable, leads to a gigantic changes down the line.
We shall vary just one parameter, rho.
And see how the behaviour changes drastically.
u0 = [1.0, 7.0, 0.0]
tspan = (0.0, 100.0)
p = [13.0, 24.0, 7/3]
prob = ODEProblem(parameterized_lorenz!, u0, tspan, p)
sol = solve(prob);plot5 = plot(sol, vars = (1,2,3), xlabel="x", ylabel="y", zlabel="z",label="x against y against z",title="ρ = 24")plot6 = plot(sol, label = ["x" "y" "z"], title="ρ = 24")plot(plot5, plot6, layout = (1, 2))
Okay, that wasn’t so bad right? The plots don’t look they differ much from the original.
Now lets change, rho from 24 to 23
u0 = [1.0, 7.0, 0.0]
tspan = (0.0, 100.0)
p = [13.0, 23.0, 7/3]
prob = ODEProblem(parameterized_lorenz!, u0, tspan, p)
sol = solve(prob);plot3 = plot(sol, vars = (1,2,3), xlabel="x", ylabel="y", zlabel="z",label="x against y against z", title="ρ = 23")plot4 = plot(sol, label = ["x" "y" "z"], title="ρ = 23")plot(plot3, plot4, layout = (1, 2))
And boom, the behaviour changes drastically.
The system changes from fluctuating widely all over the place, into a behaviour that collapses into a single state.
Lets just view them side by side for ease of comparison.
plot(plot3, plot4, plot5, plot6, layout=(2,2))
Just one small change, and the system exhibits different characteristics. From oscillating wildly, out of control, to a stable and ordered system that goes into a stable state.
Conclusion
And that is the power of chaotic dynamics.
These variables may not correspond to actual, physical conditions that exist in reality. They could very well correspond to a complex interplay of factors that occur only once in a blue moon.
And finding the key factor, is often, the holy grail in modelling dynamical systems.
