Re: Linearization of a system of differential equations
- To: mathgroup at christensen.cybernetics.net
- Subject: [mg829] Re: Linearization of a system of differential equations
- From: rknapp (Robert Knapp)
- Date: Mon, 24 Apr 1995 03:06:39 -0400
- Organization: Wolfram Research, Inc.
In article <3n1sja$avv at news0.cybernetics.net> alves at unixg.ubc.ca (Ivan Boris Alves) writes:
>
> I am in the process of exploring some systems of differential equations.
> The questions I have a fairly simple, such as determining whether a
> unique equilibrium is saddle point stable or not.
>
> Does anybody know of packages that simplify this process?
>
> ANy help would be appreciated
>
A package I wrote last year for undergraduate differential equations
projects may help. Here are some examples--the (short) package is at
the end of this post.
The project conatining the examples below and many others are
available via anonymous ftp from polar.bowdoin.edu
The examples are for the Lorenz equations.
Find the critical points for the Lorenz equations (with r a parameter)
In[1]:= Clear[x,y,z];v={x,y,z};
s=10;b=8/3;
f={s*(y-x), -x*z+r*x-y, x*y-b*z}
8 z
Out[1]= {10 (-x + y), r x - y - x z, x y - ---}
3
In[2]:= {c0,cm,cp}=ReplaceAll[v,#]& /@ Solve[f == {0,0,0},v]
Out[2]=
Sqrt[-8 + 8 r] Sqrt[-8 + 8 r]
{{0, 0, 0}, {-(--------------), -(--------------),
Sqrt[3] Sqrt[3]
Sqrt[-8 + 8 r] Sqrt[-8 + 8 r]
-1 + r}, {--------------, --------------, -1 + r}}
Sqrt[3] Sqrt[3]
LOAD the Package
In[3]:= <<Bifurcation.m;
Linearize around the critical point at {0,0,0} (the output is a matrix)
In[4]:= a0=Linearize[f,v,c0]
The linearized system is:
x' == -10 x + 10 y
y' == r x - y
-8 z
z' == ----
3
8
Out[4]= {{-10, 10, 0}, {r, -1, 0}, {0, 0, -(-)}}
3
This plots the eigenvalues of the matrix as a function of the
parameter r (for c.c. pairs, the imag parts are plotted greeen adn the
real part is plotted blue)
In[5]:= EigenPlot[a0,{r,0,100}];
Linearize about the critical point "ap"
In[30]:= ap=Linearize[f,v,cp]
The linearized system is:
x' == -10 x + 10 y
Sqrt[-8 + 8 r] z
y' == x - y - ----------------
Sqrt[3]
Sqrt[-8 + 8 r] x Sqrt[-8 + 8 r] y 8 z
z' == ---------------- + ---------------- - ---
Sqrt[3] Sqrt[3] 3
Out[30]=
Sqrt[-8 + 8 r]
{{-10, 10, 0}, {1, -1, -(--------------)},
Sqrt[3]
Sqrt[-8 + 8 r] Sqrt[-8 + 8 r] 8
{--------------, --------------, -(-)}}
Sqrt[3] Sqrt[3] 3
Again, you can plot the eigenvalues with EigenPlot. Interesting
regions to look at are for r between 1.32 and 1.36 (a bifurcation from
two real eigenvalues to a c.c. pair occurs) and r between 24.732 and
24.753 (the real part of the c.c. pair becomes positive, making this
critical point unstable-which eventually sets the stage for chaos at a
slightly higher value of r.
Here is the package:
BeginPackage["Bifurcation`",{"Utilities`FilterOptions`"}];
EigenPlot::usage="EigenPlot[a,{r,r1,r2,nr},opts] produces a plot of the
eigenvalues of the matrix a as a function of the bifurcation parameter r by
sampling at nr points between r1 and r1. Real eigenvalues are plotted in red. For complex eigenvalues, the real part is plotted in blue and the imaginary part is plotted in green. Valid options are any Graphics option and Tolerance: if Tolerance->tol, then if Abs[Im[lambda]]< tol, the eigenvalue is considered real. The default is Tolerance->10^(-12)."
Linearize::usage="Linearize[f,v,vc] produces a linearization of the system of
DE's v[[1]]'=f[[1]],etc. at the critical point vc. The output is the matrix for
the linear system."
Options[EigenPlot]={
Tolerance->10^(-12)}
Begin["`Private`"];
Linearize[fin_,vin_,vcin_]:=Module[{f=fin,v=vin,vc=vcin,i,matrix,sys},
If[Length[fin] == 0,f={fin}];
If[Length[vin] == 0,v={vin};vc->{vcin}];
matrix=Table[(D[f[[i]],#]& /@ v),{i,1,Length[f]}] /. Thread[v->vc];
sys=Simplify[f /. Thread[v->vc]]+matrix.v;
Print["The linearized system is:"];
Do[Print[v[[i]]'==sys[[i]]],{i,1,Length[v]}];
matrix]
EigenPlot[a_,{r_,r1_,r2_,nr_:100},opts___]:=
Module[{i,n=Length[a],tol,real={},cre={},cim={},points={}},
tol = Tolerance /. {opts} /. Options[EigenPlot];
Do[evals=Eigenvalues[N[a /. r->rr]];
Do[ev=evals[[i]];
If[Abs[Im[ev]] < tol,
AppendTo[real,{rr,Re[ev]}],
AppendTo[cre,{rr,Re[ev]}];AppendTo[cim,{rr,Im[ev]}]],
{i,1,n}],
{rr,r1,r2,(r2-r1)/nr}];
If[Length[real] > 0,AppendTo[points,{{RGBColor[1,0,0],(Point /@ real)}}]];
If[Length[cre] > 0,points=Join[points,{{RGBColor[0,0,1],(Point /@ cre)}},
{{RGBColor[0,1,0],(Point /@ cim)}}]];
Show[Graphics[{points}],
FilterOptions[Graphics,opts],
Axes->True,GridLines->Automatic]]
End[];
EndPackage[];