A "Basic" implementation of Euler's method

The Real Maple Way to do it

by Prof. Matt Miller, Department of Mathematics, University of South Carolina
email: miller@math.sc.edu

We look at two simple models of population growth for a bunch of (pseudo-)rabbits. To distinguish them, we use R and S to denote the respective rabbit populations. Both begin with size 10, and both have intrinsic (per capita) rate of increase r = 0.1.

> restart:

> tinitial:= 0; tfinal:= 60 ; # set up initial and final times (in months)

> t:= tinitial; # start the clock at the initial time

> R:= 10 ; S:= 10 ; # set the initial values

> r:= 0.1; K:= 100; # parameter values

> numberofsteps:= 120 ; deltat:= (tfinal-tinitial) / numberofsteps;

> Rpoints:= [tinitial, R ]; Spoints:= [tinitial, S ];# set the initial point on the graph

The rate equation for S is called a logistic equation. We use Euler's method of stepwise approximation, taking the current value of S to compute the current rate S' = rS(K-S)/K, and from this the change in S, which we call DeltaS, and an updated value of S (which leads to an updated value of S' and the cycle continues). You can use a different rate equation, of course, for instance taking constant per capita growth (as we have done for R), or other more exotic models.

> for k from 1 by 1 to numberofsteps do # set up the counter for the cycle

> Rprime:= r * R : # exponential growth -- get the current rate

> Sprime:= r * S * ( K - S ) / K : # logistic growth -- get the current rate

> deltaR:= Rprime*deltat: R:= R + deltaR: # calculate change in R and update value of R

> deltaS:= Sprime*deltat: S:= S + deltaS: # same for S

> t:= t + deltat: # advance t by one step

The following command adds the point on the graph that has just been calculated to the roster of points that are there already. If you're not interested in the graph, you can "comment out" this command by placing a # in front.

> Rpoints:= Rpoints , [ t , R ] : Spoints:= Spoints , [ t , S ] :

> od: #commands between do and od are executed repeatedly until the counter ends the loop

Let's print out the final values of t, R, and S.

> t:= t; R:= R; S:= S;

Make a note of the result of this SIMULATION; that is, note all the parameter values used (final time, number of steps, rate equation coefficients) and the final population values. We are now going to change these values one at a time, and observe the results of the modified SIMULATIONS. What you are doing in effect is experimentation on the computer, and you need to keep track of what you are doing just as you would in a lab.

We can also plot R and S as functions of t. Again, comment this out if you aren't interested in the graph.

> plot([ Rpoints ], style = LINE, title = `Rabbit population: exponential growth`);

>

> plot([ Spoints ], style = LINE, title = `Rabbit population: logistic growth` );

>

Now select the menu item Edit, and Remove Output. Then go back to the beginning, change the number of steps, and then run through the whole thing again until you are satisfied that you have a sufficiently accurate final values of R and S. Next try changing the final time (modify the number of steps to keep the stepsize the same as before). What do you observe? Select a final time and stepsize that you are happy with (say t = 72 and numberofsteps = 288) and predict what effect changing the inital values of R and S will have. Test your prediction. Reset the inital values back to 10 and predict what effect changing r will have. Test your prediction. Reset r to 0.1 and predict what effect changing K will have. Test your prediction. Remember every time you make a change you need to reset all the other changes you have made (so that comparisons can be made) and execute the worksheet from the top!

Once you understand Euler's method, and what this program is all about, you may want to graduate to using Maple's buit-in commands for solving differential equations. There are a great many of them (change the : to a ; in the with(DEtools) command and you will see the list, or investigate the Help browser), but for now let's just see how we can plot of S as a function of t. All the information has to be specified, in the correct order: the rate equation (or equations), the variables (showing how the dependent variable(s) depend on the independent variable), the domain over which the independent variable (in this case t for time) runs, the initial condition (or conditions), the stepsize, and the way we want the solution plotted (the scene). Arrows are optional and come in a variety of styles; they indicate a sampling of the locally linear pieces of the graph that satisfies the rate equation at selected points. I usually find it handy to give some of these pieces of information names and lines of their own: this makes things a lot easier to modify later on without messing too much with the DEplot line itself.

> restart; with(DEtools): with(plots):

Warning, the name changecoords has been redefined

> rate_eq:= [ diff( S(t) , t ) = 0.1 * S * ( (100 - S) / 100 ) ];

> vars:= [S(t)] ; # what depends on what?

> domain:= 0 .. 60 ; # initial t to final t values (often called the time domain)

> init1:= [S(0)=10]; init2:= [S(0)=150]; init3:= [S(0)=30]; # initial values of S at t = 0

> # init1:= [0, 10]; init2:= [0, 150]; init3:= [0, 30]; # an alternate form

> DEplot(rate_eq, vars, domain, [init1], stepsize = 0.1, scene = [t, S], arrows = NONE);

One can give a DEplot a name (BE SURE TO END IT WITH A COLON IN THIS CASE!) to save it up to plot later. Usually this is done to get several plots showing together.

> fall:= DEplot(rate_eq, vars, domain,[init2], stepsize=0.1, scene=[t, S], arrows = NONE):

> rise:= DEplot(rate_eq, vars, domain, [init3], stepsize=0.1, scene=[t, S], arrows = NONE):

> display( [fall , rise ]);

>

Another way to accomplish the same end is to put all the initial conditions in one statement.

> DEplot(rate_eq, vars, domain, [init1, init2, init3], stepsize=0.1, scene=[t,S], arrows=NONE, linecolor = [BLUE, GREEN, BLACK]);

>

Now try all the various experiments that you did earlier, changing the time domain, step size, initial conditions, parameter values, and KEEP TRACK OF YOUR OBSERVATIONS! This is what we mean by "playing" with a model. What features are stable, which are highly variable (the model is "sensitive" to what, "insensitive" to what)?