SDSU

 

Math 337 - Differential Equations
Spring Semester, 2002
Maple and Second Order Systems of Differential Equations

 © 2002, All Rights Reserved, SDSU & Joseph M. Mahaffy
San Diego State University -- This page last updated 08-Mar-02


Maple and Second Order Systems of Differential Equations

This section is designed to show how Maple can be used to numerically solve second order differential equations. It also shows how Maple provides tools to view the direction field of these equations and observe the qualitative behavior. The link here allows you to download the Maple worksheet that produced this webpage.

Consider a predator-prey system with rabbits, R(t), and foxes, F(t). The classic Lotka-Volterra system describing this interaction is given by the following system of differential equations:

> Rab := diff(R(t),t) = 2*R(t) - 1.2*R(t)*F(t);
Fox := diff(F(t),t) = -F(t) + 0.9*R(t)*F(t);

Rab := diff(R(t),t) = 2*R(t)-1.2*R(t)*F(t)

Fox := diff(F(t),t) = -F(t)+.9*R(t)*F(t)

The way you solve this system takes advantage of Maple's dsolve with the numeric option. The dsolve begins with the system of equations and initial conditions inside the set brackets. Then the variables solved are listed, which are followed by the numeric solving option and output being a procedure.

> sys := dsolve({Rab, Fox, R(0)=1,F(0) = 0.5}, {R(t),F(t)}, type = numeric, output = listprocedure);

sys := [t = proc (t) option `Copyright (c) 1993 by ...

Next we want to isolate the solutions for R(t) and F(t) for further analysis.

> pR := subs(sys, R(t));
pF := subs(sys, F(t));

pR := proc (t) local rkf45_s, outpoint, r1, r2; glo...

pF := proc (t) local rkf45_s, outpoint, r1, r2; glo...

We want to view the time varying solutions for the rabbits and foxes.

> Rabt := [seq([0.1*n, pR(0.1*n)], n = 0..200)]:
plot(Rabt, color = RED);

[Maple Plot]

> Foxt := [seq([0.1*n, pF(0.1*n)], n = 0..200)]:
plot(Foxt, color = BLUE);

[Maple Plot]

These solutions clearly oscillate in time. For both solutions together, we need the following:

> with(plots):
A := plot(Rabt, color = RED):
B := plot(Foxt, color = BLUE):
display({A, B});

Warning, the name changecoords has been redefined

[Maple Plot]

From this we see that the rabbit population grows first, followed by the fox population booming, which causes the rabbit population to crash, which later causes the fox population to fall. Then the cycle repeats itself.

Phase Portrait

Next we want to view the phase portrait. We will first view it using the numeric solution above, then use the special package in Maple that includes the direction field.

> RabFox := [seq([pR(0.1*n), pF(0.1*n)], n = 0..200)]:
plot(RabFox, color = MAGENTA);

[Maple Plot]

> with(DEtools):
DEplot([Rab,Fox],[R(t),F(t)], t=0..10, [[R(0)=1,F(0)=0.5]], R=0..4, F=0..5, title=`Lotka-Volterra model`, stepsize=.05, color = blue, linecolor = t/10);

[Maple Plot]

Maple can readily be used to find the equilibria.

> Req := 2*R - 1.2*R*F = 0;
Feq := -F + 0.9*R*F= 0;

Req := 2*R-1.2*R*F = 0

Feq := -F+.9*R*F = 0

> solve({Req,Feq},{R,F});

{R = 0., F = 0.}, {R = 1.111111111, F = 1.666666667...

Modified Lotka-Volterra Model

Consider a modification to the Lotka-Volterra model that includes logistic growth of the rabbits. We repeat most of the calculations above, and observe the differences for this system.

> RabM := diff(R(t),t) = 2*R(t)*(1 - R(t)/2) - 1.2*R(t)*F(t);
FoxM := diff(F(t),t) = -F(t) + 0.9*R(t)*F(t);

RabM := diff(R(t),t) = 2*R(t)*(1-1/2*R(t))-1.2*R(t)...

FoxM := diff(F(t),t) = -F(t)+.9*R(t)*F(t)

Once again we find the numerical solution to this system of equations.

> sysM := dsolve({RabM, FoxM, R(0)=1,F(0) = 0.5}, {R(t),F(t)}, type = numeric, output = listprocedure);

sysM := [t = proc (t) option `Copyright (c) 1993 by...

We isolate the solutions for R(t) and F(t) for further analysis.

> pRM := subs(sysM, R(t));
pFM := subs(sysM, F(t));

pRM := proc (t) local rkf45_s, outpoint, r1, r2; gl...

pFM := proc (t) local rkf45_s, outpoint, r1, r2; gl...

We view the time varying solutions for the rabbits and foxes.

> RabMt := [seq([0.1*n, pRM(0.1*n)], n = 0..200)]:

> FoxMt := [seq([0.1*n, pFM(0.1*n)], n = 0..200)]:

> AM := plot(RabMt, color = RED):
BM := plot(FoxMt, color = BLUE):
display({AM, BM});

[Maple Plot]

This shows a dramatically different behavior with the two populations leveling off at some equilibrium.

The Phase portrait is given by

> RabFoxM := [seq([pRM(0.1*n), pFM(0.1*n)], n = 0..200)]:
plot(RabFoxM, color = MAGENTA);

[Maple Plot]

> with(DEtools):
DEplot([RabM,FoxM],[R(t),F(t)], t=0..10, [[R(0)=1,F(0)=0.5]], R=0..1.5, F=0..1.2, title=`Lotka-Volterra modified model`, stepsize=.05, color = blue, linecolor = t/10);

[Maple Plot]

Maple can readily be used to find the equilibria.

> RMeq := 2*R*(1 - R/2) - 1.2*R*F = 0;
FMeq := -F + 0.9*R*F= 0;

RMeq := 2*R*(1-1/2*R)-1.2*R*F = 0

FMeq := -F+.9*R*F = 0

> solve({RMeq,FMeq},{R,F});

{F = 0., R = 0.}, {R = 2., F = 0.}, {R = 1.11111111...

There are now 3 equilibria.