Example 5

ODE and DDE solvers are intended for problems with solutions that have several continuous derivatives. It is not unusual for equations to have different forms in different circumstances, which leads to discontinuities in low-order derivatives of the solution, or even in the solution itself, when the circumstances change. Although a robust solver may be able to produce an acceptable solution, it is better practice to account for the changes and it can be necessary. There are two issues: Do we know in advance where the changes occur? Is the solution itself continuous? In this example we show how to solve problems that have a continuous solution with discontinuities in a low-order derivative at points known in advance. The history is the solution prior to the initial point and its discontinuities must also be taken into account because they propagate into the interval of integration. Discontinuities in the history are handled in the same way, but are a little simpler because discontinuities in the history itself are permitted.

Example 4.4 of [14] is an infection model due to Hoppensteadt and Waltman. The equation
y(t) =





-r  y(t)  0.4  (1 - t)
if 0 t 1-c,
-r  y(t) (0.4 (1-t) + 10 - em y(t))
if 1-c < t 1,
-r  y(t) (10 - em y(t))
if 1 < t 2 -c,
-r  em y(t) (y(t-1) - y(t))
if 2-c < t.
is solved on [0,10] with history y(t) = 10 for t 0. Here c = 1/2 and m = r/10. Oberle and Pesch solve this problem for several values of the parameter r, but we solve it only for r = 0.5. The different phases of the spread of the disease are described by different equations. In this example the phases change at times known in advance. The model requires the solution to be continuous, but the changes in the equation lead to jumps in low order derivatives. In addition to y(t), an approximation to I(t) = - y(t)/(r y(t)) is required.

dde23 deals easily with problems that have a continuous solution and discontinuities in low-order derivatives at known points. All you have to do is tell the solver where the discontinuities are by providing them as the value of the 'Jumps' option. However, you need to keep in mind that the history is the solution prior to the initial point, so you must also account for its discontinuities. For instance, the Marchuk immunology model discussed in [6] has the history max(0,t+10-6) for t 0. Its solution has a jump in the first derivative at t = -10-6 which propagates into the interval of integration. Discontinuities in the history are handled like discontinuities at known points during the integration. In one respect they are simpler; a jump in the history itself is treated the same as a jump in one of its low order derivatives. Low-order discontinuities in the history have an effect in the interval of integration because of the delays. If the initial point is a and the longest delay is T, discontinuities that occur before a - T have no effect on the integration, so there is no need to include them in 'Jumps'.

Having discussed how to deal with the discontinuities, it is straightforward to solve the problem. We compute an approximation to y(10) and compare it to an accurate value reported in [14]. This illustrates the computation of an approximation at a specific point and confirms the accuracy of the computation. We compute and plot I(t) at the points of sol.x using the fields sol.y and sol.yp. If we should want values at other t or should want a smoother graph, we would compute the necessary values with ddeval. If we treat r as a parameter, the equation can be coded as

   function v = exam5f(t,y,Z,r)
   c = 1/sqrt(2);
   mu = r/10;
   if t <= 1 - c
      v = -r*y*0.4*(1 - t);
   elseif t <= 1
      v = -r*y*(0.4*(1 - t) + 10 - exp(mu)*y);
   elseif t <= 2 - c
      v = -r*y*(10 - exp(mu)*y);
   else
      v = -r*exp(mu)*y*(Z - y);
   end
The complete program exam5.m is then
 r = 0.5;
 c = 1/sqrt(2);
 opts = ddeset('Jumps',[(1-c), 1, (2-c)],...
               'RelTol',1e-5,'AbsTol',1e-8);
 sol = dde23('exam5f',1,10,[0, 10],opts,r);

 y10 = ddeval(sol,10);
 fprintf('DDE23 computed     y(10) =%15.11f.\n',y10);
 fprintf('Reference solution y(10) =%15.11f.\n',0.06302089869);

 plot(sol.x,sol.y)
 title(['Figure 5a. Hoppensteadt-Waltman model with r = ',...
        num2str(r),'.'])
 xlabel('time t')
 ylabel('y(t)')

 Ioft = -(1/r)*(sol.yp ./ sol.y);
 figure
 plot(sol.x,Ioft)
 title(['Figure 5b. Hoppensteadt-Waltman model with r = ',...
        num2str(r),'.'])
 xlabel('time t')
 ylabel('I(t)')

This program results in the output

 DDE23 computed     y(10) =  0.06301980845.
 Reference solution y(10) =  0.06302089869.
and the two figures displayed. The accuracy of the computed result is what we might expect for the specified error tolerances.

images/5a.gif
images/5b.gif

References

[6]
E. Hairer, S.P. Nø rsett, and G. Wanner, Solving Ordinary Differential Equations I, Springer-Verlag, Berlin, 1987.

[14]
H.J. Oberle and H.J. Pesch, Numerical treatment of delay differential equations by Hermite interpolation, Numer. Math., 37 (1981) 235-255.