An Improved Method for Lyapunov Exponents Computation 
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/.

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

Compute estimate of maximum Lyapunov exponent using first method

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

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

J=JacobianMatrix[F[{(Subscript[y, 1])[t],(Subscript[y, 2])[t],(Subscript[y, 3])[t]}],{(Subscript[y, 1])[t],(Subscript[y, 2])[t],(Subscript[y, 3])[t]}]

J//MatrixForm

Y=Table[{(Subscript[y, i])[t],(Subscript[y, i+3])[t],(Subscript[y, i+6])[t]},{i,4,6}]

%//MatrixForm

J.Y//MatrixForm

E1=Flatten[Transpose[J.Y]]

EQ3=Table[D[(Subscript[y, i])[t],{t,1}]== F[{(Subscript[y, 1])[t],(Subscript[y, 2])[t],(Subscript[y, 3])[t]}][[i]],{i,1,3}]

EQ9=Table[D[(Subscript[y, i])[t],{t,1}]==E1[[i-3]],{i,4,12}]

YI9=Table[{(Subscript[y, i])[0]== If[Mod[i,4]== 0,1,0],(Subscript[y, i+3])[0]== If[Mod[i+3,4]== 0,1,0],(Subscript[y, i+6])[0]==If[Mod[i+6,4]== 0,1,0]},{i,4,6}]

YI3={(Subscript[y, 1])[0]== 19,(Subscript[y, 2])[0]== 20,(Subscript[y, 3])[0]== 50}

sol=NDSolve[Join[EQ9,EQ3,YI3,YI9],Table[(Subscript[y, i])[t],{i,1,12}],{t,0,40},MaxSteps-> Infinity,AccuracyGoal->2]

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

T=40;

PhiT=Transpose[First[Table[{(Subscript[y, i])[t],(Subscript[y, i+3])[t],(Subscript[y, i+6])[t]},{i,4,6}]/.sol/.t-> 40]]

Log[Norm[PhiT.u]]/T

Compute all three Lyapunov exponents using second method

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

YI9=Table[{(Subscript[y, i])[0]== If[Mod[i,4]== 0,1,0],(Subscript[y, i+3])[0]== If[Mod[i+3,4]== 0,1,0],(Subscript[y, i+6])[0]==If[Mod[i+6,4]== 0,1,0]},{i,4,6}]

YI3={(Subscript[y, 1])[0]== 19,(Subscript[y, 2])[0]== 20,(Subscript[y, 3])[0]== 50}

s={};

K=500;

Do[
      sol=NDSolve[Join[EQ9,EQ3,YI3,YI9],Table[(Subscript[y, i])[t],{i,1,12}],{t,0,0.1},MaxSteps-> Infinity,AccuracyGoal->2];
PhiT=Transpose[Table[{First[(Subscript[y, i])[t]/.sol/.t-> 0.1],First[(Subscript[y, i+3])[t]/.sol/.t-> 0.1],First[(Subscript[y, i+6])[t]/.sol/.t-> 0.1]},{i,4,6}]];
      W = GramSchmidt[PhiT, Normalized -> False];
      norms = Map[Norm, W];
      s = Append[s, norms];
      PhiT = W/norms;YI9=Table[{(Subscript[y, i+3])[0]== PhiT[[1,i]],(Subscript[y, i+6])[0]== PhiT[[2,i]],(Subscript[y, i+9])[0]==PhiT[[3,i]]},{i,1,3}];YI3={(Subscript[y, 1])[0]== (Subscript[y, 1])[t]/.sol/.t-> 0.1,(Subscript[y, 2])[0]==  (Subscript[y, 2])[t]/.sol/.t-> 0.1,(Subscript[y, 3])[0]== (Subscript[y, 3])[t]/.sol/.t-> 0.1 },
      {K}
    ];

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

Last[lce]

ListPlot[Transpose[lce]]