Part IV. Linear Differential Equations
Section 3. Two Dimensional Systems
A second order differential equation in normal form
y'' = f( t, y, y' )
can always be converted into an equivalent system of two first order equations
y' = v
v' = f( t, y, v )
If the function f does not depend upon t, then the second order equation (and the system) is called autonomous. More generally, a two dimensional system is a pair of first order differential equations of the form
x' = f ( x, y, t )
y' = g ( x, y, t )
The primes denote differentiation with respect to t. If neither f nor g depend upon t the system is autonomous. In this section Mathematica is used to analyze the solutions to such systems. By now you are familiar with the tools we need: DSolve, NDSolve, PlotVectorField, and ParametricPlot.
In[1]:=
<<Graphics`PlotField`
The solution to a two dimensional system is a pair of functions x(t) and y(t). Solution behavior is analyzed using time series and plots of phase plane trajectories of the form ( x(t), y(t) ) when the system is autonomous. State space trajectories ( x(t), y(t), t ) can be used when the system is not. We consider the autonomous case first
x' = f ( x, y )
y' = g ( x, y )
See Ledder, Chapter 5, Section 3.
Fields and flows
Either DSolve (or NSSolve) will output time series for a two dimensional system. Trajectories can be obtained using the phi function (or Table) with ParametricPlot. If the system is autonomous PlotVectorField will also draw a direction field. Consider the following example derived from the Lotka-Volterra predator prey model (Ledder, page 306, Example 2).
In[2]:=
DEsystem = {x'[t] == x[t]*(1 - y[t]),
y'[t] == y[t]*(x[t] - 1)}
Out[2]=
This is a non-linear system and exact solution formulas are very hard to obtain. Mathematica's DSolve function is unable to obtain solutions for the initial conditions x(0) = 0.4 and y(0) = 0.4..
In[3]:=
DSolve[ {DEsystem, x[0]==0.4, y[0]==0.4}, {x, y}, t ]
Out[3]=
We turn to NDSolve, requesting a numerical solution for the interval from t = 0 to t = 20..
In[4]:=
soln = NDSolve[ {DEsystem,x[0]==0.4,y[0]==0.4}, {x, y}, {t,0,20} ]
Out[4]=
The following plot displays the time series, x(t) and y(t). We get fancy with informative axes labels.
In[5]:=
Plot[ Evaluate[{x[t],y[t]}/.soln], {t,0,20}, AspectRatio->1/4,
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,1]}},
AxesLabel->{"Time (years)","x (red) y (blue)"} ]
The system is autonomous so PlotVectorField can be used to display the vector field:
At point ( x , y ) plot the vector < ( x (1 - y ) , y (x - 1 ) >
The plot is named "vf" so it can be displayed with the trajectory below..
In[6]:=
vf = PlotVectorField[ {x*(1-y),y*(x-1)}, {x,0,2.8}, {y,0,2.8}, Axes->True, PlotLabel->"Vector field for DEsystem" ]
And here is the trajectory, making its way through the vector field. The solution is periodic.
In[7]:=
Show[ vf,
ParametricPlot[ Evaluate[{x[t],y[t]}/.soln], {t,0,20},
PlotStyle->RGBColor[0,1,0] ],
PlotLabel->"Field and Flow" ]
Estimating the period: Numeric data
We can estimate the period using a table of numeric data generated from the NDSolve solution. It appears from the time series plot above that the solutions start to repeat at about t = 7. The following Table pins it down to 6.96.
In[57]:=
Join[ {{t,x[t],y[t]}},
Evaluate[ Table[ {t,x[t],y[t]}, {t,6.9,7,0.01}]]/.soln[[1]] ]
//MatrixForm
Out[57]//MatrixForm=
A more complete flow
A more complete picture of the flow is obtained with a family of numeric solutions. The following input defines four initial values, creates four numeric solutions,and then shows their parametric plots in the vector field. The picture is named flow so it can be shown again below.
In[8]:=
inits = {0.4, 0.6, 1.3, 1.8};
solns = Evaluate[
Table[ NDSolve[ {DEsystem,x[0]==inits[[k]],y[0]==inits[[k]]},
{x,y}, {t,0,20} ], {k,4} ] ];
Show[ vf,
ParametricPlot[ Evaluate[Table[{x[t],y[t]}/.solns[[k]], {k,4}]],
{t,0,20}, PlotStyle->RGBColor[0,1,0] ],
PlotLabel->"Field and more flow lines",
AspectRatio->1/1 ]
flow = %;
Nullclines, critical points
The nullclines for an autonomous system are the curves determined by the equations
f (x, y) = 0
g (x, y) = 0
The points where the nullcllines cross are called the critical points or the stationary points of the system.
The ImplicitPlot function is used to draw the nullclines for DEsystem. The critical point in the first quadrant is also drawn using the ListPlot function. These plots are given names (output is suppressed) and shown with the above flow plot. We first load the ImplicitPlot function
In[12]:=
<< Graphics`ImplicitPlot`
In[22]:=
nullclines = ImplicitPlot[ {x*(1-y)==0, y*(x-1)==0}, {x,0,3}, {y,0,3} ];
stationary = ListPlot[ {{1,1}}, PlotStyle→PointSize[0.03] ];
Show[ nullclines, stationary, flow, AspectRatio->1/1 ]
Symbolic solutions, integral curves
The Lotka-Volterra system is non-linear. Mathematica cannot obtain a symbolic solution for the trajectories. However, formulas for curves determined by the trajectories can be obtained by eliminating the t variable to make a first order differential equation for the y variable as a function of x. Solutions to this differential equation are called "integrals" of the system, their graphs are called integral curves. Integral curves for Lotka-Voltera can easily be found "by hand" (Ledder, page 306). Mathematica's computations are shown in the next few entries.
Create DE:
In[26]:=
DE = y'[x] == DEsystem[[2,2]]/DEsystem[[1,2]]
Out[26]=
Change y[t] and x[t] into y[x] and x respectively.
In[27]:=
DE = DE/.{x[t]->x, y[t]->y[x]}
Out[27]=
Solve DE,with the initial condition y(0,4) = 0.4. to see if DSolve is up to the job.
In[41]:=
soln = DSolve[ {DE, y[0.4]==0.4}, y, x ]
Out[41]=
Mathematica has found an explicit solution in terms of the "ProductLog" function.
In[52]:=
?ProductLog
The plot of the solution shows that DSolve has obtained the formula for the bottom half of the integral curve for the system.
In[46]:=
Plot[ y[x]/.soln, {x,0,3}, PlotRange->{0,2.5}, AspectRatio->1/1 ]
Linear systems, linear flows
Two dimensional linear systems have the form
x' = a x + b y + f
y' = c x + d y + g
where x, y, a, b, c, d, f, and g are all functions of t. The functions f and g are sometimes referred to as the forcing functions. If they are both zero, the system is called homogeneous (or unforced).
When a, b, c, and d are constants, symbolic solutions can be always be found, provided the functions f and g are elementary (products of polynomials, sine, cosine, exp, etc.). Such systems are important not only because they arise in many applications, but also because they are used to approximate non-linear systems in regions near their critical points.
Example Consider the autonomous linear system
x' = x - 2 y + 3
y' = x + y - 2
Sketch the flow, nullclines and stationary point. Find symbolic solutions.
In[53]:=
DEsys = { x'[t] == x[t] - 2*y[t] + 3,
y'[t] == x[t] + y[t] - 2 }
Out[53]=
Let's find the stationary point first. The Solve function can be used. Note that the equations and the unknowns must be entered in separate lists. The output is a set, or sets, of arrow substitutions..
In[54]:=
stationary = Solve[ {x-2*y+3==0, x+y-2==0}, {x,y} ]
Out[54]=
Use the solution to make a "point" (i.e, a list for plotting, as follows.
In[55]:=
statpoint = {x,y}/.stationary
Out[55]=
The next entries create the plots for vector field, the flowlines (using NDSolve), the nullclines, and the stationary point. They are all shown together using Show.
In[143]:=
vf = PlotVectorField[{x-2*y+3, x+y-2}, {x,-4,4}, {y,-1,4} ];
solns = Table[ NDSolve[ {DEsys, x[0]==0, y[0]==k}, {x,y}, {t,-2,3} ],
{k,0,4} ];
flow=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.solns[[k,1]],{k,1,5}] ],
{t,-2,3}, PlotStyle->RGBColor[0,1,0] ];
nullclines = ImplicitPlot[ {x-2*y+3==0,x+y-2==0}, {x,-4,4}, {y,-1,4} ];
stationarypoint = ListPlot[ statpoint, PlotStyle->PointSize[0.03] ];
Show[ {vf, flow, nullclines, stationarypoint}, Axes->True,
PlotRange->{{-4,4},{-1,4}}, AspectRatio->5/8 ]
This kind of outward spiraling autonomous flow is called a "source". Phase portraits for constant coefficient autonomous linear systems are often referred to as "linear flows".
Symbolic solution formulas for a source are given in terms of sines, cosines, and exponential functions.
In[157]:=
solnFormulas = DSolve[ DEsys, {x,y}, t ]
Out[157]=
These formulas are easier to understand if displayed in simplified form.
In[158]:=
x[t]/.solnFormulas//FullSimplify
Out[158]=
In[159]:=
y[t]/.solnFormulas//FullSimplify
Out[159]=
Observe that the "stationary solution" is obtained by setting the constants equal to 0. Linear flows are discussed in more detail in the next section.
Go to Next Section Previous Section
Created by Mathematica (October 5, 2004)