Exercise 7

The equations of the Marchuk immunology model discussed in [6] are

y₁˘(t)

=

(h₁ - h₂ y₃(t)) y₁(t)
y₂˘(t)

=

x(y₄) h₃ y₃(t-t) y₁(t-t) - h₅ (y₂(t) - 1)
y₃˘(t)

=

h₄ (y₂(t) - y₃(t)) - h₈ y₃(t) y₁(t)
y₄˘(t)

=

h₆ y₁(t) - h₇ y₄(t)

Here the coefficient

x(y₄) = ě
í
î

1
if y₄ Ł 0.1

(1 - y₄) 10/9
if 0.1 < y₄ Ł 1

is continuous, but has a jump in its first derivative when y₄(t) = 0.1, which leads to a jump in a low order derivative of y₂(t). The problem is solved on [0,60] with history y₁(t) = max(0,t+10^-6), y₂(t) = 1, y₃(t) = 1, y₄(t) = 0 for t Ł 0 . As was noted in Example 5, y₁(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 h₁ = 2, h₂ = 0.8, h₃ = 10⁴ , h₄ = 0.17, h₅ = 0.5, h₇ = 0.12, h₈ = 8 and two values of h₆, namely 10 and 300. Treat h₆ as a parameter in your program and try to reproduce the figure for h₆ = 300. For this you will have to plot the scaled components 10⁴y₁, y₂/2, y₃, 10 y₄ 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 y₄(t) - 0.1 vanishes. Use a parameter state with value +1 if y₄(t) Ł 0.1 and -1 otherwise. The problem is to be solved with y₄(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(y₄) = 1 if state is +1 and x(y₄) = (1 -y₄) 10/9 otherwise.

Reference

[6]

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