Monday, August 4, 2014

Happy Birthday Fortran (II)

In the previous article, I show a small piece of code. Let's have a look at a larger one.

In the sixties, Edward Lorenz, who was working at MIT, found that a weather model he built exhibited a very strange behavior: two sets of inputs, only different by minute quantities, would initially behave similarly and then diverge widely.

After he studied the phenomenon, he showed an example of another system that had the same property, a system which is known as Lorenz's attractor. Several texts have been published on the subject, let's just recap the main points:
  • The system oscillates between two points (attractors)
  • Two different starting points will lead to two different trajectories
  • A trajectory never repeats itself
As a computer represents real numbers with a finite precision, the second point is open to discussion: if two trajectories are so close that the difference between the closest points of each is below the machine precision, the rounding may result in the one of the trajectories being altered. This is actually what happens and started Lorenz investigations: while the internal representation of the numbers was up to 6 decimals, Lorenz's printout was only 4. When he stopped his experiment and restarted it later from an earlier point, the results started by being the same, then diverged. The field of chaos theory was created ...

My second homage to Fortran is a program that generates the coordinates for the points of the Lorenz's attractor. The code can be found in my GitHub repository.

Two methods are used to calculate the points: Explicit Euler (or Forward Euler), Implicit Euler (or Backward Euler) and a mixed method, which uses the average of a backward and forward Euler step.

The Explicit Euler is quite simple and straightforward to implement: the three differential equations are explicit and do not require solving anything special.

The Implicit Euler, on the other hand and by extension the Implicit-Explicit Euler, is a bit more challenging: each step requires solving a system of three equations with three unknowns. This is done - in my case - using Newton's method, which in turn requires computing the inverse of the Jacobian matrix, which is performed in my code using Gauss-Jordan decomposition: starting with the original matrix and the identity matrix, I apply all the operations needed to transform the original matrix into the identity matrix. As the same operations are applied to the second (starting as the identity) matrix, the result is that the original matrix is transform into the identity matrix, and the identity matrix is transformed into the inverse of the original matrix.

Both the Implicit and Explicit Euler are first order methods, the mixed Implicit-Explicit Euler is a second  order method: its convergence is faster and it is more stable. (Read here for a discussion of the Backward and Forward Euler methods)