Copyright (C) 2005 Chad E. Brassil; Last modified May 20, 2008.  Note, this has not been fully tested in Mathematica 6

Set-up

<<"Statistics`DataManipulation`" 
<<"PlotLegends`"

(* This is a special package that comes with the book "Mathematica Graphics" *)
(* It can be downloaded at http://library.wolfram.com/infocenter/Books/3753/ *)
<<ExtendGraphics`Ticks`

(* modify this for your computer *)
mypath="c:/documents and settings/Chad/Desktop/";

SaveEPS[name_,function_]:=Export[ToString[Row[{mypath,name,".eps"}]],function,"EPS",ConversionOptions->{"IncludeSpecialFonts"->True}];

tickChadY[xmin_,xmax_]:=TickFunction[xmin,xmax,MinorLength->{0,0},MajorLength->{0.007,0},MajorStyle->{AbsoluteThickness[1.5]},TickNumbers->{6,32}]
tickChad[xmin_,xmax_]:=TickFunction[xmin,xmax,MinorLength->{0,0},MajorLength->{0.007,0},MajorStyle->{AbsoluteThickness[1.5]},TickNumbers->{8,32}]

Clear[PlotStyleFunction]; SetOptions[Plot,Frame->True,FrameStyle->AbsoluteThickness[1.5`],ImageSize->72 3,BaseStyle->{FontFamily->"Swiss721BT-Roman",FontSize->8},PlotRegion->{{0,1},{0,0.995`}},FrameTicks->{tickChad,tickChad,{},{}},PlotRange->{0,Automatic},PlotPoints->50,AxesStyle->AbsoluteThickness[1.5`],PlotStyle->PlotStyleFunction]; PlotStyleFunction:=AbsoluteThickness[1.5`];

Clear[ListPlotStyleFunction]; SetOptions[ListPlot,Frame->True,FrameStyle->AbsoluteThickness[1.5`],ImageSize->72 3,BaseStyle->{FontFamily->"Swiss721BT-Roman",FontSize->8},PlotRegion->{{0,1},{0,0.995`}},FrameTicks->{tickChad,tickChad,{},{}},AxesStyle->AbsoluteThickness[1.5`],PlotStyle->ListPlotStyleFunction]; ListPlotStyleFunction:=AbsolutePointSize[2.5`];

FigureLabel[graphic_,label_String]:=ShowLegend[graphic,{GrayLevel[1]&,1.1`,"",label,LegendPosition->{-0.95`,0.5`},LegendShadow->{0,0},LegendTextSpace->100,LegendSize->{0.1`,0.1`},LegendBorder->GrayLevel[1],BaseStyle->{FontFamily->Evaluate[FontFamily/.Options[Plot,TextStyle][[1,2]]],FontSize->12}}];

TwoAxisPlot[{f_, g_}, {x_, min_, max_}, frange_, grange_, (opts___)?OptionQ] := Module[{old, new, scale, pts, gstuff, fstuff}, {fm, fM} = frange; {gm, gM} = grange; scale[var_] = ((var - gm)*(fM - fm))/(gM - gm) + fm ; old = AbsoluteOptions[ Plot[g, {x, min, max}, Frame -> True, PlotRange -> grange, DisplayFunction -> Identity], FrameTicks][[1, 2, 2]]; new = (Prepend[Rest[#1], scale[First[#1]]] & ) /@ old; Plot[{f, scale[g]}, {x, min, max}, Frame -> True, PlotRange -> frange + {0, .05 (fM - fm)}, FrameTicks -> {Automatic, Automatic, None, new}, PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 1]}}, FrameStyle -> {{}, {RGBColor[1, 0, 0]}, {}, {RGBColor[0, 0, 1]}}, opts] ]

TwoAxisPlot[{x^3, Sqrt[x]}, {x, 0, 4}, {0, 64}, {0, 2}] ;

chopToInteger[x_?NumericQ]:=If[IntegerPart[x]-x==0,IntegerPart[x],x];
chopToInteger[x_]:=x;
SetAttributes[chopToInteger,Listable];

TwoAxisPlotChad[{f_, g_}, {x_, min_, max_}, frange_, grange_, (opts___)?OptionQ] := Module[{(*old, new, scale,*) pts, gstuff, fstuff}, {fm, fM} = frange; {gm, gM} = grange; scale[var_] = ((var - gm)*(fM - fm))/(gM - gm) + fm ; old = AbsoluteOptions[ Plot[g, {x, min, max}, Frame -> True, PlotRange -> grange, DisplayFunction -> Identity,FrameTicks -> {tickChad, tickChad}], FrameTicks][[1, 2, 2]]; new = (Prepend[Rest[#1], scale[First[#1]]] & ) /@ old; 
new=chopToInteger[new];
Plot[Evaluate[Flatten[{f, scale[g]}]], {x, min, max}, Frame -> True, PlotRange -> frange + {0, .05 (fM - fm)}, FrameTicks -> {tickChad, tickChad, None, new}, opts] ]

General Differential Equation Solver

This compact code can be used to analyze differential equations.  It can handle from 1 to many variables, just keep them in the same order for deqs, vars, and start.

Equation Setup

deqs = a list of the right-hand side of the differential equations
vars = the variables, in the same order as the deqs
start = initial conditions
tmax = maximum time to plot time series

deqs={
(Subscript[δ, b])-(Subscript[μ, b])sb-(Subscript[β, bb])sb ib-(Subscript[β, db])sb id,
(Subscript[β, bb]) sb ib+(Subscript[β, db])sb id-(Subscript[α, b])ib,
(Subscript[δ, d])-(Subscript[μ, d])sd-(Subscript[β, dd])sd id-(Subscript[β, bd])sd ib,
(Subscript[β, dd])sd id+(Subscript[β, bd])sd ib-(Subscript[α, d])id
};
vars={sb,ib,sd,id};
start={0.1,0.2,0.3,0.4};
paramset={(Subscript[δ, b])->2,(Subscript[δ, d])->2,(Subscript[μ, b])->1.8,(Subscript[μ, d])->2,(Subscript[β, db])->2.2,(Subscript[β, bb])->1.8,(Subscript[β, dd])->1.8,(Subscript[β, bd])->1.6};
tmax=400;

(* run this just once to remove zero starting values from vars and deqs *)
initialZeroPositions=Position[start,0.0];
paramset=Join[paramset,Thread[Extract[vars,initialZeroPositions]->0]];
start=Delete[start,initialZeroPositions];
vars=Delete[vars,initialZeroPositions];
deqs=Delete[deqs,initialZeroPositions];

Use the following code to specify dashing, thickness, and colors for the vars

PlotStyleFunction:={{AbsoluteDashing[{1.0}],AbsoluteThickness[1.5],RGBColor[1,0,0]},{AbsoluteDashing[{1.0}],AbsoluteThickness[1.5],RGBColor[0,0,1]},{AbsoluteDashing[{0.0}],AbsoluteThickness[1.5],RGBColor[1,0,0]},{AbsoluteDashing[{0.0}],AbsoluteThickness[1.5],RGBColor[0,0,1]}};

This code can be used to create a legend which can be called throughout the notebook

lineWidth=0.6;
pLegend=Module[{i},Show[Graphics[Table[Text[vars[[i]],{i,0}],{i,Length[vars]}]],Table[Graphics[Append[PlotStyleFunction[[i]],Line[{{i-(lineWidth/2),-.4},{i+(lineWidth/2),-.4}}]]],{i,Length[vars]}],Options[Plot,BaseStyle][[1]],PlotRange->All,ImageSize->{Options[Plot,ImageSize][[1,2]],20},AspectRatio->Automatic]]

Equilibrium

equilib=Solve[Thread[deqs==0],vars]

(* Lists which equilibrium have positive densities for all variables *)
positivePositions=Module[{i},Table[Nor@@TrueQ/@(Thread[vars==0]/.equilib[[i]]),{i,Length[equilib]}]]

(* first all-positive equilibrium *)
firstPositivePosition=Position[positivePositions,True][[1,1]];
posEquilib=equilib[[firstPositivePosition]]

(* numerical *)
Solve[Thread[deqs==0],vars]/.paramset

Attempt Null Clines

This section uses a special package that can be downloaded from http://www.math.armstrong.edu/faculty/hollis/mmade/

Will start initial values and run them for tmax or until they go off the screen

<<DiffEqs`DEGraphics`

maxVar1=2;
maxVar2=2.5`;
PhasePlot[deqs/.paramset,{g,0,tmax},{vars[[1]],0,maxVar1},{vars[[2]],0,maxVar2},StepSize->0.01`,ScaleArrows->0.5`,Frame->True,FrameLabel->{vars[[1]],vars[[2]]},BaseStyle->{FontFamily->"Swiss721BT-Roman"},ImageSize->350,ShowNullclines->True,ShowEquilibriumPoints->True];

PhasePlot[deqs/.paramset,{g,0,tmax},{vars[[1]],0,maxVar1},{vars[[2]],0,maxVar2},StepSize->0.01`,InitialPoints->{start},ScaleArrows->0.5`,Frame->True,FrameLabel->{vars[[1]],vars[[2]]},BaseStyle->{FontFamily->"Swiss721BT-Roman"},ImageSize->350,ShowNullclines->True,ShowEquilibriumPoints->True];

Plot of Density vs. Time

view=100;

(* find time-solution *)
addT=Module[{i},Table[vars[[i]]->Map[vars[[i]],{t}][[1]],{i,Length[vars]}]];
outset=NDSolve[Join[Thread[D[vars/.addT,t]==(deqs/.addT)],Thread[vars==start]/.addT/.t->0]/.paramset,vars,{t,0,tmax},MaxSteps->50*tmax,AccuracyGoal->Infinity];

Show[pLegend]
pA=Plot[Evaluate[vars/.addT/.outset],{t,tmax-view,tmax},PlotRange->All,FrameLabel->{t,""}]

fig=FigureLabel[pA,"A"];

(*SaveEPS["figA",fig];*)

Attempt Analytical Stability

numberOfTerms[x_]:=If[TrueQ[Head[x]==Plus],Length[x],1];

i=firstPositivePosition; (* choose which equilibrium to examine *)
jacSub=Outer[D,deqs,vars]/.equilib[[i]];
(* Note, this syntax will also produce the Jacobian *)
(* D[deqs,{vars}] *)
Simplify[MatrixForm[jacSub]]
jacNoSub=Outer[D,deqs,vars];
Simplify[MatrixForm[jacNoSub]]

Above is the Jacobian with substitutions and without substitutions.  Below is attempt at the simpliest combination of the above.

jac=Table[
If[numberOfTerms[Expand[jacSub[[ijac,jjac]]]]<numberOfTerms[Expand[jacNoSub[[ijac,jjac]]]],Simplify[jacSub[[ijac,jjac]]],Simplify[jacNoSub[[ijac,jjac]]]]
,{ijac,Length[vars]},{jjac,Length[vars]}];
jac=jac/.Thread[(vars/.posEquilib)->vars];
MatrixForm[jac]

1) Brute force
-stable if the real part of the largest eigenvalue is negative
-may be able to simplify be hand knowing equilibrium values need to be positive

Eigenvalues[jac]

2)  Special rule for 2 x 2 matrix: stable if trace < 0 and det > 0
For stable or unstable, oscillatory if det>((trace^2)/4)

(* for a 2x2 matrix, the determinant > 0 for stability *)
Det[jac]//Simplify

(* for a 2x2 matrix, the trace < 0 for stability *)
Tr[jac]//Simplify

3) Routh-Hurrwitz criteria for a 3 x 3 matrix.
characteristic equation: (λ^3)+(Subscript[A, 1])(λ^2)+(Subscript[A, 2])λ+(Subscript[A, 3])=0
stable if:
	1) (Subscript[A, 1])>0
	2) (Subscript[A, 3])>0
	3) (Subscript[A, 1])(Subscript[A, 2])>(Subscript[A, 3])

(* characterisic polynomial *)
Collect[Det[jac-λ IdentityMatrix[Length[jac]]],λ]

4) Some matrices can be reduced to a 2x2 matrix.  
A matrix in the form (({{a, x}, {0, b}})) is called an upper triangular matrix and has the eigenvalues {a,b}.  I can also replace each entry in the matrix with another matrix, represented here by a capital letter, M = (({{A, X}, {0, B}})).  If the matrix A and B are each square matrices, then the eigenvalues for the matrix M will be the eigenvalues for both A and B.  Here is an example.  The matrix (({{a, b, c}, {d, e, f}, {0, 0, g}}))can put into the form of matrix M with A = (({{a, b}, {d, e}})), X = (({{c}, {f}})), B = (g), and 0 = (({{0, 0}})).  I can now examine A as a 2x2 matrix and g is the third eigenvalue.  This technique can be useful when looking at invasion analysis.

Numerical Stability

Numerically we can always examine the eigenvalues of the Jacobian matrix at the equilibrium point.  Sometimes these can involve complex numbers.  If the real part of the largest eigenvalue is negative, the equilibrium is stable.

jacSolved=jacSub/.paramset;
MatrixForm[jacSolved]

Eigenvalues[jacSolved]

(* stable if < 0 *)
Max[Re[Eigenvalues[jacSolved]]]

Lyapunov Exponent

The largest Lyapunov exponent can be used to characterize a system as chaotic.  If the exponent is > 0, then the system is chaotic.

Example Equation

Same format as above

deqs={
((α p n1)/(1+β n1+β n2))+((α p n2)/(1+β n1+β n2))-d p,
n1(1-n1)-((n1 p)/(1+β n1+β n2)),
n2 (1-n2)-((n2 p)/(1+β n1+β n2))
};
vars={p,n1,n2};
start={0.1,0.05,0.02};
paramset={α->4.5+4.05 Sin[0.490987 t],β->0.02,d->0.2};
tmax=10000;

transtime is the the time units to ignore so that the long-term dynamics are examined

transtime=3000;

Find largest Lyapunov using built-in integrator

Use the same time-solution code from "Plot of Density vs. Time" section

(* find time-solution *)
addT=Table[vars[[i]]->Map[vars[[i]],{t}][[1]],{i,Length[vars]}];
outset=NDSolve[Join[Thread[D[vars/.addT,t]==(deqs/.addT)],Thread[vars==start]/.addT/.t->0]/.paramset,vars,{t,0,tmax},MaxSteps->50*tmax,AccuracyGoal->Infinity];

The following code assumes you have already defined the equations and found the time-solution, as illustrated just above

Module[{n,Phi0,DPhi,PhiNoT},
n=Length[start];
PhiNoT=Transpose[Table[Unique["b"],{n},{n}]];
Phi=Table[PhiNoT[[i,j]][t],{i,n},{j,n}];
Phi0=IdentityMatrix[n];
DPhi=Outer[D,deqs/.addT,vars/.addT].Phi;
derivhead=Flatten[D[Phi,t]];
initialconditions=Thread[Flatten[Phi]==Flatten[Phi0]]/.t->transtime;
eqns=Union[Thread[derivhead==Flatten[DPhi]]/.outset,initialconditions];
solLyap=NDSolve[eqns/.paramset,Flatten[PhiNoT],{t,transtime,tmax},MaxSteps->100000][[1]];
PhiT[tsub_]:=Phi/.solLyap/.t->tsub;
Log[Norm[PhiT[tmax]]]/(tmax-transtime)
]

Plot[((Log[Norm[PhiT[t]]])/(t-transtime)),{t,transtime+1,tmax},PlotRange->{-0.01`,0.01`}]

Lyapunovs using Runga-Kuta and LCE package

Package is written by Marco Sandri.  Can be downloaded from http://www.msandri.it/soft.html

<<Lyapunov`LCE`

F[{x_,y_,z_}]:={σ(y-x),
x(r-z)-y,
x y-b z}/.{σ->16.0,r->45.92,b->4.0};
initials={19,20,50};

LCEsC[F,initials,0.1,800,20,0.02,LCEsPlot->True]

Find Period or Max/Min

setup

(* This version may have problems with solutions in which the starting point happens to sit at the bottom of the cycle and the difference never changes sign *)
(* I'm not using a distance measure here but rather a summed difference so that the EventLocator can monitor when the function switches sign *)
NFindPeriod[deqs_,vars_,paramset_,tmax_,outset_,distanceThreshold_]:=
Block[{startAgain,timeFirstGuess,addT,perOutSet,},
addT=Table[vars[[i]]->Map[vars[[i]],{t}][[1]],{i,Length[vars]}];
startAgain=vars/.addT/.outset[[1]]/.t->tmax;
perOutSet=NDSolve[Join[Thread[D[vars/.addT,t]==(deqs/.addT)],Thread[vars==startAgain]/.addT/.t->tmax]/.paramset,vars,{t,tmax,2*tmax},MaxSteps->∞,Method->{EventLocator, "Event"->Plus@@(startAgain-vars/.addT),"EventAction":>If[Norm[startAgain-vars/.addT]<distanceThreshold,Throw[Null,"StopIntegration"]], "Direction"->All}];timeFirstGuess=perOutSet[[1,1,2,1,1,2]]-tmax
];
NFindPeriod::usage="NFindPeriod[ deqs, vars, paramset, tmax, outset, distanceThreshold ] finds the period for a set of differential equations deqs with variable vars and parameters paramset.  The interpolating function solution evaluated up to tmax needs be included as outset.  If no periodic solution is found, the function will return tmax.";

use

NFindPeriod[deqs,vars,paramset,tmax,outset, 0.01]

Another set of useful code finds the max and min of a the solution of NDSolve

(* find max and min *)
(* currently set to examine the second variabile *)
rangeNVarsForced=Block[{endPoint,d0List,minmaxOutset},
endPoint=vars/.addT/.outset[[1]]/.t->tmax;
d0List={};
minmaxOutset=NDSolve[Join[Thread[D[vars/.addT,t]==(deqs/.addT)],Thread[vars==endPoint]/.addT/.t->tmax]/.paramset,vars,{t,tmax,tmax+per},MaxSteps->∞,Method->{EventLocator, "Event"->deqs[[2]]/.addT/.paramset,"EventAction":>AppendTo[d0List,t], "Direction"->All}];
d0List=vars[[2]]/@d0List/.minmaxOutset[[1]];
{Max[d0List],Min[d0List]}
];