Subsections

• 9.1 General remarks
• 9.2 To track or not to track discontinuities of derivatives of DDE solutions?
• 9.3 Runge-Kutta-like numerical methods

9 Numerical algorithms of the toolbox

9.1 General remarks

At present theoretical aspects of DDE theory are elaborated with almost the same completeness as the corresponding parts of ordinary differential equations (ODE) theory. However, unlike ODE, even for linear DDE there are no general methods of finding solutions in explicit forms.

Various specific numerical methods are constructed for solving specific delay differential equations. Most of investigations are devoted to numerical methods for systems with discrete delays and Volterra integro-differential equations. Consequent development of DDE numerical analysis and the corresponding bibliography is reflected in [11,6,2,4,5,3] and the corresponding chapters of the books [16,15].

For spesific classes of DDE there were elaborated special codes: [1,9,12,17,29,30,34].

Numerical algorithms, relized in Toolbox, are based on an approach to DDE numerical methods elaborated in [19,20,22].

The distinguishing features of the described numerical methods are the following:

1)

the numerical methods for DDE are direct analogies of the corresponding classical numerical methods of ODE theory, i.e., if delays disappear, then the methods coincide with ODE methods ¹ ;

2)

the proposed numerical algorithms do not depend on specific forms of delays, so the corresponding numerical schemes are the same for different classes of DDE.

In the frameworks of the approach one can construct for DDE analogies of all known numerical methods of ODE case.

Of course, the numerical algorithms relized in Toolbox are substantiated under some assumptions on general form of DDE and initial data, however these assumptions are satisfied for many realistic initial value problems.

Specific numerical algorithms of solving some specific delay differential equations can be more efficient than the algorithms realized in Toolbox. Nevertheless in the frameworks of the proposed approach the numerical schemes are the same for different classes of DDE and this fact allows us to elaborate sufficiently simple structure of software algorithms and packages.

An approach to overcoming the problem of handling the discontinuities of derivatives of DDE solutions will be discussed in the next subsection.

9.2 To track or not to track discontinuities of derivatives of DDE solutions?

Exact (analytical) solutions of DDE can have discontinuities of derivatives which can affect the numerical algorithms used for their approximate solving. So at present a big attention is paid to the problem of constructing special numerical algorithms in the regions of discontinuities. Actually, on our opinion, there can be (at lest) two points of view on this problem:

(A)

to construct special algorithms that handle possible discontinuities of solutions;

(B)

to construct numerical algorithms assuming smoothness of DDE solutions on whole interval.

In Toolbox we choose the second way (B) by the following reasons

$\bullet$

for a specific DDE, an initial function can be approximated, as rule, by a sequence of (initial) functions which generate smooth solutions ²,

$\bullet$

numerical experiments showed that numerical algorithms, realized in Toolbox, are robast with respect to discontinuities of derivatives of DDE solutions.

So, if the construction of DDE numerical methods is not the end in itself but only the mean of investigating DDE then the assumption of smoothness of solutions is acceptable in many cases and is a suitable compromise. Assumption of smoothness of solutions is one of the key points of our approach to realization of the numerical algorithms in the present version of Toolbox.

Of course, there are problems in which one cannot avoid discontinuities of derivatives of solutions. We reffer readers to the works [2,4,5,3] for effective algorithms and furthere references.

9.3 Runge-Kutta-like numerical methods

In this subsection we describe the numerical methods [19,20,22] realized in Toolbox.

For the sake of simplicity let us consider the uniform partition t_k = t₀ + k$\Delta$, k = 0, 1,..., N, of the interval [t₀, t₀ + $\theta$] (here $\Delta$ = $\theta$/N), and assume that the ratio $\tau$/$\Delta$ = N_$\tau$ is a positive integer.

Let u_k $\in$ $\bf R^{n}_{}$ be the approximation of the true solution x(t_k) = x_k of system (5.7), (5.15), (5.16) at time t_k.

The set of values {u_i : i$\ge$0, k - N_$\tau$$\le$i$\le$k} ³will be denoted by {u_i}_k. This set define at time t_k future behavior of the discrete model.

In order to calculate u_{k + 1} by means of {u_i}_k it is necessary at time t_k to construct :

1)

an interpolation u_tk( ^. ) = {u(t_k + s), - $\tau$ $\leq$ s < 0 } of the pre-history {u_i}_k of the discrete model on the interval [t_k - $\tau$, t_k];

2)

an extrapolation u_{tk + a}( ^. ) (for a constant a > 0) of the pre-history {u_i}_k of the discrete model on the interval [t_k, t_k + a].

Definition 9.1. For a constant a > 0 a mapping I : {u_i}_k$\to$u( ^. ) $\in$ Q[t_k - $\tau$, t_k + a) is called the interpolational-extrapolational operator (IE-operator) of the model pre-history.

$\Box$

We will make the following assumption: IE-operator I({u_i}_k) satisfies the Lipschitz condition

$$\le \max_{k - N_\tau \le i \le k}^{}$$ (9.1)

for t $\in$ [t_k - $\tau$, t_k + a), pre-histories {u⁽¹⁾_i}_k and {u⁽²⁾_i}_k, u⁽¹⁾( ^. ) = I({u⁽¹⁾_i}_k), u⁽²⁾( ^. ) = I({u⁽²⁾_i}_k).

Definition 9.2. The interpolational-interpolational operator I has an order p on a class of functions X $\subset$ Q[t_k - $\tau$, t_k + a) if there exist a constant C₁ such that for all x( ^. ) $\in$ X, discrete pre-history {x_i = x(t_i), k - N_$\tau$$\le$i$\le$k}, $\bar{x}$( ^. ) = I({u_i}_k) and t $\in$ [t_k - $\tau$, t_k + a) the following inequality is satisfied

$$\bar{x} \le \Delta^{p}_{}$$ (9.2)

$\Box$

Interpolation by spline of (p - 1)-order is one of methods of pre-history interpolation [20,22]. In this case one can realize the extrapolation by continuation of the last polynomial to the right side of point t_n. For this method of interpolation and extrapolation the Lipschitz condition (9.1) is satisfyed; moreover if the set X consists of p-times continuous differentiable functions then the operator I has the order p on the class X.

m-Stage explicit Runge-Kutta method (with an IE-operator I) is the discrete model

u₀ = x₀; (9.3)

u_i = y⁰(t₀ - t_i), i = - m,..., - 1, (9.4)

$$\Delta \sum_{i=1}^{m} \sigma_{i}^{}$$ (9.5)

$$\Delta \Delta \sum_{j=1}^{i-1} \Delta$$ (9.6)

h₁(u_k, u_tk( ^. )) = f (t_k, u_k, u_tk( ^. )) . (9.7)

The pre-history of the model is defined as u_t(s) = {I({u_i}_k) for  t_k - $\tau$$\le$t + s < t_k + a$\Delta$}, a = max{| a_i|, 1$\le$i$\le$m}.

Let us consider the functions z(t) = u(t) - x(t) and z_k = z(t_k) = u_k - x_k. One can derive

$$\Delta \psi^{(1)}_{k} \psi^{(2)}_{k}$$ (9.8)

where

$$\psi^{(1)}_{k} {\frac{x_k - x_{k+1}}{\Delta }} \sum_{i=1}^{m} \sigma_{i}^{}$$ (9.9)

$$\psi^{(2)}_{k} \sum_{i=1}^{k} \sigma_{i}^{}$$ (9.10)

The function $\psi^{(1)}_{k}$, k = 0, 1,..., N, is called the residual or the approximation error of method (9.3) - (9.7).

Definition 9.3. Method (9.3) - (9.7) :
1) has the approximation order p, if |$\psi^{(1)}_{k}$|$\le$C$\Delta^{p}_{}$ for all k = 0,..., N ;
2) converges, if z_k$\to$ 0 as N$\to$$\infty$ for all k = 1,..., N;
3) has the convergence order p if | z_k|$\le$C$\Delta^{p}_{}$ for all k = 0, 1,..., N. $\Box$

Theorem 9.1. Let method (9.3) - (9.7) have an approximation order p > 0 and an IE-operator I has the order p on a class of functions X (which contains the true solution of problem (5.7), (5.15), (5.16)), then Runge-Kutta-like method (9.3) - (9.7) converges and has the convergence order p. $\Box$

Conditions that determine approximation orders of numerical methods can be obtained by the expansion of the exact solution x(t) and the functionals h_i into the Taylor series at point t_k [22].

It is important to note that methods (9.3) - (9.7) are complete analogies of corresponding numerical methods of ODEs, that is approximation orders are the same for ODEs and DDEs with the same coefficients of numerical methods. For example, for DDE (5.7) the Runge-Kutta-like numerical schemes of the order 4 has the form

u_{n + 1} = u_n + $${\textstyle\frac{1}{6}}$$$$\Delta$$$$\left(\vphantom{ h_1 + 2 h_2 + 2 h_3 + h_4 }\right.$$h₁ + 2h₂ + 2h₃ + h₄$$\left.\vphantom{ h_1 + 2 h_2 + 2 h_3 + h_4 }\right)$$,

h₁ = f (t_n, u_n, u_tn( ^. )) ,

h₂ = f (t_n + $${\frac{\Delta}{2}}$$, u_n + $${\frac{\Delta}{2}}$$h₁, u_{tn + ${\frac{\Delta}{2}}$}( ^. )) ,

h₃ = f (t_n + $${\frac{\Delta}{2}}$$, u_n + $${\frac{\Delta}{2}}$$h₂, u_{tn + ${\frac{\Delta}{2}}$}( ^. )) ,

h₄ = f (t_n + $$\Delta$$, u_n + $$\Delta$$h₃, u_{tn + $\Delta$}( ^. )) .

Note, all results are valid for non-uniform partition of the time interval [t₀, t₀ + $\theta$]. It allows to use variable steps for constructing numerical models.

In Toolbox the Runge-Kutta-Fehlberg-like numerical schemes of the order 4(5) with automatic step size control and the interpolation and extrapolation by cubic splines are realized.