## An improved PARI/gp script for solving the Lorenz ODE system

Today, I’m posting a variant of the PARI/gp script

to solve the Lorenz system of differential equations.

The use of factorials and binomial coefficients has been

eliminated, so I’d expect it to run somewhat faster than

the script from an earlier post.

The lorenzy4.gp PARI/gp script:

Lorenz(X0, Y0, Z0) =

{

order = 300;

fa=vector(order);

xp=vector(order+1);

yp=vector(order+1);

zp=vector(order+1);

dt = 1.0/250 ;

for(J=1,order, fa[J] = factorial(J));

x0 = X0; y0 = Y0; z0 = Z0;

xc=x0;

yc = y0;

zc = z0;

sig = 10.0;

rho = 28.0;

beta = 8.0/3 ;

for(Q = 1,1000,

for(L = 1, 250,

xp[1]=xc;

yp[1] = yc;

zp[1] = zc;

for(J=2, order+1,

xp[J] = (sig*(yp[J-1]-xp[J-1]))/(J-1) ;

yp[J] = (rho*xp[J-1] – sum(K=0,J-2, xp[K+1]*zp[J-1-K]) -yp[J-1])/(J-1);

zp[J] = (sum(K=0,J-2, xp[K+1]*yp[J-1-K]) – beta*zp[J-1])/(J-1); );

xn = xc + sum(J=1, order, xp[J+1]*dt^J);

yn = yc + sum(J=1, order, yp[J+1]*dt^J);

zn = zc + sum(J=1, order, zp[J+1]*dt^J) ;

xc = xn;

yc = yn;

zc = zn; );

print(Q, ” “, precision(xc, 20)))

}

Results for x(100), x(200), x(300) and x(400), which agree with

the Kehlett and Logg reference solution for those values:

? Lorenz(1.0, 0.0, 0.0)

100 12.2860066770189848088523959400428665893

200 -0.428878362621396885476845988868367367420

300 -7.9465362714960512478195581280017205595

400 5.9432455860803396010549930121329814024

Advancing the numerical solution by one time unit,

for example from t = 100.0 to t = 101.0, takes approximately

50 seconds. In the inner loop “for L = 1,250”,

one can see that there is no factorial function and

there are no binomial coefficients. That’s why I consider

this script an improvement over the earlier one.

I believe there is a trade-off that can be made, for the

same precision, between a smaller time-step and a smaller

value of the order of the Taylor series.

So, to minimize the total time, one would look at

pairs of time-step and Taylor series order giving the

same precision , e.g. at t = 1000, or even t = 250,

to find the pair which runs in the least amount of time.

Personal note:

the lorenzy4.gp file is in the directory :

/home/david/graphs/TABU_H_homebrew/CHECK/KEEP6/JULY24

on my PC.