Essentially I am trying to calculate the Lyapunov exponent the naive way; i.e., using the method described here by Vaggelis_Z.
In each pass of the loop over time a new trajectory is spawned by adding a small displacement vector to the state and running it for a single time step. This is repeated at each time step after renormalisation of the displacement vector.
While invoking the complete NDSolve within the loop works OK, I can see that the costs for a 10+ dimensional system of ODEs will be prohibitive.
I want to speed up things by accessing and changing the state variables and resubmitting them in NDSolve`Iterate
inside the loop. How can I do that?
MarcoB, here's the code I have got
(* original ODE's for the Sprott system *)
f1 = (y1[t]*z1[t] + 0.01);
f2 = x1[t]*x1[t] - y1[t];
f3 = 1 - 4*x1[t] ;
(* cloned ODE's for shadowing the system *)
f4 = (y2[t]*z2[t] + 0.01);
f5 = x2[t]*x2[t] - y2[t];
f6 = 1 - 4*x2[t] ;
(* starting initial conditions for the original trajectory *)
x10 = 1; y10 = 1; z10 = 0;
dx0 = 10^-8;
(* starting initial conditions for the shadowing trajectory *)
x20 = x10 + dx0; y20 = y10; z20 = z10;
tin = 0; tfin = 1000;
dt = .01;
iterat = tfin/dt;
acc = 12;
lexp = {};
sum = 0;
d0 = Sqrt[(x10 - x20)^2 + (y10 - y20)^2 + (z10 - z20)^2];
Timing[For[i = 1, i < iterat, i++,
eqns = {x1'[t] == f1, y1'[t] == f2, z1'[t] == f3, x2'[t] == f4,
y2'[t] == f5, z2'[t] == f6, x1[0] == x10, y1[0] == y10,
z1[0] == z10, x2[0] == x20, y2[0] == y20, z2[0] == z20};
sol = NDSolve[
eqns, {x1[t], y1[t], z1[t], x2[t], y2[t], z2[t]}, {t, 0, dt},
MaxSteps -> Infinity, PrecisionGoal -> acc, AccuracyGoal -> acc];
(* solution functions *)
xx1[t_] = x1[t] /. sol[[1]];
yy1[t_] = y1[t] /. sol[[1]];
zz1[t_] = z1[t] /. sol[[1]];
xx2[t_] = x2[t] /. sol[[1]];
yy2[t_] = y2[t] /. sol[[1]];
zz2[t_] = z2[t] /. sol[[1]];
d1 = Sqrt[(xx1[dt] - xx2[dt])^2 + (yy1[dt] -
yy2[dt])^2 + (zz1[dt] - zz2[dt])^2];
sum += Log[d1/d0];
lyap = sum/(dt*i);
AppendTo[lexp, {dt*i, Log10[lyap]}];
(* components of displacement vector normalised to d0 *)
vc1 = (xx1[dt] - xx2[dt])*(d0/d1);
vc2 = (yy1[dt] - yy2[dt])*(d0/d1);
vc3 = (zz1[dt] - zz2[dt])*(d0/d1);
(*new initial conditions for the original trajectory *)
x10 = xx1[dt];
y10 = yy1[dt];
z10 = zz1[dt];
(* new initial conditions spawned for the shadowing trajectory *)
x20 = x10 + vc1;
y20 = y10 + vc2;
z20 = z10 + vc3;
i = i++;
If[Mod[dt*i, 10] == 0,
Print[" t = ", dt*i, " , ", "largest Lyapunov = ", lyap]]]]
Answer
The run time of the code in the question is 202 seconds on my PC. Even before employing the component structure of NDSolve
, it is helpful to eliminate unnecessary code to both increase speed and increase readability. Doing so yields,
f1 = (y1[t]*z1[t] + 0.01);
f2 = x1[t]*x1[t] - y1[t];
f3 = 1 - 4*x1[t] ;
f4 = (y2[t]*z2[t] + 0.01);
f5 = x2[t]*x2[t] - y2[t];
f6 = 1 - 4*x2[t] ;
x10 = 1; y10 = 1; z10 = 0; dx0 = 10^-8;
x20 = x10 + dx0; y20 = y10; z20 = z10;
tin = 0; tfin = 1000; dt = .01; iterat = Round[tfin/dt]; acc = 12;
eqns = {x1'[t] == f1, y1'[t] == f2, z1'[t] == f3,
x2'[t] == f4, y2'[t] == f5, z2'[t] == f6};
lexp = ConstantArray[0, iterat - 1]; sum = 0;
d0 = Norm[{x10 , y10, z10} - {x20 , y20, z20}];
Timing[Do [
ic = {x1[0] == x10, y1[0] == y10, z1[0] == z10,
x2[0] == x20, y2[0] == y20, z2[0] == z20};
sol = Last /@ Flatten@NDSolve[{eqns, ic}, {x1[dt], y1[dt], z1[dt],
x2[dt], y2[dt], z2[dt]}, {t, 0, dt},
MaxSteps -> Infinity, PrecisionGoal -> acc, AccuracyGoal -> acc];
d1 = Norm[sol[[1 ;; 3]] - sol[[4 ;; 6]]];
sum += Log[d1/d0];
lyap = sum/(dt*i);
lexp[[i]] = {dt*i, Log10[lyap]};
{x10 , y10, z10} = sol[[1 ;; 3]];
{x20 , y20, z20} = {x10 , y10, z10} + (sol[[1 ;; 3]] - sol[[4 ;; 6]]) (d0/d1);
If[Mod[dt*i, 10] == 0,
Print[" t = ", dt*i, " , ", "largest Lyapunov = ", lyap]], {i, iterat - 1}]]
The run time now is only 114 seconds. Although many changes were made, most of the savings were achieved by replacing AppendTo[lexp, {dt*i, Log10[lyap]}]
by lexp[[i]] = {dt*i, Log10[lyap]}
. AppendTo
is notoriously slow for large arrays, here 100 000
elements. Asking NDSolve
to return only the final values, {x1[dt], y1[dt], z1[dt], x2[dt], y2[dt], z2[dt]}
, instead of {x1[t], y1[t], z1[t], x2[t], y2[t], z2[t]}
also helps.
Now, employ the NDSolve
Component Structure, which requires only minor modifications to the preceding code.
f1 = (y1[t]*z1[t] + 0.01);
f2 = x1[t]*x1[t] - y1[t];
f3 = 1 - 4*x1[t] ;
f4 = (y2[t]*z2[t] + 0.01);
f5 = x2[t]*x2[t] - y2[t];
f6 = 1 - 4*x2[t] ;
x10 = 1; y10 = 1; z10 = 0; dx0 = 10^-8;
x20 = x10 + dx0; y20 = y10; z20 = z10;
tin = 0; tfin = 1000; dt = .01; iterat = Round[tfin/dt]; acc = 12;
eqns = {x1'[t] == f1, y1'[t] == f2, z1'[t] == f3,
x2'[t] == f4, y2'[t] == f5, z2'[t] == f6};
ic = {x1[0] == x10, y1[0] == y10, z1[0] == z10, x2[0] == x20, y2[0] == y20, z2[0] == z20};
state = First@NDSolve`ProcessEquations[{eqns, ic},
{x1[dt], y1[dt], z1[dt], x2[dt], y2[dt], z2[dt]}, t ,
MaxSteps -> Infinity, PrecisionGoal -> acc, AccuracyGoal -> acc];
lexp = ConstantArray[0, iterat - 1]; sum = 0;
d0 = Norm[{x10 , y10, z10} - {x20 , y20, z20}];
Timing[Do [
ic = {x1[0] == x10, y1[0] == y10, z1[0] == z10,
x2[0] == x20, y2[0] == y20, z2[0] == z20};
newstate = First@NDSolve`Reinitialize[state, ic];
NDSolve`Iterate[newstate, dt];
sol = Last /@ NDSolve`ProcessSolutions[newstate];
d1 = Norm[sol[[1 ;; 3]] - sol[[4 ;; 6]]];
sum += Log[d1/d0];
lyap = sum/(dt*i);
lexp[[i]] = {dt*i, Log10[lyap]};
{x10 , y10, z10} = sol[[1 ;; 3]];
{x20 , y20, z20} = {x10 , y10, z10} + (sol[[1 ;; 3]] - sol[[4 ;; 6]]) (d0/d1);
If[Mod[dt*i, 10] == 0, Print[" t = ", dt*i, " , ", "largest Lyapunov = ", lyap]],
{i, iterat - 1}]]
Run time now is only 34 seconds. This large additional savings results from moving NDSolve`ProcessEquations
out of the loop
Comments
Post a Comment