a collection of lengthy notes on how I am usually coward to do what is wrong and how I am sometimes brave to do what is right. Bits of management and self management, .NET programming, Clojure parenthesis and lots of design mistakes are the things one will find here.

воскресенье, 25 августа 2013 г.

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.

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.

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:

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.

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.

Комментариев нет:

Отправить комментарий

About

I work as a team-lead at Acumatica, the Cloud ERP. Besides, I hack code on my own for the sake of learning, fun and whatever else that programming can bring. To reflect on these activities and structure my thoughts and ideas better I also maintain a little blog.