Project Help Notebook by Keith Stroyan

Help NoteBook Copyright 1998 by Academic Press, Inc. - All rights reserved.
Modify this Notebook to prepare your project report.

Project 1001.4: The Stable Limit Cycle of the Van der Pol Equation

by A. Student and B. Pupil

with technical help from Da LabMonitor and Help Notebook

Special Mathematica Functions: flow2D[.] & vectorFlow2D[.]

Off[General::spell];
Off[General::spell1];
Off[InterpolatingFunction::dmval];
Off[ParametricPlot::pptr];
Off[Graphics::gprim];
Off[Graphics::gptn];
Off[NDSolve::ndsz];

Needs["Graphics`Colors`"]
	
colorQ[a_]:=(Head[a]===RGBColor);
grayQ[a_]:=(Head[a]===GrayLevel);
thickQ[a_]:=(Head[a]===Thickness);
absthickQ[a_]:=(Head[a]===AbsoluteThickness);
notruleQ[u_] := !(Head[u]===Rule);
ruleQ[u_] := (Head[u]===Rule);

Norm::usage = "Norm[expr] returns the Euclidean length of the
		vector, 'expr'.";
Norm[a_?VectorQ] := Sqrt[a.a];

tol = .5*10^-6;

Vector2D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics[Point[{0,0}]]];
	,						
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = (.075 Norm[ntip-ntail]/Norm[ntip-ntail]) (ntip-ntail) ;
	up = 0.15 {-Punit[[2]],Punit[[1]]};
	arrow = Line[{ntail,ntip,ntip-Punit+up,ntip,ntip-Punit-up}];
	]];
Return[Graphics[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

		
uArrow[vec_,opts___]:= 
	Module[{a,b},
		If[Norm[N[vec]]<1,b = vec,b=vec/Norm[vec]];
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0}};
				Return[Vector2D[a,a+b,opts]]
				];

flow2D[{f_,g_},{x_,y_,t_},initVals_,tFinal_,dt_,opns___] :=
 Module[{soln,arr,pt,crv,tt},
soln = Flatten[Map[
	NDSolve[{x[0]==#[[1]],y[0]==#[[2]],
	x'[t] == f[x[t],y[t]],y'[t] == g[x[t],y[t]]},
	{x[t],y[t]},{t,0,tFinal}]&,
										initVals],1];
											
pt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[4],Point[{#[[1]],#[[2]]}]}]&,
{x[t],y[t]}/.soln/.t->tt];
										
crv[tt_] := ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,tt},opns,
			AspectRatio->1,DisplayFunction->Identity];
			
Show[{crv[dt/100],pt[0]},opns,AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],pt[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
	Return[soln];	
];

vectorFlow2D[{f_,g_},{x_,y_,t_},initVals_,tFinal_,dt_,opns___] :=
 Module[{soln,arr,pt,crv,tt},
soln = Flatten[Map[
	NDSolve[{x[0]==#[[1]],y[0]==#[[2]],
	x'[t] == f[x[t],y[t]],y'[t] == g[x[t],y[t]]},
	{x[t],y[t]},{t,0,tFinal}]&,
										initVals],1];
								
	arr[tt_] := Map[uArrow[{f[#[[1]],#[[2]]],g[#[[1]],#[[2]]]},Tail->#,Red]&,
	{x[t],y[t]}/.soln/.t->tt];
	
pt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[4],Point[{#[[1]],#[[2]]}]}]&,
{x[t],y[t]}/.soln/.t->tt];
										
crv[tt_] := ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,tt},opns,
			AspectRatio->1,DisplayFunction->Identity];
			
Show[{crv[dt/100],pt[0],arr[dt/100]},opns,AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],arr[tt],pt[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
	Return[soln];	
];

My Project Introduction

Write this last.  Begin with the basic question and state your main finding clearly, but with few technical details.  You might even modify the title of this section to make it more interesting.

Limit Cycles

Chapters 20-24 of the main text study equilibrium points of dynamical systems.  An attracting equilibrium point is a "stable" point of the dynamics.  More complicated dynamics are discussed in the Mathematical Background section on stable limit cycles.  Following is a dynamical system with one of these attracting cycles.  This system tends to a single oscillation.

Van der Pol's Equations, f[x,y] & g[x,y] in terms of the parameter a:

Clear[x,y,a,f,g];
a = 3;
f[x_,y_] := y - a x(x^2/3 - 1);
g[x_,y_] := -x;
f[x,y]
g[x,y]

The initial x and y values

initXYs = Table[{y,-y},{y,-5,5,2}];

Solution using  f & g, initial values = initxys,  up to time = tfinal.

tFinal = 8;
dt = .1 ;	(*Time step*)
solutions=
  flow2D[{f,g},{x,y,t},initXYs,tFinal,dt,AxesLabel->{"x","y"},PlotRange->{{-6,6},{-6,6}}];

Highlite the cell with the graphs and press the <Command> & <y> keys together to see an animation of the flow associated with the differential equations

				x' = f[x,y]		dx = f[x,y] dt
				y' = g[x,y]       or        dy = g[x,y] dt

Microscopic Analysis at the Equilibrium Point

Clear[a,x,y];
aX = D[f[x,y],x];
aY = D[f[x,y],y];
bX = D[g[x,y],x];
bY = 0;
x = 0;
y = 0;
charPoly = (aX - r)(bY - r) - aY bX;
Expand[charPoly]

Solve[ charPoly == 0, r]