Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Nov 2016
Lecture topic: Dynamical Systems

Runge Kutta for Systems (Begininnings of Dynamical Systems)

How to Extend RK4 to Higher Dimensions

It is relatively easy to extend the numerical routines from a single unknown funtion to multiple unknown functions. All we have to do is label the pieces to keep track of which function goes with which variable.

For example, to solve

y'(t)=f(t,y)=y-(t^2)+1, 0≤t≤4, y(0)=0.5

we can use Runge-Kutta Order 4 (RK4):

f[t_,y_] = y-(t^2)+1;

h = 0.2;
num=20;

t=0.0;
w=0.5;

data={{t,w}};
Do[
k1=h*f[t,w];
k2=h*f[t+h/2,w+k1/2];
k3=h*f[t+h/2,w+k2/2];
k4=h*f[t+h,w+k3];
w=w+(k1+2k2+2k3+k4)/6.0;
t=t+h;
data=Append[data,{t,w}],
{i,1,num}]

TableForm[data,TableHeadings->{None,{"(Subscript[t, i])","(Subscript[y, i])"}}]

To extend this process to N unknown functions, we will just need to keep track of (Subscript[w, 1]),(Subscript[w, 2]),...(Subscript[w, N]) which means we will need (Subscript[k1, 1]),(Subscript[k1, 2]),...(Subscript[k1, N]) and (Subscript[k2, 1]),(Subscript[k2, 2]),...(Subscript[k2, N]) and (Subscript[k3, 1]),(Subscript[k3, 2]),...(Subscript[k3, N]) and (Subscript[k4, 1]),(Subscript[k4, 2]),...(Subscript[k4, N]). Let’s start by looking at an N=2 predator-prey situation to see how the works. It’s actually not too complicated.

Predator-Prey Systems

Consider the nonlinear system initial value problem (this is the form of typical predator-prey systems, where two species compete for resources in the same area, and the population of each species is given by (Subscript[y, i]))

(Subscript[y, 1])'(t)=(Subscript[f, 1])(t,(Subscript[y, 1]),(Subscript[y, 2]))=(Subscript[y, 1])-(1/2)(Subscript[y, 1])(Subscript[y, 2]),
 (Subscript[y, 2])'(t)=(Subscript[f, 2])(t,(Subscript[y, 1]),(Subscript[y, 2]))=-(3/4)(Subscript[y, 2])+(1/4)(Subscript[y, 1])(Subscript[y, 2]), 
0≤t≤8, (Subscript[y, 1])(0)=(Subscript[α, 1])=2, (Subscript[y, 2])(0)=(Subscript[α, 2])=3

We can use Runga Kutta to solve this system!

Mathematica Solution

First, let's try to solve this using Mathematica:

Clear[f,t,y]
(Subscript[f, 1])[t_,x_,y_]=x-0.5x*y;
(Subscript[f, 2])[t_,x_,y_]=-0.75y+0.25x*y;
DSolve[{x'[t]==(Subscript[f, 1])[t,x[t],y[t]],y'[t]==(Subscript[f, 2])[t,x[t],y[t]],x[0]==2,y[0]==3},{x[t],y[t]},t]

Couldn't do it.

Runge Kutta Method (RK4)

Now, let's solve this using Runge Kutta Order Four. Notice that RK4 has a nice, easy to understand structure when extended to N=2.

Clear[t,y,f]
(Subscript[f, 1])[t_,x_,y_]=x-0.5x*y;
(Subscript[f, 2])[t_,x_,y_]=-0.75y+0.25x*y;

h = 0.1;
num=80;

t=0.0;
(Subscript[w, 1])=2.0;
(Subscript[w, 2])=3.0;

data={{t,(Subscript[w, 1]),(Subscript[w, 2])}};
Do[
(Subscript[k1, 1])=h*(Subscript[f, 1])[t,(Subscript[w, 1]),(Subscript[w, 2])];
(Subscript[k1, 2])=h*(Subscript[f, 2])[t,(Subscript[w, 1]),(Subscript[w, 2])];
(Subscript[k2, 1])=h*(Subscript[f, 1])[t+h/2,(Subscript[w, 1])+(Subscript[k1, 1])/2,(Subscript[w, 2])+(Subscript[k1, 2])/2];
(Subscript[k2, 2])=h*(Subscript[f, 2])[t+h/2,(Subscript[w, 1])+(Subscript[k1, 1])/2,(Subscript[w, 2])+(Subscript[k1, 2])/2];
(Subscript[k3, 1])=h*(Subscript[f, 1])[t+h/2,(Subscript[w, 1])+(Subscript[k2, 1])/2,(Subscript[w, 2])+(Subscript[k2, 2])/2];
(Subscript[k3, 2])=h*(Subscript[f, 2])[t+h/2,(Subscript[w, 1])+(Subscript[k2, 1])/2,(Subscript[w, 2])+(Subscript[k2, 2])/2];
(Subscript[k4, 1])=h*(Subscript[f, 1])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2])];
(Subscript[k4, 2])=h*(Subscript[f, 2])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2])];
(Subscript[w, 1])=(Subscript[w, 1])+((Subscript[k1, 1])+2(Subscript[k2, 1])+2(Subscript[k3, 1])+(Subscript[k4, 1]))/6.0;
(Subscript[w, 2])=(Subscript[w, 2])+((Subscript[k1, 2])+2(Subscript[k2, 2])+2(Subscript[k3, 2])+(Subscript[k4, 2]))/6.0;
t=t+h;
data=Append[data,{t,(Subscript[w, 1]),(Subscript[w, 2])}],
{i,1,num}]

TableForm[data,TableHeadings->{None,{"(Subscript[t, i])","((Subscript[y, 1]))_i","((Subscript[y, 2]))_i"}}]

No problem. Let's graph the solution:

phase=Table[{data[[i,2]],data[[i,3]]},{i,1,Length[data]}];

rk4plot = ListPlot[phase,Frame->True,FrameLabel->{"(Subscript[y, 1])","(Subscript[y, 2])"},PlotStyle->{Blue,PointSize[0.02]},ImageSize->{500,Automatic},BaseStyle->20]

rk4plot=ListPlot[phase,Joined->True,Frame->True,FrameLabel->{"(Subscript[y, 1])","(Subscript[y, 2])"},PlotStyle->{Blue},ImageSize->{500,Automatic},BaseStyle->20]

The populations depend on each other (obviously), but neither population is able to increase without bound.

Multiple Initial Points

Here we slap a loop around the RK4 to get a bunch of solutions for different initial values, using the Corput sequence to sample the initial starting values.

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[N[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))]];
]

Clear[t,y,f]
(Subscript[f, 1])[t_,x_,y_]=x-0.5 x y;
(Subscript[f, 2])[t_,x_,y_]=-0.75y+0.25x y;
h=0.1;
num=100;
numicond=10;
plots={};
Do[t=0.;
xa=0;
xb=8;
ya=0;
yb=5;
(Subscript[w, 1])=N[xa+(xb-xa) Corput[ic,2]];
(Subscript[w, 2])=N[ya+(yb-ya) Corput[ic,3]];
data={{t,(Subscript[w, 1]),(Subscript[w, 2])}};
Do[
(Subscript[k1, 1])=h (Subscript[f, 1])[t,(Subscript[w, 1]),(Subscript[w, 2])];
(Subscript[k1, 2])=h (Subscript[f, 2])[t,(Subscript[w, 1]),(Subscript[w, 2])];
(Subscript[k2, 1])=h (Subscript[f, 1])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2)];
(Subscript[k2, 2])=h (Subscript[f, 2])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2)];
(Subscript[k3, 1])=h (Subscript[f, 1])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2)];
(Subscript[k3, 2])=h (Subscript[f, 2])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2)];
(Subscript[k4, 1])=h (Subscript[f, 1])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2])];
(Subscript[k4, 2])=h (Subscript[f, 2])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2])];
(Subscript[w, 1])=(Subscript[w, 1])+(((Subscript[k1, 1])+2 (Subscript[k2, 1])+2 (Subscript[k3, 1])+(Subscript[k4, 1]))/6.);
(Subscript[w, 2])=(Subscript[w, 2])+(((Subscript[k1, 2])+2 (Subscript[k2, 2])+2 (Subscript[k3, 2])+(Subscript[k4, 2]))/6.);
t=t+h;
data=Append[data,{t,(Subscript[w, 1]),(Subscript[w, 2])}],{i,1,num}];
 data;

phase=Table[{data[[i,2]],data[[i,3]]},{i,1,Length[data]}];

rk4plot=ListPlot[phase,Joined->True,Frame->True,FrameLabel->{"(Subscript[y, 1])","(Subscript[y, 2])"},PlotStyle->{Blue},ImageSize->{500,Automatic}];

plots=Append[plots,rk4plot];
,{ic,1,numicond}]

Show[plots,PlotRange->{All,All},Axes->None,BaseStyle->20]

This is the sketch one usually sees when looking at predator-prey population models. Notice that there is a bit of blockiness in some of the lines, since we are just connecting points to get the plot.

New Example: Limit Cycles

Here we slap a loop around the RK4 to get a bunch of solutions for different initial values, using the Corput sequence to sample the initial starting values. For this example (simply a mathematical example, it isn't a predator-prey example anymore), we see that the solutions all seem to be approaching a curve, which is known as a limit cycle. A limit cycle is an isolated closed trajectory, where isolated means that neighboring trajectories are not closed--they spiral either toward or away from the limit cycle.

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[N[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))]];
]

Clear[t,y,f]
(Subscript[f, 1])[t_,x_,y_]=x+y-x ((x^2)+2 (y^2));
(Subscript[f, 2])[t_,x_,y_]=-x+y-y ((x^2)+(y^2));
h=0.1;
num=100;
numicond=10;
plots={};
Do[t=0.;
xa=-2;
xb=2;
ya=-2;
yb=2;
(Subscript[w, 1])=N[xa+(xb-xa) Corput[ic,2]];
(Subscript[w, 2])=N[ya+(yb-ya) Corput[ic,3]];
data={{t,(Subscript[w, 1]),(Subscript[w, 2])}};
Do[(Subscript[k1, 1])=h (Subscript[f, 1])[t,(Subscript[w, 1]),(Subscript[w, 2])];
(Subscript[k1, 2])=h (Subscript[f, 2])[t,(Subscript[w, 1]),(Subscript[w, 2])];
(Subscript[k2, 1])=h (Subscript[f, 1])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2)];
(Subscript[k2, 2])=h (Subscript[f, 2])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2)];
(Subscript[k3, 1])=h (Subscript[f, 1])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2)];
(Subscript[k3, 2])=h (Subscript[f, 2])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2)];
(Subscript[k4, 1])=h (Subscript[f, 1])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2])];
(Subscript[k4, 2])=h (Subscript[f, 2])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2])];
(Subscript[w, 1])=(Subscript[w, 1])+(((Subscript[k1, 1])+2 (Subscript[k2, 1])+2 (Subscript[k3, 1])+(Subscript[k4, 1]))/6.);
(Subscript[w, 2])=(Subscript[w, 2])+(((Subscript[k1, 2])+2 (Subscript[k2, 2])+2 (Subscript[k3, 2])+(Subscript[k4, 2]))/6.);
t=t+h;
data=Append[data,{t,(Subscript[w, 1]),(Subscript[w, 2])}],{i,1,num}] 
data;

phase=Table[{data[[i,2]],data[[i,3]]},{i,1,Length[data]}];

rk4plot=ListPlot[phase,Joined->True,Frame->True,FrameLabel->{"(Subscript[y, 1])","(Subscript[y, 2])"},PlotStyle->{Blue},ImageSize->{500,Automatic},DisplayFunction->Identity];

plots=Append[plots,rk4plot];
,{ic,1,numicond}]

Show[plots,PlotRange->{All,All},Axes->None,BaseStyle->20]

In this case, the solution is asymptotically stable, since solutions on both sides of the limit cycle approach the limit cycle. “Asymptotic” just means in the limit as t->∞.

The Lorenz Equations: Fluid Dynamics and The Butterfly Curve

The Lorenz equations are used to model the motion of the earth's atmosphere (it is basically a fluid dynamics problem). Let's take a look at the solution to these equations, which are in three dimensions! So here we have N=3 and have to modify our RK4 again.

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[N[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))]];
]

Clear[t,x,y,z,f]
(Subscript[f, 1])[t_,x_,y_,z_]=10. (-x+y);
(Subscript[f, 2])[t_,x_,y_,z_]=28 x-y-x z;
(Subscript[f, 3])[t_,x_,y_,z_]=-((8 z)/3)+x y;
h=0.01;
num=2000;
t=0.;
(Subscript[w, 1])=5.;
(Subscript[w, 2])=5.;
(Subscript[w, 3])=5.;
data={{t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])}};
Do[(Subscript[k1, 1])=h (Subscript[f, 1])[t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])];
(Subscript[k1, 2])=h (Subscript[f, 2])[t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])];
(Subscript[k1, 3])=h (Subscript[f, 3])[t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])];
(Subscript[k2, 1])=h (Subscript[f, 1])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2),(Subscript[w, 3])+((Subscript[k1, 3])/2)];(Subscript[k2, 2])=h (Subscript[f, 2])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2),(Subscript[w, 3])+((Subscript[k1, 3])/2)];(Subscript[k2, 3])=h (Subscript[f, 3])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2),(Subscript[w, 3])+((Subscript[k1, 3])/2)];(Subscript[k3, 1])=h (Subscript[f, 1])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2),(Subscript[w, 3])+((Subscript[k2, 3])/2)];(Subscript[k3, 2])=h (Subscript[f, 2])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2),(Subscript[w, 3])+((Subscript[k2, 3])/2)];(Subscript[k3, 3])=h (Subscript[f, 3])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2),(Subscript[w, 3])+((Subscript[k2, 3])/2)];(Subscript[k4, 1])=h (Subscript[f, 1])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2]),(Subscript[w, 3])+(Subscript[k3, 3])];(Subscript[k4, 2])=h (Subscript[f, 2])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2]),(Subscript[w, 3])+(Subscript[k3, 3])];(Subscript[k4, 3])=h (Subscript[f, 3])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2]),(Subscript[w, 3])+(Subscript[k3, 3])];
(Subscript[w, 1])=(Subscript[w, 1])+(((Subscript[k1, 1])+2 (Subscript[k2, 1])+2 (Subscript[k3, 1])+(Subscript[k4, 1]))/6.);
(Subscript[w, 2])=(Subscript[w, 2])+(((Subscript[k1, 2])+2 (Subscript[k2, 2])+2 (Subscript[k3, 2])+(Subscript[k4, 2]))/6.);
(Subscript[w, 3])=(Subscript[w, 3])+(((Subscript[k1, 3])+2 (Subscript[k2, 3])+2 (Subscript[k3, 3])+(Subscript[k4, 3]))/6.);
t=t+h;
data=Append[data,{t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])}],{i,1,num}]

phasexy=Table[{data[[i,2]],data[[i,3]]},{i,1,Length[data]}];
rk4plot=ListPlot[phasexy,Joined->True,Frame->True,FrameLabel->{"(Subscript[y, 1])","(Subscript[y, 2])"},PlotStyle->{Blue},ImageSize->{500,Automatic},BaseStyle->20]

phasexz=Table[{data[[i,2]],data[[i,4]]},{i,1,Length[data]}];
rk4plot=ListPlot[phasexz,Joined->True,Frame->True,FrameLabel->{"(Subscript[y, 1])","(Subscript[y, 2])"},PlotStyle->{Blue},ImageSize->{500,Automatic},BaseStyle->20]

This is the famous "butterfly curve", which is also important in chaos theory. It is very pretty, don't you think?

If we wanted to plot it in three dimensions, we can.

phasexyz=Table[{data[[i,2]],data[[i,3]],data[[i,4]]},{i,1,Length[data]}];
rk4plot=ListPointPlot3D[phasexyz,PlotStyle->{Blue},ImageSize->{500,Automatic},BaseStyle->20]

The Lorenz Equations: Chaos

Let's look at the solution to the Lorenz equations for two trajectories which start very close together.

Clear[t,x,y,z,f]
(Subscript[f, 1])[t_,x_,y_,z_]=10.(-x+y);
(Subscript[f, 2])[t_,x_,y_,z_]=28 x-y-x z;
(Subscript[f, 3])[t_,x_,y_,z_]=-((8 z)/3)+x y;
h=0.01;
num=2000;
t=0.;
(Subscript[w, 1])=5.;
(Subscript[w, 2])=5.;
(Subscript[w, 3])=5.;
data={{t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])}};
Do[(Subscript[k1, 1])=h (Subscript[f, 1])[t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])];
(Subscript[k1, 2])=h (Subscript[f, 2])[t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])];
(Subscript[k1, 3])=h (Subscript[f, 3])[t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])];
(Subscript[k2, 1])=h (Subscript[f, 1])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2),(Subscript[w, 3])+((Subscript[k1, 3])/2)];(Subscript[k2, 2])=h (Subscript[f, 2])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2),(Subscript[w, 3])+((Subscript[k1, 3])/2)];(Subscript[k2, 3])=h (Subscript[f, 3])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2),(Subscript[w, 3])+((Subscript[k1, 3])/2)];(Subscript[k3, 1])=h (Subscript[f, 1])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2),(Subscript[w, 3])+((Subscript[k2, 3])/2)];(Subscript[k3, 2])=h (Subscript[f, 2])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2),(Subscript[w, 3])+((Subscript[k2, 3])/2)];(Subscript[k3, 3])=h (Subscript[f, 3])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2),(Subscript[w, 3])+((Subscript[k2, 3])/2)];(Subscript[k4, 1])=h (Subscript[f, 1])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2]),(Subscript[w, 3])+(Subscript[k3, 3])];(Subscript[k4, 2])=h (Subscript[f, 2])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2]),(Subscript[w, 3])+(Subscript[k3, 3])];(Subscript[k4, 3])=h (Subscript[f, 3])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2]),(Subscript[w, 3])+(Subscript[k3, 3])];
(Subscript[w, 1])=(Subscript[w, 1])+(((Subscript[k1, 1])+2 (Subscript[k2, 1])+2 (Subscript[k3, 1])+(Subscript[k4, 1]))/6.);
(Subscript[w, 2])=(Subscript[w, 2])+(((Subscript[k1, 2])+2 (Subscript[k2, 2])+2 (Subscript[k3, 2])+(Subscript[k4, 2]))/6.);
(Subscript[w, 3])=(Subscript[w, 3])+(((Subscript[k1, 3])+2 (Subscript[k2, 3])+2 (Subscript[k3, 3])+(Subscript[k4, 3]))/6.);t=t+h;data=Append[data,{t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])}],{i,1,num}]

phasexy=Table[{data[[i,2]],data[[i,3]]},{i,1,Length[data]}];
rk4plot=ListPlot[phasexy,Joined->True,Frame->True,FrameLabel->{"(Subscript[y, 1])","(Subscript[y, 2])"},PlotStyle->{Blue},ImageSize->{500,Automatic}]

phasexz=Table[{data[[i,2]],data[[i,4]]},{i,1,Length[data]}];
rk4plot=ListPlot[phasexz,Joined->True,Frame->True,FrameLabel->{"(Subscript[y, 1])","(Subscript[y, 2])"},PlotStyle->{Blue},ImageSize->{500,Automatic}]

phasexyz=Table[{data[[i,2]],data[[i,3]],data[[i,4]]},{i,1,Length[data]}];

The only change from the above is in the initial onditions (highlighted in red below):

Clear[t,x,y,z,f]
(Subscript[f, 1])[t_,x_,y_,z_]=10.(-x+y);
(Subscript[f, 2])[t_,x_,y_,z_]=28 x-y-x z;
(Subscript[f, 3])[t_,x_,y_,z_]=-((8 z)/3)+x y;
h=0.01;
num=2000;
t=0.;
(Subscript[w, 1])=5.-0.1;
(Subscript[w, 2])=5.;
(Subscript[w, 3])=5.;
data={{t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])}};
Do[(Subscript[k1, 1])=h (Subscript[f, 1])[t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])];
(Subscript[k1, 2])=h (Subscript[f, 2])[t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])];
(Subscript[k1, 3])=h (Subscript[f, 3])[t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])];
(Subscript[k2, 1])=h (Subscript[f, 1])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2),(Subscript[w, 3])+((Subscript[k1, 3])/2)];(Subscript[k2, 2])=h (Subscript[f, 2])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2),(Subscript[w, 3])+((Subscript[k1, 3])/2)];(Subscript[k2, 3])=h (Subscript[f, 3])[t+(h/2),(Subscript[w, 1])+((Subscript[k1, 1])/2),(Subscript[w, 2])+((Subscript[k1, 2])/2),(Subscript[w, 3])+((Subscript[k1, 3])/2)];(Subscript[k3, 1])=h (Subscript[f, 1])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2),(Subscript[w, 3])+((Subscript[k2, 3])/2)];(Subscript[k3, 2])=h (Subscript[f, 2])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2),(Subscript[w, 3])+((Subscript[k2, 3])/2)];(Subscript[k3, 3])=h (Subscript[f, 3])[t+(h/2),(Subscript[w, 1])+((Subscript[k2, 1])/2),(Subscript[w, 2])+((Subscript[k2, 2])/2),(Subscript[w, 3])+((Subscript[k2, 3])/2)];(Subscript[k4, 1])=h (Subscript[f, 1])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2]),(Subscript[w, 3])+(Subscript[k3, 3])];(Subscript[k4, 2])=h (Subscript[f, 2])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2]),(Subscript[w, 3])+(Subscript[k3, 3])];(Subscript[k4, 3])=h (Subscript[f, 3])[t+h,(Subscript[w, 1])+(Subscript[k3, 1]),(Subscript[w, 2])+(Subscript[k3, 2]),(Subscript[w, 3])+(Subscript[k3, 3])];
(Subscript[w, 1])=(Subscript[w, 1])+(((Subscript[k1, 1])+2 (Subscript[k2, 1])+2 (Subscript[k3, 1])+(Subscript[k4, 1]))/6.);
(Subscript[w, 2])=(Subscript[w, 2])+(((Subscript[k1, 2])+2 (Subscript[k2, 2])+2 (Subscript[k3, 2])+(Subscript[k4, 2]))/6.);
(Subscript[w, 3])=(Subscript[w, 3])+(((Subscript[k1, 3])+2 (Subscript[k2, 3])+2 (Subscript[k3, 3])+(Subscript[k4, 3]))/6.);t=t+h;data=Append[data,{t,(Subscript[w, 1]),(Subscript[w, 2]),(Subscript[w, 3])}],{i,1,num}]

phasexy2=Table[{data[[i,2]],data[[i,3]]},{i,1,Length[data]}];
rk4plot2=ListPlot[phasexy2,Joined->True,Frame->True,FrameLabel->{"(Subscript[y, 1])","(Subscript[y, 2])"},PlotStyle->{Red},ImageSize->{500,Automatic}]

phasexz2=Table[{data[[i,2]],data[[i,4]]},{i,1,Length[data]}];
rk4plot2=ListPlot[phasexz2,Joined->True,Frame->True,FrameLabel->{"(Subscript[y, 1])","(Subscript[y, 2])"},PlotStyle->{Red},ImageSize->{500,Automatic}]

phasexyz2=Table[{data[[i,2]],data[[i,3]],data[[i,4]]},{i,1,Length[data]}];

Show[rk4plot,rk4plot2]

We can see that the two curves are starting to be in different regions after some time, but it is a bit hard to see the divergence of the two trajectories. Let’s make the divergence more apparent through an animation of points following the two curves. You might have to slow this animation down to make it easier to see. It is a beautiful dance, and shows the concept of chaos, or sensititivity to initial conditions. A small change in initial conditions can lead to a large change in behaviour at later times. This is known as the “butterfly effect”: you might get a hurricane or sunshine depending on whether a buterfly had flapped its wings several weeks before.

point1[i_]:=Graphics[{Blue,Disk[phasexz[[i]],0.5]}]
point2[i_]:=Graphics[{Red,Disk[phasexz2[[i]],0.5]}]

Manipulate[Show[point1[i],point2[i],ImageSize->{400,Automatic},PlotRange->{{-20,20},{0,50}},AspectRatio->50/40,Frame->True],{i,1,num,3}]

Here’s a plot in three dimensions.

point1b[i_]:=Graphics3D[{Blue,Sphere[phasexyz[[i]],1.0]}]
point2b[i_]:=Graphics3D[{Red,Sphere[phasexyz2[[i]],1.0]}]

Manipulate[Show[point1b[i],point2b[i],ImageSize->{400,Automatic},PlotRange->{{-20,20},{-20,20},{0,50}},AspectRatio->50/40,Frame->True],{i,1,num,3}]

We can quantify the distance between the two points using the distance formula, and get a nice static plot that represents the chaos in the system:

dist=Table[√((((phasexyz[[i,1]]-phasexyz2[[i,1]]))^2)+(((phasexyz[[i,2]]-phasexyz2[[i,2]]))^2)+(((phasexyz[[i,3]]-phasexyz2[[i,3]]))^2)),{i,1,num}];

ListPlot[dist,BaseStyle->20,AxesLabel->{"n","Distance"}]

If we want to discuss the amount of chaos in the system in more concrete terms, we turn our attention to Lyapunov exponents.

Lyapunov Exponents

I love Lyapunov (sometimes spelled Liapunov) exponents! Don’t worry about all the details, but I want you to see an application of systems and this is a nice one. Lyapunov exponents attempt to measure the exponential divergence of nearby trajectories. Since we are attempting to measure sensitivity to initial conditions, we can compute quantities like:

(( ∂x(t))/(∂y(0)))

which will measure the rate of change of the coordinate x with respect to the initial point y(0) as a function of time t. Of course, we have a bunch of these quantities, which are usually written in terms of a Jacobi Matrix (here is the Jacobi matrix for N=3)

J=(({{(( ∂x(t))/(∂x(0))), (( ∂x(t))/(∂y(0))), (( ∂x(t))/(∂z(0)))}, {(( ∂y(t))/(∂x(0))), (( ∂y(t))/(∂y(0))), (( ∂y(t))/(∂z(0)))}, {(( ∂z(t))/(∂x(0))), (( ∂z(t))/(∂y(0))), (( ∂z(t))/(∂z(0)))}}))

Each element of the Jacobin matrix can be thought of an a new unknown function to be determined. This is especially important since if we are hoping to solve a system numerically, knowing the Lyapunov exponent can help us gauge how numerical errors will propogate through the system (since sensitivity to initial conditions can be thought of as round-off error entering the approximation to the solution). There will be a time horizon beyond which prediction breaks down.

So, the set of differential equations to be solved has expanded from 3 to 3+(3^2)=12 if we want to compute the Jacobi matrix, where we will solve for the variables:

(Subscript[x, 1])=x
(Subscript[x, 2])=y
(Subscript[x, 3])=z
(Subscript[x, 4])=(( ∂x(t))/(∂x(0)))
(Subscript[x, 5])=(( ∂y(t))/(∂x(0)))
(Subscript[x, 6])=(( ∂z(t))/(∂x(0)))
(Subscript[x, 7])=(( ∂x(t))/(∂y(0)))
(Subscript[x, 8])=(( ∂y(t))/(∂y(0)))
(Subscript[x, 9])=(( ∂z(t))/(∂y(0)))
(Subscript[x, 10])=(( ∂x(t))/(∂z(0)))
(Subscript[x, 11])=(( ∂y(t))/(∂z(0)))
(Subscript[x, 12])=(( ∂z(t))/(∂z(0)))

We get equations for these extra variables by differentiating the differential equations we started with with respect to the intial conditions. For example,

(dx/dt)=10 (-x+y)
(∂/(∂x(0)))(dx/dt)=(∂/(∂x(0)))10 (-x+y)
(d/dt)(((∂x(t))/(∂x(0))))=10(-((∂x(t))/(∂x(0)))+((∂y(t))/(∂x(0))))

Using the new notation, these equations become
((Subscript[dx, 1])/dt)=(Subscript[f, 1])=10 (-(Subscript[x, 1])+(Subscript[x, 2]))
((Subscript[dx, 4])/dt)=(Subscript[f, 4])=10(-(Subscript[x, 4])+(Subscript[x, 5])), with initial condition (Subscript[x, 4])(0)=1

OK, now we need to upgrade our RK4 routine to be more flexible. Here we have restructured the Runge-Kutta code to be more flexible with the number of equations we have to solve. Let’s check it before putting in the Jacobi matrix:

Clear[t,x,y,z,τ,f]
(Subscript[f, 1])=10 (-(Subscript[x, 1])+(Subscript[x, 2]));
(Subscript[f, 2])=28 (Subscript[x, 1])-(Subscript[x, 2])-(Subscript[x, 1]) (Subscript[x, 3]);
(Subscript[f, 3])=-((8 (Subscript[x, 3]))/3)+(Subscript[x, 1]) (Subscript[x, 2]);
h=0.01;
num=2000;
neq=3;
τ=0.;
(Subscript[w, 1])=5.;
(Subscript[w, 2])=5.;
(Subscript[w, 3])=5.;
data={};
datapoint={τ};
Do[datapoint=Append[datapoint,(Subscript[w, i])],{i,1,neq}];
data=Append[data,datapoint];
Do[subs={t->τ};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])],{i,1,neq}];Do[(Subscript[k1, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+(h/2)};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+((Subscript[k1, i])/2)],{i,1,neq}];Do[(Subscript[k2, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+(h/2)};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+((Subscript[k2, i])/2)],{i,1,neq}];Do[(Subscript[k3, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+h};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+(Subscript[k3, i])],{i,1,neq}];Do[(Subscript[k4, i])=h (Subscript[f, i])/.subs,{i,1,neq}];Do[(Subscript[w, i])=(Subscript[w, i])+(((Subscript[k1, i])+2 (Subscript[k2, i])+2 (Subscript[k3, i])+(Subscript[k4, i]))/6.),{i,1,neq}];
τ=τ+h;
datapoint={τ};
Do[datapoint=Append[datapoint,(Subscript[w, i])],{i,1,neq}];data=Append[data,datapoint];
,{i,1,num}]
phasexz=Table[{data[[i,2]],data[[i,4]]},{i,1,Length[data]}];
rk4plot=ListPlot[phasexz,Joined->True,Frame->True,FrameLabel->{"(Subscript[y, 1])","(Subscript[y, 3])"},PlotStyle->{Blue},ImageSize->{500,Automatic},BaseStyle->20]

If we want to calculate the Jacobi Matrix, we need to add derivatives with respect to the initial conditions:

Clear[t,x,y,z,τ,f]
(Subscript[f, 1])=10 (-(Subscript[x, 1])+(Subscript[x, 2]));
(Subscript[f, 2])=28 (Subscript[x, 1])-(Subscript[x, 2])-(Subscript[x, 1]) (Subscript[x, 3]);
(Subscript[f, 3])=-((8 (Subscript[x, 3]))/3)+(Subscript[x, 1]) (Subscript[x, 2]);
(Subscript[f, 4])=10 (-(Subscript[x, 4])+(Subscript[x, 5]));
(Subscript[f, 5])=10 (-(Subscript[x, 7])+(Subscript[x, 8]));
(Subscript[f, 6])=10 (-(Subscript[x, 10])+(Subscript[x, 11]));
(Subscript[f, 7])=28 (Subscript[x, 4])-(Subscript[x, 5])-(Subscript[x, 4]) (Subscript[x, 3])-(Subscript[x, 1]) (Subscript[x, 6]);
(Subscript[f, 8])=28 (Subscript[x, 7])-(Subscript[x, 8])-(Subscript[x, 7]) (Subscript[x, 3])-(Subscript[x, 1]) (Subscript[x, 9]);
(Subscript[f, 9])=28 (Subscript[x, 10])-(Subscript[x, 11])-(Subscript[x, 10]) (Subscript[x, 3])-(Subscript[x, 1]) (Subscript[x, 12]);
(Subscript[f, 10])=-((8 (Subscript[x, 6]))/3)+(Subscript[x, 4]) (Subscript[x, 1])+(Subscript[x, 2]) (Subscript[x, 5]);
(Subscript[f, 11])=-((8 (Subscript[x, 9]))/3)+(Subscript[x, 7]) (Subscript[x, 1])+(Subscript[x, 2]) (Subscript[x, 8]);
(Subscript[f, 12])=-((8 (Subscript[x, 12]))/3)+(Subscript[x, 10]) (Subscript[x, 1])+(Subscript[x, 2]) (Subscript[x, 11]);
h=0.01;
num=2000;
neq=12;
τ=0.;
(Subscript[w, 1])=5.;
(Subscript[w, 2])=5.;
(Subscript[w, 3])=5.;
(Subscript[w, 4])=1.;
(Subscript[w, 5])=0.;
(Subscript[w, 6])=0.;
(Subscript[w, 7])=0.;
(Subscript[w, 8])=1.;
(Subscript[w, 9])=0.;
(Subscript[w, 10])=0.;
(Subscript[w, 11])=0.;
(Subscript[w, 12])=1.;
data={};
datapoint={τ};
Do[datapoint=Append[datapoint,(Subscript[w, i])],{i,1,neq}];
data=Append[data,datapoint];
Do[subs={t->τ};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])],{i,1,neq}];Do[(Subscript[k1, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+(h/2)};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+((Subscript[k1, i])/2)],{i,1,neq}];Do[(Subscript[k2, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+(h/2)};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+((Subscript[k2, i])/2)],{i,1,neq}];Do[(Subscript[k3, i])=h (Subscript[f, i])/.subs,{i,1,neq}];subs={t->τ+h};Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+(Subscript[k3, i])],{i,1,neq}];Do[(Subscript[k4, i])=h (Subscript[f, i])/.subs,{i,1,neq}];Do[(Subscript[w, i])=(Subscript[w, i])+(((Subscript[k1, i])+2 (Subscript[k2, i])+2 (Subscript[k3, i])+(Subscript[k4, i]))/(6.`)),{i,1,neq}];
τ=τ+h;datapoint={τ};Do[datapoint=Append[datapoint,(Subscript[w, i])],{i,1,neq}];data=Append[data,datapoint];
,{i,1,num}]

phasexz=Table[{data[[i,2]],data[[i,4]]},{i,1,Length[data]}];

rk4plot=ListPlot[phasexz,Joined->True,Frame->True,FrameLabel->{"(Subscript[y, 1])","(Subscript[y, 3])"},PlotStyle->{Blue},ImageSize->{500,Automatic},BaseStyle->20]

Notice that it took a little longer to solve, since we are solving a lot more equations. Here is one of the elements of the Jacobi matrix:

Jacobi11=Table[data[[i,5]],{i,1,num,4}];
ListPlot[Jacobi11,BaseStyle->20,AxesLabel->{"n","λ"}]

Notice that it is growing exponentially, which is an indication of chaos. The Lyapunov exponents are defined as (you get a Lyapunov exponent for each element of the Jacobi matrix, so you get 9 Lyapunov exponents in this case):

|(Subscript[x, 4])|=|((∂x(t))/(∂x(0)))|~(e^(λ t))
ln|(Subscript[x, 4])|~λ t
λ~(1/t)ln|(Subscript[x, 4])|
Typically we would use
λ=(lim_(t-> ∞))((1/t)ln|(Subscript[x, 4])|)

lyopunov=Table[(1/(data[[i,1]]))Log[Abs[data[[i,5]]]],{i,2,num}];
ListPlot[lyopunov,BaseStyle->20,AxesLabel->{"n","λ"}]

Here are all the Lyapunov exponents. Typically, we are interested in the largest Lyapunov exponent since it will dominant any computation and indicate the largest degree of chaos in the system.

alllyopunov=Table[((Log[Abs[data[[i,j]]]])/(data[[i,1]])),{j,5,neq+1},{i,2,num,10}];
ListPlot[alllyopunov]
ListPlot[alllyopunov,Joined->True,AxesLabel->{"n","λ"},BaseStyle->20]

So this system has a Lypunov exponent around 8, which is not a tremendous amount of chaos.

Lyapunov Exponents: Further Example

Here is a different example (different initial conditions) where we compute the Lyapunov exponent.

Clear[t,x,y,z,τ,f]
(Subscript[f, 1])=10 (-(Subscript[x, 1])+(Subscript[x, 2]));
(Subscript[f, 2])=28 (Subscript[x, 1])-(Subscript[x, 2])-(Subscript[x, 1]) (Subscript[x, 3]);
(Subscript[f, 3])=-((8 (Subscript[x, 3]))/3)+(Subscript[x, 1]) (Subscript[x, 2]);
(Subscript[f, 4])=10 (-(Subscript[x, 4])+(Subscript[x, 5]));
(Subscript[f, 5])=10 (-(Subscript[x, 7])+(Subscript[x, 8]));
(Subscript[f, 6])=10 (-(Subscript[x, 10])+(Subscript[x, 11]));
(Subscript[f, 7])=28 (Subscript[x, 4])-(Subscript[x, 5])-(Subscript[x, 4]) (Subscript[x, 3])-(Subscript[x, 1]) (Subscript[x, 6]);
(Subscript[f, 8])=28 (Subscript[x, 7])-(Subscript[x, 8])-(Subscript[x, 7]) (Subscript[x, 3])-(Subscript[x, 1]) (Subscript[x, 9]);
(Subscript[f, 9])=28 (Subscript[x, 10])-(Subscript[x, 11])-(Subscript[x, 10]) (Subscript[x, 3])-(Subscript[x, 1]) (Subscript[x, 12]);
(Subscript[f, 10])=-((8 (Subscript[x, 6]))/3)+(Subscript[x, 4]) (Subscript[x, 1])+(Subscript[x, 2]) (Subscript[x, 5]);
(Subscript[f, 11])=-((8 (Subscript[x, 9]))/3)+(Subscript[x, 7]) (Subscript[x, 1])+(Subscript[x, 2]) (Subscript[x, 8]);
(Subscript[f, 12])=-((8 (Subscript[x, 12]))/3)+(Subscript[x, 10]) (Subscript[x, 1])+(Subscript[x, 2]) (Subscript[x, 11]);
h=0.01;
num=2000;
neq=12;
τ=0.;
(Subscript[w, 1])=-3.;
(Subscript[w, 2])=2.;
(Subscript[w, 3])=5.;
(Subscript[w, 4])=1.;
(Subscript[w, 5])=0.;
(Subscript[w, 6])=0.;
(Subscript[w, 7])=0.;
(Subscript[w, 8])=1.;
(Subscript[w, 9])=0.;
(Subscript[w, 10])=0.;
(Subscript[w, 11])=0.;
(Subscript[w, 12])=1.;
data={};
datapoint={τ};
Do[datapoint=Append[datapoint,(Subscript[w, i])],{i,1,neq}];
data=Append[data,datapoint];
Do[subs={t->τ};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])],{i,1,neq}];Do[(Subscript[k1, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+(h/2)};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+((Subscript[k1, i])/2)],{i,1,neq}];Do[(Subscript[k2, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+(h/2)};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+((Subscript[k2, i])/2)],{i,1,neq}];Do[(Subscript[k3, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+h};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+(Subscript[k3, i])],{i,1,neq}];Do[(Subscript[k4, i])=h (Subscript[f, i])/.subs,{i,1,neq}];Do[(Subscript[w, i])=(Subscript[w, i])+(((Subscript[k1, i])+2 (Subscript[k2, i])+2 (Subscript[k3, i])+(Subscript[k4, i]))/(6.`)),{i,1,neq}];
τ=τ+h;
datapoint={τ};
Do[datapoint=Append[datapoint,(Subscript[w, i])],{i,1,neq}];data=Append[data,datapoint];
,{i,1,num}]

phasexz=Table[{data[[i,2]],data[[i,4]]},{i,1,Length[data]}];
rk4plot=ListPlot[phasexz,Joined->True,Frame->True,FrameLabel->{"(Subscript[w, 1])","(Subscript[w, 3])"},PlotStyle->{Blue},ImageSize->{500,Automatic},BaseStyle->20]

alllyopunov=Table[((Log[Abs[data[[i,j]]]])/(data[[i,1]])),{j,(1/2) (-1+(Sqrt[1+4 neq])),neq+1},{i,2,num,10}];

ListPlot[alllyopunov,AxesLabel->{"n","λ"},BaseStyle->20]

Here the degree of chaos associated with the initial point is a little larger than the previous example.

Lyapunov Exponents: 2D Anharmonic Oscillator (Hamiltonian Mechanics)

For the final example of computing Lyapunov exponents, we will look at the two dimensional Anharmonic oscillator (if b=0 this would be the Harmonic oscillator).

Clear[b]

Hamiltonian=(1/2)((Subscript[x, 2])^2)+(1/2)((Subscript[x, 1])^2)+(b/4) ((Subscript[x, 1])^4)
(Subscript[f, 1])=D[Hamiltonian,(Subscript[x, 2])]
(Subscript[f, 2])=-D[Hamiltonian,(Subscript[x, 1])]

If we want to calculate the Jacobi Matrix, we need to add derivatives with respect to the initial conditions:

Clear[t,x,y,z,τ,f]
b=10.;
Hamiltonian=((Subscript[x, 2]^2)/2)+((Subscript[x, 1]^2)/2)+((b (Subscript[x, 1]^4))/4);
(Subscript[f, 1])=(∂_(Subscript[x, 2]))Hamiltonian;
(Subscript[f, 2])=-(∂_(Subscript[x, 1]))Hamiltonian;
(Subscript[f, 3])=(Subscript[x, 5]);
(Subscript[f, 4])=(Subscript[x, 6]);
(Subscript[f, 5])=-(Subscript[x, 3])-3 b (Subscript[x, 1]^2) (Subscript[x, 3]);
(Subscript[f, 6])=-(Subscript[x, 4])-3 b (Subscript[x, 1]^2) (Subscript[x, 4]);
h=0.01;
num=2000;
neq=6;
τ=0.;
(Subscript[w, 1])=1.;
(Subscript[w, 2])=(Sqrt[2.]);
(Subscript[w, 3])=1.;
(Subscript[w, 4])=0.;
(Subscript[w, 5])=0.;
(Subscript[w, 6])=1.;
data={};
datapoint={τ};
Do[datapoint=Append[datapoint,(Subscript[w, i])],{i,1,neq}];
data=Append[data,datapoint];
Do[subs={t->τ};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])],{i,1,neq}];Do[(Subscript[k1, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+(h/2)};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+((Subscript[k1, i])/2)],{i,1,neq}];Do[(Subscript[k2, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+(h/2)};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+((Subscript[k2, i])/2)],{i,1,neq}];Do[(Subscript[k3, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+h};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+(Subscript[k3, i])],{i,1,neq}];Do[(Subscript[k4, i])=h (Subscript[f, i])/.subs,{i,1,neq}];Do[(Subscript[w, i])=(Subscript[w, i])+(((Subscript[k1, i])+2 (Subscript[k2, i])+2 (Subscript[k3, i])+(Subscript[k4, i]))/(6.`)),{i,1,neq}];
τ=τ+h;datapoint={τ};Do[datapoint=Append[datapoint,(Subscript[w, i])],{i,1,neq}];data=Append[data,datapoint];
,{i,1,num}]
phasexy=Table[{data[[i,2]],data[[i,3]]},{i,1,Length[data]}];
rk4plot=ListPlot[phasexy,Joined->True,Frame->True,FrameLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotStyle->{Blue},ImageSize->{500,Automatic},BaseStyle->20]

alllyopunov=Table[((Log[Abs[data[[i,j]]]])/(data[[i,1]])),{j,(1/2) (-1+(Sqrt[1+4 neq])),neq+1},{i,2,num,10}];

ListPlot[alllyopunov,AxesLabel->{"n","λ"},BaseStyle->20]

Compare to the Harmonic Oscillator (b=0):

Clear[t,x,y,z,τ,f]
b=0.0;
Hamiltonian=((Subscript[x, 2]^2)/2)+((Subscript[x, 1]^2)/2)+((b (Subscript[x, 1]^4))/4);
(Subscript[f, 1])=(∂_(Subscript[x, 2]))Hamiltonian;
(Subscript[f, 2])=-(∂_(Subscript[x, 1]))Hamiltonian;
(Subscript[f, 3])=(Subscript[x, 5]);
(Subscript[f, 4])=(Subscript[x, 6]);
(Subscript[f, 5])=-(Subscript[x, 3])-3 b (Subscript[x, 1]^2) (Subscript[x, 3]);
(Subscript[f, 6])=-(Subscript[x, 4])-3 b (Subscript[x, 1]^2) (Subscript[x, 4]);
h=0.01;
num=2000;
neq=6;
τ=0.;
(Subscript[w, 1])=1.;
(Subscript[w, 2])=(Sqrt[2.]);
(Subscript[w, 3])=1.;
(Subscript[w, 4])=0.;
(Subscript[w, 5])=0.;
(Subscript[w, 6])=1.;
data={};
datapoint={τ};
Do[datapoint=Append[datapoint,(Subscript[w, i])],{i,1,neq}];
data=Append[data,datapoint];
Do[subs={t->τ};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])],{i,1,neq}];Do[(Subscript[k1, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+(h/2)};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+((Subscript[k1, i])/2)],{i,1,neq}];Do[(Subscript[k2, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+(h/2)};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+((Subscript[k2, i])/2)],{i,1,neq}];Do[(Subscript[k3, i])=h (Subscript[f, i])/.subs,{i,1,neq}];
subs={t->τ+h};
Do[subs=Append[subs,(Subscript[x, i])->(Subscript[w, i])+(Subscript[k3, i])],{i,1,neq}];Do[(Subscript[k4, i])=h (Subscript[f, i])/.subs,{i,1,neq}];Do[(Subscript[w, i])=(Subscript[w, i])+(((Subscript[k1, i])+2 (Subscript[k2, i])+2 (Subscript[k3, i])+(Subscript[k4, i]))/(6.`)),{i,1,neq}];
τ=τ+h;datapoint={τ};Do[datapoint=Append[datapoint,(Subscript[w, i])],{i,1,neq}];data=Append[data,datapoint];
,{i,1,num}]
phasexy=Table[{data[[i,2]],data[[i,3]]},{i,1,Length[data]}];
rk4plot=ListPlot[phasexy,Joined->True,Frame->True,FrameLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotStyle->{Blue},ImageSize->{500,Automatic},BaseStyle->20]

alllyopunov=Table[((Log[Abs[data[[i,j]]]])/(data[[i,1]])),{j,(1/2) (-1+(Sqrt[1+4 neq])),neq+1},{i,2,num,10}];

ListPlot[alllyopunov,AxesLabel->{"n","λ"},BaseStyle->20]

References

You can learn more about the importance of predator prey models, limit cycles, and the Lorenz equation in Chapter 9 of Elementary Differential Equations by Boyce and Diprima.