Lyapunov Exponents Computation 
for Various Non-Linear Dynamic Problems
Authors' Data: Housam Binous & Nasri Zakia
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowlegment : 

1- Sandri, M. "Numerical Calculation of Lyapunov Exponents." Mathematica J. 6, 78-84, 1996. http://library.wolfram.com/infocenter/Articles/2902/. 

2- Govorukhin V.N. E-Mail: vgov@math.rsu.ru, matds@nm.ru 
MATDS is a MATLAB-based program for dynamical system investigation. It works with version 6.5* of MATLAB. 
http://kvm.math.rsu.ru/matds/

One can compare the results obtained in the present notebook (using Mathematica' s methods by M. Sandir) with the Matlab files in the compressed folder lyapunov2.zip (using the method by Govorukhin which is inpired from earlier work by A.Wolf, J.B.Swift, H.L.Swinney, and J.A.Vastano, "Determining Lyapunov Exponents from a Time Series," Physica D, Vol .16, pp.285 - 317, 1985).

Test with Lorentz system. We get a chaotic behavior for this choice of parametres. Max  lyapunov exponent is positive = 1.41

F[{x_,y_,z_}]:={16 (y-x),x (45.92-z)-y,x y -4z}

x0={19,20,50};

RKStep[f_,y_,y0_,dt_]:=Module[{k1,k2,k3,k4},k1=dt N[f/.Thread[y->y0]];
k2=dt N[f/.Thread[y->y0+k1/2]];
k3=dt N[f/.Thread[y->y0+k2/2]];
k4=dt N[f/.Thread[y->y0+k3]];
y0+(k1+2 k2+2 k3+k4)/6]

IntVarEq[F_List, DPhi_List, x_List, Phi_List, x0_List, 
        Phi0_List, {t1_, dt_}] :=
Module[{n, f, y, y0, yt},
    n = Length[x0];
    f = Flatten[Join[F, DPhi]];
    y = Flatten[Join[x, Phi]];
    y0 = Flatten[Join[x0, Phi0]];
    yt = Nest[ RKStep[f, y, #, N[dt]]&, N[y0], Round[N[t1/dt]] ];
    {First[#], Rest[#]}& @ Partition[yt, n] ]

JacobianMatrix[funs_List, vars_List] := Outer[D, funs, vars]

T=20;

stepsize=0.02;

n=Length[x0];

x=Array[a,n];

Phi=Array[b,{n,n}];

DPhi = Phi.Transpose[JacobianMatrix[F[x], x]];

Phi0=IdentityMatrix[n];

{xT,PhiT}=IntVarEq[F[x],DPhi,x,Phi,x0,Phi0,{T,stepsize}]

u=Table[Random[],{n}]

Log[Norm[PhiT.u]]/T

T=0.5;

K=800;

s={};

PhiT=Phi0;

<<"C:\\Program Files\\Wolfram Research\\Mathematica\\6.0\\AddOns\\LegacyPackages\\LinearAlgebra\\Orthogonalization.m"

Do[
      {xT, PhiT} = IntVarEq[F[x], DPhi, x, Phi, xT, PhiT, {T, stepsize}];
      W = GramSchmidt[PhiT, Normalized -> False];
      norms = Map[Norm, W];
      s = Append[s, norms];
      PhiT = W/norms,
      {K}
    ];

lce=Rest[FoldList[Plus,0,Log[s]]]/(T Range[K]);

Last[lce]

ListPlot[Transpose[lce]]

Test with Rosler system. We get a chaotic behavior for this choice of parametres. Extimate for the Max  lyapunov exponent is positive = 0.089.

F[{x_,y_,z_}]:={-y-z,x + a y,b+ z (x-c)}

x0={.2,.3,.5};

RKStep[f_,y_,y0_,dt_]:=Module[{k1,k2,k3,k4},k1=dt N[f/.Thread[y->y0]];
k2=dt N[f/.Thread[y->y0+k1/2]];
k3=dt N[f/.Thread[y->y0+k2/2]];
k4=dt N[f/.Thread[y->y0+k3]];
y0+(k1+2 k2+2 k3+k4)/6]

IntVarEq[F_List, DPhi_List, x_List, Phi_List, x0_List, 
        Phi0_List, {t1_, dt_}] :=
Module[{n, f, y, y0, yt},
    n = Length[x0];
    f = Flatten[Join[F, DPhi]];
    y = Flatten[Join[x, Phi]];
    y0 = Flatten[Join[x0, Phi0]];
    yt = Nest[ RKStep[f, y, #, N[dt]]&, N[y0], Round[N[t1/dt]] ];
    {First[#], Rest[#]}& @ Partition[yt, n] ]

JacobianMatrix[funs_List, vars_List] := Outer[D, funs, vars]

T=500;

stepsize=0.02;

a=0.15;b=0.20;c=10.0;

n=Length[x0];

x=Array[a,n];

Phi=Array[b,{n,n}];

DPhi = Phi.Transpose[JacobianMatrix[F[x], x]];

Phi0=IdentityMatrix[n];

{xT,PhiT}=IntVarEq[F[x],DPhi,x,Phi,x0,Phi0,{T,stepsize}]

u=Table[Random[],{n}]

Log[Norm[PhiT.u]]/T

T=1;

K=1000;

s={};

PhiT=Phi0;

<<"C:\\Program Files\\Wolfram Research\\Mathematica\\6.0\\AddOns\\LegacyPackages\\LinearAlgebra\\Orthogonalization.m"

Do[
      {xT, PhiT} = IntVarEq[F[x], DPhi, x, Phi, xT, PhiT, {T, stepsize}];
      W = GramSchmidt[PhiT, Normalized -> False];
      norms = Map[Norm, W];
      s = Append[s, norms];
      PhiT = W/norms,
      {K}
    ];

lce=Rest[FoldList[Plus,0,Log[s]]]/(T Range[K]);

Last[lce]

ListPlot[Transpose[lce],Joined->True]

Test with Rosler system. We get a Period 2 behavior for this choice of parametres. Estimate for the Max  lyapunov exponent is almost equal to 0. Thus all exponents are less than zero.  Periodic behavior.

x0={.2,.3,.5};

RKStep[f_,y_,y0_,dt_]:=Module[{k1,k2,k3,k4},k1=dt N[f/.Thread[y->y0]];
k2=dt N[f/.Thread[y->y0+k1/2]];
k3=dt N[f/.Thread[y->y0+k2/2]];
k4=dt N[f/.Thread[y->y0+k3]];
y0+(k1+2 k2+2 k3+k4)/6]

IntVarEq[F_List, DPhi_List, x_List, Phi_List, x0_List, 
        Phi0_List, {t1_, dt_}] :=
Module[{n, f, y, y0, yt},
    n = Length[x0];
    f = Flatten[Join[F, DPhi]];
    y = Flatten[Join[x, Phi]];
    y0 = Flatten[Join[x0, Phi0]];
    yt = Nest[ RKStep[f, y, #, N[dt]]&, N[y0], Round[N[t1/dt]] ];
    {First[#], Rest[#]}& @ Partition[yt, n] ]

JacobianMatrix[funs_List, vars_List] := Outer[D, funs, vars]

T=500;

stepsize=0.1;

a=0.2;b=0.2;c=3.0;

F[{x1_,y1_,z1_}]:={-y1-z1,x1 + a y1,b+ z1 (x1-c)}

n=Length[x0];

x=Array[a,n];

Phi=Array[b,{n,n}];

DPhi = Phi.Transpose[JacobianMatrix[F[x], x]];

Phi0=IdentityMatrix[n];

{xT,PhiT}=IntVarEq[F[x],DPhi,x,Phi,x0,Phi0,{T,stepsize}]

u=Table[Random[],{n}]

Log[Norm[PhiT.u]]/T

T=1;

K=1000;

s={};

PhiT=Phi0;

<<"C:\\Program Files\\Wolfram Research\\Mathematica\\6.0\\AddOns\\LegacyPackages\\LinearAlgebra\\Orthogonalization.m"

Do[
      {xT, PhiT} = IntVarEq[F[x], DPhi, x, Phi, xT, PhiT, {T, stepsize}];
      W = GramSchmidt[PhiT, Normalized -> False];
      norms = Map[Norm, W];
      s = Append[s, norms];
      PhiT = W/norms,
      {K}
    ];

lce=Rest[FoldList[Plus,0,Log[s]]]/(T Range[K]);

Last[lce]

ListPlot[Transpose[lce],Joined->True]

Lyapunov Exponents for the Dynamic behavior of non - isothermal chemical system - Reference : S. K. Scott and A. S. Tomlin, Period doubling and other complex bifurcations in non-isothermal chemical systems, Phil. Trans. R. Soc. Lond. A (1990) 332, 51-68.

For μ = 0.498. Estimate for the Max  lyapunov exponent is almost  zero (-0.008).  Thus all exponents are less than zero. 
Confirmation of  Period 1 behavior.

κ=10^-3;γ=0.05;δ=0.025;

F[{x1_,y1_,z1_}]:={μ Exp[δ z1]-x1 y1^2-κ x1,
x1 y1^2+κ x1-y1,y1-γ z1}

μ=0.498;

F[{x1,y1,z1}]

x0={1,0.5,0};

DPhi = Phi.Transpose[JacobianMatrix[F[x], x]];

T=100;

stepsize=0.01;

{xT,PhiT}=IntVarEq[F[x],DPhi,x,Phi,x0,Phi0,{T,stepsize}]

u=Table[Random[],{n}]

Log[Norm[PhiT.u]]/T

T=0.1;

PhiT=Phi0;

K=2000;

s={};

Do[
      {xT, PhiT} = IntVarEq[F[x], DPhi, x, Phi, xT, PhiT, {T, stepsize}];
      W = GramSchmidt[PhiT, Normalized -> False];
      norms = Map[Norm, W];
      s = Append[s, norms];
      PhiT = W/norms,
      {K}
    ];

lce=Rest[FoldList[Plus,0,Log[s]]]/(T Range[K]);

Last[lce]

ListPlot[Transpose[lce],Joined->True]

For μ = 0.5586. Max  lyapunov exponent is positive (0.035). 
Confirmation of Chaotic behavior.

μ=0.5586;

F[{x1,y1,z1}]

x0={1,0.5,0};

DPhi = Phi.Transpose[JacobianMatrix[F[x], x]];

T=500;

stepsize=0.01;

{xT,PhiT}=IntVarEq[F[x],DPhi,x,Phi,x0,Phi0,{T,stepsize}]

u=Table[Random[],{n}]

Log[Norm[PhiT.u]]/T

T=0.1;

PhiT=Phi0;

K=2000;

s={};

Do[
      {xT, PhiT} = IntVarEq[F[x], DPhi, x, Phi, xT, PhiT, {T, stepsize}];
      W = GramSchmidt[PhiT, Normalized -> False];
      norms = Map[Norm, W];
      s = Append[s, norms];
      PhiT = W/norms,
      {K}
    ];

lce=Rest[FoldList[Plus,0,Log[s]]]/(T Range[K]);

Last[lce]

ListPlot[Transpose[lce],Joined->True]