Example 6

Discontinuities in the solution itself complicate matters. If nothing else, we must specify the jumps. In this example we show how to deal with the special case of a jump in the solution at the initial point. We also show how to plot the solution in a phase plane.

We solve a model of the infamous four year life cycle of a population of lemmings [22]. The equation
y(t)
=
r y(t)

1  -   y(t-0.74)
m


is solved on [0,40] with history y(t) = 19 for t < 0. The parameters r and m have the values 3.5 and 19, respectively. Notice that with these values, the equation has a constant (steady-state) solution, y(t) = 19. Tavernini uses this solution as history and perturbs the initial value to y(0) = 19.00001 so that y(t) will move away from the steady state. Here we use y(0) = 19.001 so as to see the cyclic behavior sooner.

Because the solution settles into a cyclic behavior, it is interesting to plot y(t) against y(t). This is easily done because in addition to sol.y, dde23 also returns a field sol.yp with values of the first derivative. For this example, these values provide an acceptable graph, but if they did not, we could get one by using ddeval to obtain more values for both y(t) and y(t).

Most DDE problems have solutions that are continuous at the initial point, so there is no need to supply the solver with an initial value in addition to a history function. However, if you should want to use a different initial value, all you have to do is provide it as the value of the 'InitialY' option. The solver deals automatically with the discontinuity in the first derivative that is ordinarily present at the initial point, so you need act only if the solution itself is discontinuous. Here the solution has a small jump at the initial point, indeed small enough that we must use error tolerances smaller than the default values so that the solver ``sees'' the jump.

Using the capability of passing parameters through dde23, exam6f.m can be coded as

   function v = exam6f(t,y,Z,r,m)
   v = r*y*(1 - Z/m);
The complete program exam6.m to compute and plot the solution is then
 r = 3.5;
 m = 19;
 options = ddeset('RelTol',1e-4,'AbsTol',1e-7,...
                  'InitialY',19.001);
 sol = dde23('exam6f',0.74,19,[0 40],options,r,m);

 plot(sol.x,sol.y);
 title('Figure 6a. Population of Lemmings--Time Series')
 xlabel('time t');
 ylabel('y(t)');
 figure
 plot(sol.y,sol.yp);
 title('Figure 6b. Population of Lemmings--Phase Plane')
 xlabel('y(t)');
 ylabel('y''(t)');

images/6a.gif
Clearly the normalized population gets quite small, but how small? A reasonably accurate answer is obtained easily: The smallest value of y(t) is approximately min(sol.y), namely 0.0116. If we wanted a better answer, we could obtain it by introducing event functions as in the last example.

images/6b.gif

Reference

[22]
L. Tavernini, Continuous-Time Modeling and Simulation, Gordon and Breach, Amsterdam, 1996.