1

I have a piece of code which increments a time-step for the so-called Lorenz95 model (invented by Ed Lorenz in 1995). It's normally implemented as a 40-variable model, and displays chaotic behaviour. I have coded up the time-stepping for the algorithm as follows:

class Lorenz:
    '''Lorenz-95 equation'''

    global F, dt, SIZE
    F = 8
    dt = 0.01
    SIZE = 40

    def __init__(self):
        self.x = [random.random() for i in range(SIZE)]

    def euler(self):
        '''Euler time stepping'''
        newvals = [0]*SIZE
        for i in range(SIZE-1):
            newvals[i] = self.x[i] + dt * (self.x[i-1] * (self.x[i+1] - self.x[i-2]) - self.x[i] + F)
        newvals[SIZE-1] = self.x[SIZE-1] + dt * (self.x[SIZE-2] * (self.x[0] - self.x[SIZE-3]) - self.x[SIZE-1] + F)
        self.x = newvals

This function euler is not slow, but unfortunately, my code needs to make a very large number of calls to it. Is there a way I could code the time-stepping to make it run faster?

Many thanks.

Improve this question
3
  • 5
    What the heck is global doing in there? Commented May 8, 2013 at 15:44
  • 1
    You could unroll the for loop and spell out the series of computations it does -- but there's no getting around the fact that they each have to be done. Might want to consider writing a C extension to do the actual number crunching. Commented May 8, 2013 at 16:14
  • As we all know the use of global is a sign of lazy programming. ;-) I want to make sure that all instances of the class Lorenz have the same size for their data array. So, forcing SIZE to be a global in this class is the easiest way of ensuring this. Commented May 9, 2013 at 8:20

1 Answer 1

Reset to default
3

There are at least two kinds of possible optimizations: working in a smarter way (algorithmic improvements) and working faster.

  • In the algorithmic side, you're using the Euler method, which is a first-order method (so global error is proportional to the step size) and has a smallish stability region. That is, it's not very efficient.

  • In the other side, if you're using the standard CPython implementation this kind of code is going to be quite slow. To get around that, you could simply try running it under PyPy. Its Just-in-Time compiler can make numerical code run maybe 100x faster. You could also write a custom C or Cython extension.

But there's a better way. Solving systems of ordinary differential equations is quite common, so scipy, one of the core scientific libraries in Python wraps fast, battle-tested Fortran libraries to solve them. By using scipy, you get both the algorithmic improvemnts (as integrators will have a higher order) and a fast implementation.

Solving the Lorenz 95 model for a set of perturbated initial conditions looks like this:

import numpy as np


def lorenz95(x, t):
    return np.roll(x, 1) * (np.roll(x, -1) - np.roll(x, 2)) - x + F

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    from scipy.integrate import odeint
    SIZE = 40
    F = 8
    t = np.linspace(0, 10, 1001)
    x0 = np.random.random(SIZE)
    for perturbation in 0.1 * np.random.randn(5):
        x0i = x0.copy()
        x0i[0] += perturbation
        x = odeint(lorenz95, x0i, t)
        plt.plot(t, x[:, 0])
    plt.show()

And the output (setting np.random.seed(7), yours can be different) is nicely chaotic. Small perturbations in the initial conditions (in just one of he coordinates!) produce very different solutions: Lorenz-95 dynamical system

But, is it really faster than Euler time stepping? For dt = 0.01 it seems almost three times faster, but the solutions don't match except at the very beginning. Euler vs odeint

If dt is reduced, the solution provided by the Euler method gets increasingly similar to the odeint solution, but it takes much longer. Notice how the smaller dt, the later Euler solutions loose track of the odeint solution. The most precise Euler solution took 600x longer to computed the solution up to t=6 than odeint up to t=10. See the full script here. Euler vs odeint

In the end, this system is so unstable that I guess not even the odeint solution is accurate along all plotted time.

Improve this answer
Sign up to request clarification or add additional context in comments.

5 Comments

Thanks @jorgeca for the suggestions. Your code makes a couple of important changes: 1. The use of np.roll means that the whole array can be updated in a single sweep, so I now have a new Euler time-stepping scheme: def euler2(self): self.x = self.x + dt * (np.roll(self.x, 1) * (np.roll(self.x, -1) - np.roll(self.x, 2)) - self.x + F). This is around 3 times faster than the original implementation. 2. The use of odeint, rather than a crude Euler time-step. This maybe more stable, but is quite a bit slower (2 times slower than the Euler routine). I guess I shall have to use compiled code.
@NeillBowler The faster results from your Euler routine are spurious except at the very very beginning. I'll expand my answer maybe tomorrow to show you that, but the summary is "to achieve the same error in the solution, Euler time stepping will take way many orders of magnitude more time than odeint".
@NeillBowler I've expanded the answer.
Many thanks for your great response. I've looked into the accuracy of the solver. Most papers using this model (including the [journals.ametsoc.org/doi/abs/10.1175/… paper)) use a Runge-Kutta fourth-order scheme with a time-step of 0.05. I've tested this with compiled code and found that RK4 with time-step 0.01 is pretty accurate. odeint can be more accurate, but slower and reducing the solver tolerance doesn't make it much faster, esp when using many calls (as I have to) rather than one large call (as in your code)
@NeillBowler You're welcome! It's a rather cool system and I didn't know about it. Maybe you are interested in scipy.integrate.ode, which lets you choose a RK45 for instance, and doesn't return the whole solution. Instead, you get to call .integrate(time), which may be better suited for whatever you're doing. (PS: there's a small typo in your link, I'm not sure if you'll be able to edit it).

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.