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);

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);

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));

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);

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

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
```

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);

> 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 can readily be used to find the equilibria.

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

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

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);

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);

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

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

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});

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);

> 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 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;

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

There are now 3 equilibria.