Exercise 7

The equations of the Marchuk immunology model discussed in [6] are
y1(t)
=
(h1 - h2 y3(t)) y1(t)
y2(t)
=
x(y4) h3  y3(t-t) y1(t-t) - h5 (y2(t) - 1)
y3(t)
=
h4 (y2(t) - y3(t)) - h8 y3(t) y1(t)
y4(t)
=
h6 y1(t) - h7  y4(t)
Here the coefficient
x(y4) =

1
if y4 0.1
(1 - y4) 10/9
if 0.1 < y4 1
is continuous, but has a jump in its first derivative when y4(t) = 0.1, which leads to a jump in a low order derivative of y2(t). The problem is solved on [0,60] with history y1(t) = max(0,t+10-6), y2(t) = 1, y3(t) = 1, y4(t) = 0 for t 0 . As was noted in Example 5, y1(t) has a jump in its first derivative at t = -10-6 that propagates into the interval of integration. Figure 15.8 of [6] presents plots for h1 = 2, h2 = 0.8, h3 = 104 , h4 = 0.17, h5 = 0.5, h7 = 0.12, h8 = 8 and two values of h6, namely 10 and 300. Treat h6 as a parameter in your program and try to reproduce the figure for h6 = 300. For this you will have to plot the scaled components 104y1, y2/2, y3, 10 y4 with axis([0 60 -1 15.5]) . An array yplot of scaled values for plotting can be formed easily by

 yplot = sol.y;
 yplot(1,:) = 1e4*yplot(1,:);
and so forth. To solve this problem accurately over the whole interval, you will need to reduce the tolerances to, say, a relative tolerance of 10-5 and an absolute tolerance of 10-8.

dde23 is sufficiently robust that it can solve this problem in a straightforward way. However, you can be much more confident of the results if you ensure that the solver is always working with equations that have smooth coefficients. There are two kinds of discontinuities in this problem. As in Example 5, use the 'Jumps' option to tell the solver about the discontinuity at t = -10-6. As in Example 7, terminate the integration when the event function y4(t) - 0.1 vanishes. Use a parameter state with value +1 if y4(t) 0.1 and -1 otherwise. The problem is to be solved with y4(0) = 0, so initialize state to +1. Thereafter, each time that the solver returns, check whether you have reached the end of the interval. If sol.x(end) < 60, change the sign of state and call dde23 again with the previous solution as history. In the function for evaluating the DDEs, set x(y4) = 1 if state is +1 and x(y4) = (1 -y4) 10/9 otherwise.

Reference

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