Interval methods of Adams-Bashforth type with variable step sizes

In a number of our previous papers, we have proposed interval versions of multistep methods (explicit and implicit), including interval predictor-corrector methods, in which the step size was constant. In this paper, we present interval versions of Adams-Bashforth methods with a possibility to change step sizes. This possibility can be used to obtain interval enclosures of the exact solution with a width given beforehand.

Interval arithmetic (see, e.g., [12,31,32,37]) realized in floating-point computer arithmetic is a way to estimate two first kinds of errors. Applying interval methods to approximate the solution of the initial value problem in floating-point interval arithmetic (see, e.g., [11]) we can obtain enclosures of the solution in the form of intervals which contain both these errors and also the truncation error.
We can distinguish three main kinds of interval methods for solving the initial value problem: methods based on high-order Taylor series (see, e.g., [1,2,6,15,33,34]), explicit and implicit methods of Runge-Kutta type ( [9,10,22,28,30,37]), and explicit and implicit multistep methods ( [9, 17-19, 22, 25, 26]), including interval predictor-corrector methods [27]. In interval methods based on Runge-Kutta and in interval multistep methods considered so far, a constant step size has been used. The methods based on high-order Taylor series use variable step sizes and seem to be the most universal. But it should be noted that in [25][26][27] and [30] we have shown that in some cases the interval methods of the second and third kinds give better enclosures of the exact solutions. Therefore, it is worth to take into account also such methods.
In conventional multistep methods, a problem for changing step size was described firstly by F. Ceschino in [4], and developed further, among others, by C. V. D. Forrington [8] and F. T. Krogh [20]. In these methods by appropriate selection of step size, one can minimize the truncation error, i.e., obtain an approximate solution at some points (mesh points of a grid) with a tolerance given beforehand. Since in interval methods the truncation errors are included in interval enclosures of the exact solutions, the only reason to change a given step size is to decrease the widths of these enclosures. Although we consider only interval methods of Adams-Bashforth type, it seems that the presented procedure can be also applied to other interval multistep methods.
The paper is divided into seven sections. In Section 2, we shortly recall conventional Adams-Bashforth methods with variable step sizes. Interval versions of these methods are proposed in Section 3. Section 4 is the main section of this paper, in which we describe how to match the step size to obtain interval enclosures of the exact solution with a desired width. Since for simplicity in Sections 2-4 we consider only one-dimensional initial value problem, in Section 5, we announce some remarks on solving systems of differential equations by interval Adams-Bashforth methods using our algorithm. In Section 6, we present some numerical examples, and in the last section, some conclusions are given.

Adams-Bashforth methods
Below we present briefly (on the basis of [11,p. III.5]) the methods of Adams-Bashforth (explicit multistep methods) 1 with variable step sizes. 1 The term "Adams-Bashforth method" is often used only for the four-step explicit linear method, which is a special case of the class of methods that we consider in this paper. In [11], the methods considered here are called "explicit Adams methods," but J. C. Butcher in [3, p. 143] used the same term as we use, since there are other similar linear explicit methods known as Adams-Nyström or simply Nyström methods (see, e.g., [3] or [13]).
( 1 ) Let the set of points {t 0 , t 1 , . . . , t k , . . . , t m }, such that 0 = t 0 < t 1 < . . . < t k < . . . < t m ≤ a, be a grid on the interval [0, a] (the points t i are called the mesh points), and h i = t i −t i−1 (i = 1, 2, . . . , m) denote the step sizes. On the interval t k−1 , t k , the differential equation (1) is equivalent to The integrand in (2) we approximate by the interpolation polynomial P (τ ) of degree n − 1 in such a way that P t k−j = f t k−j , y t k−j , j = 1, 2, . . . , n.
Substituting τ = t k−1 + th k and using Newton's interpolation formula, we can write the polynomial P (τ ) in the form where t k−1 , t k−2 , . . . , t k−j −1 ; f denote the divided differences given by which may be defined recursively by 2 The case when (1) presents a system of differential equations is discussed in Section 5 Let us denote Then, the formula (2) may be written as where r n (τ ) denotes the interpolation error, and the integral over this error is the local truncation error. We have (see, e.g., [14,35]) where ξ ∈ t k−n , t k and y (n+1) (ξ ) = f (n) (ξ, y (ξ )). Thus, the formula (6) may be written in the form where the values ϕ j (k) and g j (k) for j = 0, 1, . . . , n − 1 are given by (4) and (5), respectively, and The values ϕ j (k) and g j (k) for j = 0, 1, . . . , n − 1 can be computed efficiently with recurrence relations [11,13]. As an immediate consequence of (3), we have where the coefficients can be calculated by The calculation of the coefficients g j (k) is more tricky [13,21]. We introduce the q-fold integral It appears that where for the quantities c jq (t k ) we have the following recurrence relations: If the approximations y k−1 , y k−2 , . . . , y k−n of the exact values y (t k−1 ) , y (t k−2 ) , . . . , y (t k−n ) are known, we denote f i = f (t i , y i ) for i = k − n, k − n + 1, . . . , k − 1, and in the formula (7) we omit the truncation error, then we obtain the Adams-Bashforth methods with variable step sizes given by In particular, from (9) we get: where where where If for each i = k, k − 1, . . . , k − n + 1 we have h i = h (a constant step size), then from the above formulas we obtain the following well-known conventional Adams-Bashforth methods (see, e.g., [3,5,11,16,38]): -F (T , Y )-an interval extension of f (t, y), where an interval extension of the function where IR denotes the space of real intervals.
and let us assume that: -the function F (T , Y ) is defined and continuous 3 for all T ⊂ t and Y ⊂ y , -the function F (T , Y ) is monotonic with respect to inclusion, i.e., -for each T ⊂ t and for each Y ⊂ y , there exists a constant > 0 such that where w (A) denotes the width of the interval A, -the function (T , Y ) is defined for all T ⊂ t and Y ⊂ y , -the function (T , Y ) is monotonic with respect to inclusion.
The interval methods of Adams-Bashforth type with variable step sizes we define as follows (compare (7)): where g j (k) for j = 1, 2, . . . , n − 1 is given by (5), g n (k)-by (8), and In particular, for a given n, we get the following methods (compare (10)-(13)): where where 3k = 2k For the methods (14), we can prove that the exact solution of the initial value problem (1) belongs to the intervals (enclosures) obtained by these methods. Before that, it is convenient to present the following This implies that From (15) and (21), the relation (20) follows immediately.
. . , n − 1, then for the exact solution y(t) of the initial value problem (1) we have are obtained from the methods (14).
Proof Let us consider the formula (7) for k = n. We have where ξ ∈ [t 0 , t n ]. From the assumption we get y (t n−1 ) ∈ Y n−1 , and from the Lemma it follows that Applying Taylor's formula, we get where Moreover, Taking into account the above considerations, from (24), we get As we assumed, is an interval extension of f (n) (t, y). Hence, applying (25) and (26), we have y , but-according to (14)-this is the interval Y n . This conclusion ends the proof for k = n. In a similar way we can show the thesis of this theorem for k = n + 1, n + 2, . . . , m.
, and the intervals Y k for k = n, n + 1, . . . , m are obtained from (14), then and the constants A, B, and C are independent of h.
Proof Because of g j (k) > 0 for j = 0, 1, . . . , n − 1 (see (5)) and g n (k) > 0 (see (8)), from (14), we get We assumed that is monotonic with respect to inclusion. Moreover, if step sizes h k−j for j = 0, 1, . . . , n − 1 are such that satisfy the conditions From the above inclusion, it follows that If we denote then from (15) we have But we also assumed that for the function F there exists a constant > 0 such that Therefore, from (29)-(32), we get then from the above inequality it follows that we have But Thus, from (34), we get where From (35) for k = n, we obtain and using mathematical induction, one can prove that for each where P ni = 1 + nς n (h n+i ) 1 + P n,i−1 , Q ni = 1 + nς n (h n+i ) Q n,i−1 + nς n (h n+i ) , R ni = 1 + nς n (h n+i ) R n,i−1 + σ n (h n+i ) , P n0 = 1 + nς n (h n ) , Q n0 = nς n (h n ) , R n0 = σ n (h n ) .
It is self-evident that It means that from (38), we have Since i ≤ m − n, then from (40), we finally get

Using variable step sizes
In conventional one-and multistep methods, the change of step size is used to obtain approximate solutions with a given tolerance (to decrease the truncation errors). In the case of multistep methods, it is worth to recommend the methods of Adams with variable step size of predictor-corrector type, in which the explicit methods of Adams-Bashforth are predictors and the implicit methods of Adams-Moulton are correctors. A usable algorithm based on these methods has been proposed by F. T. Krogh in [21], and appropriate procedures in Fortran programming language to implement this algorithm on computers have been published in the book [36] written by L. F. Shampine and M. K. Gordon.
In interval methods, all errors (representation, rounding, and truncation errors) are included in obtained enclosures of solutions. Thus, the only aim to choose an adequate step size is to decrease the width of interval enclosures.
To obtain interval enclosures of Y k with a width eps giving beforehand, one can consider the inequality (27) and compare the right-hand side of it. But this inequality overestimates the width of Y k , and a better way to do it is to take into consideration the relation (35).
Let us assume that the right-hand side (34) is equal to eps, i.e., Since T k−j (j = 1, 2, . . . , n) is only an interval representation of t k−j , we have w T k−j ≈ 0. If we assume w T k−j = 0 and take into account (42), then we have The equation (43) presents a polynomial (with respect to h) of degree n+1. Denoting this polynomial by p n+1 (h k ), we look for h k > 0 such that p n+1 (h k ) = 0. To find such an h k , one can use any method for finding roots of polynomials. Applying, for example, the Newton iterations we have where h (0) k is given (we can take h (0) k = h k−1 ). The process (44) is stopped when where ε denotes an accuracy given beforehand. Note that in general, the real h k > 0, satisfying (43), may not exist. Since in (43) all coefficients near by any power of h k are positive, such a case occurs, for instance, when w (Y k−1 ) − eps > 0.
If this relation holds, we can try to increase eps or simply stop the calculations. Of course, if the widths of successive Y k increase, then the successive step sizes h k must decrease (see Example 2 in Section 6). Moreover, if t m is the last t k , for which t k < a, and we want to achieve the end a of integration interval, then for the last integration step we should take h m+1 = a − t m .
According to the definitions of g n (k) and ρ n (k) (see (8) and (33), respectively), it is difficult to write the general form of p n+1 (h k ). For n ≤ 4, we have the following formulas: Note that in the case of n > 1, the form of p n+1 (h k ) depends on which α j (k) (j = 0, 1, . . . , n − 1) is chosen (see (33)).

A note on a system of differential equations
For simplicity, we have considered only one-dimensional initial value problem. The methods (14) and the algorithm presented in Section 4 can be easily extended to systems of ordinary differential equations, i.e., to the initial value problem of the form where y ∈ R N and f : [0, a] × R N → R N . In this case the formula (14) and the relation (22)

Numerical examples
In the examples presented below, we have used our own implementation of floatingpoint interval arithmetic in Delphi Pascal. This implementation has been written in the form of a unit called IntervalArithmetic32and64, which current version is presented in [23]. This unit take advantage of the Delphi Pascal floating-point Extended type and makes it possible to represent any input numerical data in the form of a machine interval, perform all calculations in floating-point interval arithmetic, use some standard interval functions and give results in the form of proper intervals (if the ends of an interval are not the same machine numbers, one can see the difference in the output). All programs written in Delphi Pascal for the examples presented one can find in [24]. We have run these programs on Lenovo ® Z51 computer with Intel ® Core i7 2.4 GHz processor.
with the exact solution y = exp(0.5t).
(46) Let us assume , where x denotes the smallest machine number greater or equal to x (similarly, by x, we will denote the largest machine number less or equal to x). From (46), we obtain the results presented in Table 1 Table 2. The methods are very fast and for all of them the time of calculations (CPU time) is less then 1 s. Of course, for each method, the exact solution y(t 20 ) ≈ 2.7182818284590452 is within the interval obtained. It should be added that if we use the interval methods of Adams-Bashforth type for greater n and very small step sizes we can obtain intervals with greater widths then presented in Table 2 (see [22,Example 4.5 in Sec. 4.7]). This is caused by a great number of calculations in these methods and by a significant increase of rounding errors following from that, which is not compensated for the method orders.
Example 2 Again, let us take into account the problem (45), but now let us apply the procedure for finding step sizes described in Section 4. For the same starting data as in Example 1, requiring 10 −8 for interval widths and taking = 0.5, ε = 10 −18 , h  Table 3. The CPU times have been equal to 36.019 s for the method (16), 1.152 s for (17), 0.344 s  Fig. 1 Step size changes in problem (45) for the methods (16)- (19) for (18), and 0.415 s for the method (19). The changes of step sizes are shown in Fig. 1a-d. Let us note that it is not possible to achieve the assumed end t = 2 of integration interval. In each method, the step size is decreased for following integration steps (excluding initial steps in two-and more step methods, where some fluctuations can be caused by initial data), but there exists a point t m , after which further calculations are impossible. To counteract such a situation, one can increase the required widths of intervals or consider a limitation of integration interval. For instance, if for the problem (45) we assume t ≤ 0.6, then-for the same remaining initial data as previously-we obtain results presented in Table 4 (we present only intervals obtained at the last step for which t m < 0.6 and at the final t = 0.6, where the last step size is equal to 0.6 − t m ). In order to achieved the point t = 0.6, the methods (16), (17), (18), and (19) have needed 10.274 s, 0.574 s, 0.138 s, and 0.073 s, respectively. As we could expect, the greater order of method has permitted to obtain results in smaller numbers of integration steps.
Example 3 For the initial value problem of the form (the problem A5 from [7, p. 23 the solution is unknown. In [7], it is taken a = 20. In our methods, we cannot take this value, because for 0 ≤ t ≤ 20 we have −1 < y < 6.3 and in the interval extension F t , y of f (x, y) for t = {t ∈ R : 0 ≤ t ≤ 20} and y = {y ∈ R : −1 ≤ y ≤ 6.3} we have a division by an interval containing zero (such an operation is not defined in proper interval arithmetic). Thus, we restrict these regions to t = {t ∈ R : 0 ≤ t ≤ 10} , y = y ∈ R : 4 ≤ y ≤ 6.3 .
Let us take the fourth-order method (19) with additional starting intervals given in Table 5. These intervals have been obtained by an interval version of conventional Table 4 Intervals obtained for the problem (45) with t = t ∈ R : 0 ≤ t ≤ 0.6 and y = y ∈ R : 1 ≤ y ≤ 2.72  Runge-Kutta method (of fourth order) with an optimal step size for the first integration step (see [29] for details). For the method (19), we also need an interval extension of y (5) (t) = f (4) (t, y). Since for the problem (47), we have = h 3 as the initial approximation of step size for the first integration step, after 2.295 s, we have obtained the results presented in Table 6. Figure 2 shows the step size changes. Note that we are able to execute 464 integration steps and obtain the last interval for t ≈ 1.5476, but step sizes are small (and tend to zero) for t approximately greater than 1.54. If we assume eps = 10 −12 , then we achieve t ≈ 1.3991 (smaller t) after 4.758 s, and for eps = 10 −4 we can achieve t ≈ 2.0843 (greater t) after 3.333 s (see Tables 7 and 8 for details).
Example 4 Finally, consider the motion of a simple pendulum described by the equation where ϕ = ϕ (t), u = √ g/L, and where g is the gravitational acceleration at the Earth surface and L denotes the pendulum length. If we assume that the angle ϕ is small, i.e., sin ϕ ≈ ϕ, then the equation (48) can be reduced to the equation of simple     harmonic motion ϕ + u 2 ϕ = 0 (49) with the solution ϕ = ϕ 0 cos (ut), where ϕ 0 is an initial angle. Denoting y 1 = ϕ , y 2 = ϕ and assuming that ϕ (0) = 0, ϕ (0) = ϕ 0 , we can transform (49) into the following systems of differential equations of first order: with initial conditions y 1 (0) = 0, y 2 (0) = ϕ 0 .

Fig. 3
Step size changes in problem (50)-(51) for the method (18) 0.740 s. The changes of step size are shown in Fig. 3, from which it follows that for t approximately greater than 0.125 the step size is too small (tends to zero) to satisfactory continue the calculations. On the other hand, one can observe thataccording to Theorem 1-the exact solution belongs to interval enclosures obtained. The examples show that using the procedure for step size changing we can obtain interval enclosures of the exact solution with widths given beforehand for t max ≤ a. Unfortunately, very often we have t max << a. Of course, to reach a, we can always use a multistep method with constant step size and theorems similar to the Theorem 1 guarantee that the exact solution will be inside interval enclosures obtained. 5 We can also use a combination of methods with variable and constant step sizes. If we assume that the step size h cannot be less than some h min , we can start calculations using a method which matches step sizes with a width of intervals given beforehand, and when the step size will be decreased to h min -continue calculations with h = h min . Of course, such a procedure does not guarantee that at t ≈ a we will obtain intervals with a given width, but it can significantly reduce the number of integration steps.
Some remarks concerning computing time for our methods should be given at the end. First of all, let us note that the Adams-Bashforth methods with step size changes are very fast (usually in the process (44), only a few iterations are needed to achieve a desired width eps of interval). Obviously, for a greater n (the number of method steps and the order of method), the computing time is less for a given t ≤ a. Moreover, for decreasing desired widths of intervals, we can observe decreasing values of t max ≤ a, but computing times may not be decreasing (compare the values of t max and computing times in Example 3 for eps = 10 −4 , 10 −8 and 10 −12 ). To confirm this observation, let us consider  (17), (18), and (19) and different values of eps Example 5 For the same problem as in Example 3 and for the same t and y , let us take different values of eps (from 10 −2 to 10 −13 ) in the method (19). In Table 11, we present the number of steps k and CPU times to achieve t max for each value of eps. These times are also presented in Fig. 4, and we can see that for eps = 10 −8 the computing time is the smallest. It may be generalized: for each problem considered, there exists the smallest computing time for some eps.
On the other hand, for some t < T ≤ a, where T denotes the smallest t max common for different n and eps, the smaller computing times are obtained for greater values of n. For the problem considered, the methods (17)

Conclusions
Until now, interval multistep methods (explicit and implicit) have been used only with a constant step size. Using variable step sizes and the procedure presented in Section 4 we can obtain interval enclosure of the exact solution with a width given beforehand. Although we have presented our algorithm only for interval methods of Adams-Bashforth type, it seems that the same procedure one can apply to other multistep methods, including interval predictor-corrector ones [27].