optimization - How to optimize time-stepping algorithm in python? -

i have piece of code increments time-step so-called lorenz95 model (invented ed lorenz in 1995). it's implemented 40-variable model, , displays chaotic behaviour. have coded time-stepping algorithm follows:

class lorenz:     '''lorenz-95 equation'''      global f, dt, size     f = 8     dt = 0.01     size = 40      def __init__(self):         self.x = [random.random() in range(size)]      def euler(self):         '''euler time stepping'''         newvals = [0]*size         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 not slow, unfortunately, code needs make large number of calls it. there way code time-stepping make run faster?

many thanks.

there @ least 2 kinds of possible optimizations: working in smarter way (algorithmic improvements) , working faster.

• in algorithmic side, you're using euler method, first-order method (so global error proportional step size) , has smallish stability region. is, it's not efficient.

• in other side, if you're using standard cpython implementation kind of code going quite slow. around that, try running under pypy. just-in-time compiler can make numerical code run maybe 100x faster. write custom c or cython extension.

but there's better way. solving systems of ordinary differential equations quite common, scipy, 1 of core scientific libraries in python wraps fast, battle-tested fortran libraries solve them. using scipy, both algorithmic improvemnts (as integrators have higher order) , fast implementation.

solving lorenz 95 model set of perturbated initial conditions looks this:

import numpy 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 plt     scipy.integrate import odeint     size = 40     f = 8     t = np.linspace(0, 10, 1001)     x0 = np.random.random(size)     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 output (setting np.random.seed(7), yours can different) nicely chaotic. small perturbations in initial conditions (in 1 of coordinates!) produce different solutions:

but, faster euler time stepping? dt = 0.01 seems 3 times faster, solutions don't match except @ beginning.

if dt reduced, solution provided euler method gets increasingly similar odeint solution, takes longer. notice how smaller dt, later euler solutions loose track of odeint solution. precise euler solution took 600x longer computed solution t=6 odeint t=10. see full script here.

in end, system unstable guess not odeint solution accurate along plotted time.