Using linear regression of the linearized data between time

Using linear regression of the (linearized) data between time units 5.0 and 15.0, compute the numeric value of the Lyapunov exponent for the plot on slide 15 of Lecture 26. (Note that the rate of separation between two trajectories can be different for different realizations, also the time interval 5.0-15.0 we use here is somewhat subjective. Thus, for any system there is typically a spectrum of Lyapunov exponents, depending on the integration parameters.)

The Lyapunov exponent typically has a unit 1/time, so it is perhaps more intuitive to consider the inverse of the exponent called Lyapunov time or predictability time. It is the characteristic timescale on which a dynamical system becomes unpredictable due to deterministic chaos. As with the exponent, the Lyapunov time might show some variation based on the particular parameters. (In chaos theory studies, many realizations would be run and the maximum exponent or minimum time would be taken as characteristic of the system). Does the Lyapunov time we observe depend on the integrator, or is it rather an intrinsic property of the dynamic system? Test your hypothesis by using the MATLAB integrators ode23 and ode15s instead of ode45 and compare the three values of the Lyapunov time. Note that in the Lorenz equations (as we used them in class) the time is dimensionless, so you should obtain three dimensionless time values (from the inverse exponents). You can use the same interval 5.0- 15.0 for the three linear regression calculations (as a sanity check for all three cases, verify that the logarithm of the separation appears to be growing linearly in that interval). The code from lecture 26 is as followed:

function yp = lorenz_attractor(t,y)

yp = [10 * (y(2)-y(1)); 28*y(1) - y(1)*y(3) - y(2); ...

    y(1)*y(2) - (8/3)*y(3)];

clear

clc

tspan = [0 30]; y0 = [5 5 5]; y00 = [5.00001 5 5];

[t1, y1] = ode45(@lorenz_attractor,tspan,y0);

[t2, y2] = ode45(@lorenz_attractor,tspan,y00);

figure(1)

plot(t1,y1(:,3),t2,y2(:,3),\'--\')

title(\'y_3(t) for y(0)=[5 5 5] (solid) and y(0)=[5.00001 5 5] (dashed)\')

y21 = interp1(t2,y2,t1);

figure(2)

plot(t1,log(abs(y21(:,3)- y1(:,3))))

title(\'logarithm of separation of the two y_3(t) solutions\')

Lecture 26 slide 15 attached:

Butterfly Effect and Chaos logarithm of separation of the two y3(t) solutions -10 -15 exponential growth -20 10 The slope in the growth region of the logarithmic plot defines the so called Lyapunov exponent (it can be determined by linear regression, see Lecture 15.

Solution

The equations of motion in dimensionless variables are

q^° = w

w^° = -sinq + Acos Hw_DtL

Initialize the parameters A = driving amplitude, w_D = driving angular frequency, then set up the initial conditions for the fiducial orbit:

         A = 0.07; wD = 0.75;

q0 = 3.00;

w0 = 0.0;

tend = 2000;

s1 = NDSolve@8q^£@tD ã w@tD, w^£@tD == -Sin@q@tDD + A Sin@wD tD, q@0D ã q0, w@0D ã w0<,

8q, w<, 8t, 0, tend<, MaxSteps Ø 100 000D;

p1 = ParametricPlot@Evaluate@8q@tD, w@tD< ê. s1D, 8t, 0, tend<D

Now set up initial conditions for the nearby comparator orbit:

d0 = 0.001;

q0c = q0 + d0;

w0c = w0;

t = 2.;

kmax = Floor@tend ê tD

term1 = 8<;

1000This choice of \"pull-back\" intervals gives us k_max = 1000 pull-back segments. Set up a loop to go through each pull-back segment, pull back the test orbit, and re-start the calculation. Add each new value of logHd_k ê d₀L into the \"term1\" list for later use:

DoBs2 = NDSolve@8q^£@tD ã w@tD, w^£@tD == -Sin@q@tDD + A Sin@wD tD, q@Hk - 1L tD ã q0c, w@Hk - 1L tD ã w0c<, 8q, w<, 8t, Hk - 1L t, k t<, MaxSteps Ø 100 000D;

H* Solve comparator orbit for kth interval *L q2 = Hq@k tD ê. s2L@@1DD; q1 = Hq@k tD ê. s1L@@1DD; w2 = Hw@k tD ê. s2L@@1DD; w1 = Hw@k tD ê. s1L@@1DD;

d1v = 8q2 - q1, w2 - w1<; H* Difference vector at t=kt *L

d1 = Norm@d1vD ê d0;            H* Difference magnitude ratio with d0 *L

q0c = q1 +

q2 - q1

; w0c = w1 +

w2 - w1

; H* Pull back q and w *L

d1

d1

AppendTo@term1, Log@d1DD, 8k, 1, kmax<F

Form a list of partial sums l_k for k = 1, ..., k_max - 1

^k      term1@@jDD

          lk = TableB       ^j=1                           , 8k, 1, kmax - 1<F;

k t

Make a quick plot of the l_k \'s:

          ListPlot@lk, PlotRange Ø All, Frame Ø True,

FrameLabel Ø 8\"k\", \"l_k\"<, PlotRange Ø 80, 0.85<D

k

This looks it levels off to a nearly constant value after about t = 400. Let\'s check by expanding the axes out with a log-log plot:

          ListLogLogPlot@lk, PlotRange Ø 88100, 1100<, Automatic<,

Frame Ø True, FrameLabel Ø 8k, l_k<, GridLines Ø AutomaticD

k

Here l_k decreases, then clearly levels off for large k to a non-zero value. You can estimate the asymptotic value

by eyeball at about l_LCE º 0.085. To get a more accurate value, take the mean of the last 200 values and use the standard deviation as the uncertainty:

          Framed@Row@8\"l_LCE= \", Round@Mean@lk@@800 ;; Length@lkDDDD, 0.001D,

\"±\", Round@StandardDeviation@lk@@800 ;; Length@lkDDDD, 0.0001D<DD

             l_LCE= 0.087±0.0014

Initialize the parameters as follows for a regular orbit:

w0 = 0.0; tend = 1000;

s1 = NDSolve@8q^£@tD ã w@tD, w^£@tD == -Sin@q@tDD + A Sin@wD tD, q@0D ã q0, w@0D ã w0<,

8q, w<, 8t, 0, tend<, MaxSteps Ø 100 000D;

p1 = ParametricPlot@Evaluate@8q@tD, w@tD< ê. s1D, 8t, 0, tend<, Frame Ø True, FrameLabel Ø 8q, w<, ImageSize Ø SmallD

Here we have a bounded, oscillatory orbit (quasi-periodic). Let\'s set up the comparator orbit and do the

Benettin algorithm for it:

d0 = 0.001;

q0c = q0 + d0;

w0c = 0.0;

t = 2.;

kmax = Floor@tend ê tD

term1 = 8<;

ListPlot@ln, PlotRange Ø All, Frame Ø True, FrameLabel Ø 8k, l_k<D

k

          ListLogLogPlot@ln, PlotRange Ø 8810, 500<, Automatic<,

Frame Ø True, FrameLabel Ø 8k, l_k<, GridLines Ø AutomaticD

Here l_k decreases even more rapidly than power law and it never reaches an asymptotically constant value.

q0c = q1 +	q2 - q1	; w0c = w1 +	w2 - w1	; H* Pull back q and w *L
	d1	d1