meditationatae

Just another WordPress.com site

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.

Advertisements

Written by meditationatae

August 11, 2017 at 3:53 pm

Posted in History

%d bloggers like this: