Python: Lorenz

If you follow my blog or have devoted at least a minute to run through the headings, you should know that I am a fanatic of Chaos Theory. Or rather, I am crazy about messing with very simple systems, which are related to the ideas of sensitive dependence on initial conditions, unexpectedly complex behavior and other stuff that the words ‘mathematical chaos’ bring to one’s mind. The key word here is ‘simple’, since, because of very modest math skills, during my exploration I kept avoiding things involving any non-trivial mathematics. The most notable here is the fact that I didn’t even try to approach the Lorenz system, whose phase portrait, along with the plots of Mandelbrot set, may be considered the symbol of the Chaos Theory. However, beside playing with the Barnsley fern the last month I also made up enough guts to try and code the most well-known chaotic system.

What made me avoid the Lorenz equations becomes evident upon inspecting the narrow range of the chaotic systems shown on my webpage. All of them are unified against the Lorenz system by their discrete nature. On the other side, the latter is a continuous system, represented by three differential equations – something that will easily frighten one as dumb as I am. Modelling differential equations in a program might feel – no surprise! – different and at least seems more difficult than coding iterated systems of plain algebraic (linear or quadratic) equations. I was particularly afraid that the discretization methods, which are available to me, won’t be good enough to produce the desired behavior and to allow for observing the double loop of the attractor. This is another illustration of two obvious ideas. First, I have consistently failed learning any math at the university. The second notion is more widely applicable and much more important: one should always try to approach a problem and make at least the first steps in some direction instead of merely thinking over all the possible obstacles standing in their way. With the Lorenz equations it took me about a year to stop complaining and start acting.

Because I am absolutely not good at applied math, the set of methods, which I could use to turn the Lorenz equations into discrete – and thus easily programmable – form, was quite limited. Still, it turned out that the simplest one that I know suits my needs fairly well. This means, I have chosen the explicit finite difference method. Despite the fact that it is very simple, I will take a few lines to explain the way I use it.

x' = σ (y - x)

y' = x (ρ - z) - y

z' = x y - β z.

Having the above equations constituting the Lorenz system we can easily transform them into finite difference equations. All we have to do is substitute each derivation term with the corresponding finite difference term, that is Δx/Δt goes instead of x' – and the same way for the remaining two variables. This basically means that instead of using time derivatives we assume that our Δt is small enough to allow for accurate enough approximation of the derivative. The Δx term, in turn, might be expanded as (xn - xn-1) – the change in the x coordinate over the corresponding Δt time frame.

(xn - xn-1) / Δt = σ (y - x)

(yn - yn-1) / Δt = x (ρ - z) - y

(zn - zn-1) / Δt = x y - β z.

Due to the simplicity of the initial equations this step instantly gives us a discrete system. When iterated it produces a trajectory, which fairly approximates that of the continuous Lorenz system. To illustrate this let’s turn to the lovely Python. Our simple program need not deal with differential equations – we can start directly with the discretized system. However, having minor functional programming experience I could not miss a chance to introduce additional complexity. Because experimenting with the parameters of the system may be interesting we won’t simply hardcode the equations in the form of Python function. Instead, we write a function that takes the parameters and return us the actual function representing the system.

def lorenz(dt,sigma=10.0,beta=2.66667,ro=28.):
 def l(x,y,z):
  xn = y*dt*sigma + x*(1 - dt*sigma)
  yn = x*dt*(ro-z) + y*(1-dt)
  zn = x*y*dt + z*(1 - dt*beta)
  return (xn,yn,zn)
 return l

As you can see, in addition to the three σ, β and ρ arguments, which are the coefficients of the Lorenz equations, our function also takes a Δt argument that determines the size of the time interval used for discretization. The lesser it is, the more accurate approximation we get – at the cost of having to compute for more discrete time steps. As for the generated function, it fully mimics the above finite difference equations. Taking three arguments – the coordinates of the system at some moment – it outputs a 3-tuple, which represents the coordinates after Δt time units. With this function we can quickly produce a trajectory of the Lorenz system. All we need to have is the tuple of initial coordinates, the number of steps to go and… another function.

def trajectory(system,origin,steps):
 t = [origin]
 for _ in range(steps):
  t.append(system(*t[-1]))
 return t

The trajectory function is so simple that it’s hard to come up with any words to say about it. We merely get a vector of 3-tuples out of it, which is the trajectory that we wanted so much. What can one do with a sequence of 3-dimensional points? Right – one should plot it. Fortunately, Python has the matplotlib entirely available to us:

import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def plot(trajectory):
 fig = plt.figure()
 ax = fig.gca(projection='3d')
 ax.plot([x for (x,_,_) in trajectory],
  [y for (_,y,_) in trajectory],
  [z for (_,_,z) in trajectory])
 plt.show()


steps = 100000
dt = 1e-3

l = lorenz(dt)
t = trajectory(l,(1.0,0.0,0.0),steps)
plot(t)

That’s it! We have the Lorenz attractor right in our computer. What more can one desire? Definitely not much, but there is at least one more thing that I would like to do with these simple tools. Let us make sure that the system we have just created is actually one demonstrating sensitive dependence on initial conditions. We won’t show that with mathematical rigor – just do enough to see that the change in the starting point actually leads to an observably different trajectory. That is, we simply plot two trajectories and look whether the plots are easily distinguishable.

def plotMany(trajectories):
 fig = plt.figure()
 ax = fig.gca(projection='3d')
 for trajectory in trajectories:
  ax.plot([x for (x,_,_) in trajectory],
   [y for (_,y,_) in trajectory],
   [z for (_,_,z) in trajectory],
   label='trajectory from ' + 
   '({0:.3},{1:.3},{2:.3})'.format(*trajectory[0]))
 ax.legend()

 plt.show()


dx = 1e-4
tprime = trajectory(l,(1.0+dx,0.0,0.0),steps)
plotMany([t,tprime])

Precisely as one would expect we see two trajectories, which are clearly different. One could go further and calculate how much they differ. Even more, as papers (e.g. this one) by Nigel Crook et al suggest, it is possible to create a pattern recognition algorithm based on the Lorenz attractor – that’s one of the things I would like to try to accomplish. Still, these topics form the scope for potential future posts, while what we have now is enough for today.

Feel free to download the code and play with it. Keep in mind that in order to run it and get the plots, you need to have matplotlib.

CodeProject