See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/332118548

Interval Runge-Kutta Methods with Variable Step Sizes

Article · January 2019

DOI: 10.12921/cmst.2019.0000006

CITATIONS READS

2 89

2 authors:

Andrzej Marciniak Barbara Szyszka

Poznan University of Technology Poznan University of Technology
54 PUBLICATIONS   272 CITATIONS    17 PUBLICATIONS   88 CITATIONS   

SEE PROFILE SEE PROFILE

All content following this page was uploaded by Andrzej Marciniak on 09 July 2019.

The user has requested enhancement of the downloaded file.

CMST 25(1) 17–29 (2019) DOI:10.12921/cmst.2019.0000006

Interval Runge-Kutta Methods with Variable Step Sizes

A. Marciniak1,2 , B. Szyszka3

1
Institute of Computing Science
Poznan University of Technology
Piotrowo 2, 60-965 Poznan, Poland
E-mail: [email protected]
2
Department of Computer Science
State University of Applied Sciences in Kalisz
Poznanska 201-205, 62-800 Kalisz, Poland
3
Institute of Mathematics
Poznan University of Technology
Piotrowo 3A, 60-965 Poznan, Poland
E-mail: [email protected]

Received: 25 February 2019; revised: 28 March 2019; accepted: 28 March 2019; published online: 31 March 2019

Abstract: In a number of our previous papers we have presented interval versions of Runge-Kutta methods (explicit and
implicit) in which the step size was constant. Such an approach has required to choose manually the step size in order
to ensure an interval enclosure to the solution with the smallest width. In this paper we propose an algorithm for choosing
automatically the step size which guarantees the best (i.e., the tiniest) interval enclosure. This step size is determined with
machine accuracy.
Key words: initial value problem, Runge-Kutta methods, interval Runge-Kutta methods, variable step size, floating-point
interval arithmetic

I. INTRODUCTION of the solution in the form of intervals which contain all pos-
sible numerical errors.
Interval arithmetic (see, e.g., [10, 25, 26, 29]) realized Three main kinds of interval methods for solving the ini-
in floating-point computer arithmetic is a way to estimate tial value problem are known: the methods based on high-
two kinds of errors: representation errors of data (some order Taylor series (see, e.g., [1, 2, 4, 11, 28]), explicit and
real numbers cannot be represented exactly in floating-point implicit methods of Runge-Kutta type [6, 7, 15, 17, 23,
arithmetic) and rounding errors (the difference between 24, 29], and explicit and implicit multistep methods [12–
the calculated approximation of a number and its exact math- 14, 16, 17, 29]. So far, only in the first kind of interval
ematical value). Using approximate methods we introduce methods the step size correction has been applied. In inter-
the third kind of errors – the error of method (often called val methods based on Runge-Kutta methods and in interval
the truncation error, usually defined as the error that results multistep methods a constant step size has been used. Al-
from using an approximation in place of an exact mathemat- though the methods based on high-order Taylor series seem
ical procedure). Applying interval methods to approximate to be most universal, in [20–22, 24] we have shown that
the solution of the initial value problem in floating-point in- in some cases the interval methods of the second and third
terval arithmetic (see, e.g., [9]), we can obtain enclosures kind give better enclosures of solutions.
18 A. Marciniak, B. Szyszka

In this paper we propose an algorithm for finding the op- The local truncation error of step k + 1 for any Runge-
timal step size (taking into account the widths of inter- Kutta method of order p can be written in the form
vals obtained) for interval methods of Runge-Kutta type. !
m
We present an application of this algorithm taking into ac- X
count an interval method based on the classical Runge-Kutta rk+1 (h) = y(tk + h) − y (tk ) + h wi κik (h)
i=1
method of fourth order as an example. But the presented pro-
cedure can also be applied to other kinds of interval methods = ψ(tk , y(tk ))h p+1
+ O(hp+2
)
of Runge-Kutta type. (5)
p+1 p+2
The paper is divided into seven sections. In Sec. II. (p+1) h (p+2) h
= rk+1 (0) + rk+1 (θh) ,
and III. we shortly recall the well-known conventional (p + 1)! (p + 2)!
Runge-Kutta methods and the classical Runge-Kutta method
of fourth order, respectively. An interval version of fourth or- 0 < θ < 1,
der Runge-Kutta method [17, 29] is recalled in Sec. IV.. Sec.
V. is the main section of this paper, in which we describe where y(tk +h) and y(tk ) denote the exact solutions at tk +h
an algorithm for finding the optimal step size. In Sec. VI. and tk , respectively, and κik (h) is given by (3) for the ex-
(l)
we present some numerical examples which show the appli- act value y(tk ). From the conditions rk+1 (0) = 0 (for
cation of this algorithm. In the last section some conclusions l = 1, 2, ..., p) follow the equations for determining the co-
are given. efficients wi , ci and aij . Unfortunately, there are fewer equa-
tions than the number of unknowns, and usually we con-
sider some special cases. It is known [3, 8] that for each m
II. RUNGE-KUTTA METHODS there exists a method with maximum order pmax (m) = m
for m = 1, 2, 3, 4, pmax (m) = m − 1 for m = 5, 6, 7,
Let us consider the initial value problem pmax (m) = m − 2 for m = 8, 9, and pmax (m) ≤ m − 2 for
y 0 = f (t, y(t)), y(0) = y0 , (1) m ≥ 10.

where t ∈ [0, a], y ∈ R and f : [0, a] × R → R.

As it is well-known, the (explicit) m-stage Runge-Kutta III. THE CLASSICAL RUNGE-KUTTA METHOD
methods for solving the problem (1) are given by the for-
mula [3, 8] One of the most popular Runge-Kutta methods (simply
m
X called the Runge-Kutta method) is the four-stage fourth or-
yk+1 = yk + h wi κik , k = 0, 1, . . . , (2) der method of the form
i=1
h
where yk+1 = yk + (κ1k + 2κ2k + 2κ3k + κ4k ) , (6)
6
i−1
X
κik = f (tk + ci h, yk + h aij κjk ), where
j=1 (3)
κ1k = f (tk , yk ) ,
i = 1, 2, . . . , m,
 
h h
and κ2k = f tk + , yk + κ1k ,
2 2
i−1   (7)
h h
X
ci = aij , (4) κ3k = f tk + , yk + κ2k ,
j=1 2 2

h = tk+1 − tk is a step-size, and where the coefficients wi , κ4k = f (tk + h, yk + hκ3k ) .

ci and aij are some parameters. It is convenient to present
these coefficients in a form of a (triangular) array, called The Butcher table for this method is as follows:
the Butcher table: 0
0 1 1
2 2
c2 a21 1 1
c3 a31 a32 2 0 2
.. .. .. 1 0 0 1
. . . 
cm am1 am2 ··· amm 
1 2 2 1
w1 w2 ··· wm 6 6 6 6
 Interval Runge-Kutta Methods with Variable Step Sizes 19

The local truncation error is given by the following formula:

(3) (3) (3) (3)
rk+1 (h) = ψ(tk , y(tk ))h5 + O(h6 ). y (4) =ft3 + 3ft2 y f + 3fty2 f 2 + fy3 f 3
 
(2) (2)
The form of ψ(t, y) is rather complicated. Since this form + 3 fty + fy2 f y (2) + fy(1) y (3) ,
is very important from the point of view of the interval
(4) (4) (4) (4) (4)
method developed in the next section, below we present y (5) =ft4 + 4ft3 y f + 6ft2 y2 f 2 + 4fty3 f 3 + fy4 f 4
the adequate formulas.  
Denoting (to short the notations) (3) (3) (3)
+ 6 ft2 y + 2fty2 f + fy3 f 2 y (2)

(l) ∂lf 
(2) (2)

f = f (t, y), ftp yq = , + 4 fty + fy2 f y (3) + fy(1) y (4)
∂tp ∂y q
 2
(2)
where l = p + q, and + 3fy2 y (2) .
i−1
(l) (l) (l) (l)
X
y (l) = y (l) (t), κi = κi (0), λi = aij κj ,
j=1
IV. AN INTERVAL VERSION
OF RUNGE-KUTTA METHOD
we have
Let us denote:
• ∆t and ∆y – bounded sets in which the function
4
!
1 X (4)
ψ(t, y) = y (5) − 5 wi κi , (8) f (t, y), occurring in (1), is defined, i.e.,
120 i=1
∆t = {t ∈ R : 0 ≤ t ≤ a} ,
where 
∆y = y ∈ R : b ≤ y ≤ b ,
 
(1) (1)
κi =ci ft + fy(1) f , • F (T, Y ) – an interval extension of f (t, y), where
  an interval extension of the function
(2) (2) (2) (2) (1)
κi =c2i ft2 + 2fty f + fy2 f 2 + 2fy(1) λi ,
f : R × R ⊃ ∆t × ∆y → R
 
(3) (3) (3) (3) (3)
κi =c3i ft3 + 3ft2 y f + 3fty2 f 2 + fy3 f 3 we call a function

(2) (2)

(1) (2)
F : IR × IR ⊃ I∆t × I∆y → IR
+ 6ci fty + fy2 f λi + 3fy(1) λi ,
such that
 
(4) (4) (4) (4) (4) (4)
κi =c4i ft4 + 4ft3 y f + 6ft2 y2 f 2 + 4fty3 f 3 + fy 4 f 4 (t, y) ∈ (T, Y ) ⇒ f (t, y) ∈ F (T, Y ) ,
 
(3) (3)
+ 12c2i ft2 y + 2fty2 f + fy3 f 2 λi
(3) (1) and where IR denotes the space of real intervals,
• Ψ (T, Y ) – an interval extension of ψ (t, y) (see (8)),

(2) (2)

(2) (2)

(1)
2 and let us assume that:
+ 12ci fty + fy2 f λi + 12fy2 λi • the function F (T, Y ) is defined and continuous for all
T ⊂ ∆t and Y ⊂ ∆y 1 ,
(3)
+ 4fy(1) λi . • the function F (T, Y ) is monotonic with respect to in-
clusion, i.e.,
The derivatives of y with respect to t can be written
in the following forms: T1 ⊂ T2 ∧ Y1 ⊂ Y2 ⇒ F (T1 , Y1 ) ⊂ F (T2 , Y2 ) ,

• for each T ⊂ ∆t and for each Y ⊂ ∆y there exists

(1) a constant Λ > 0 such that
y (2) =ft + fy(1) f,
w (F (T, Y )) ≤ Λ (w (T ) + w (Y )) , (9)
(2) (2) (2)
y (3) =ft2 + 2fty f + fy2 f 2 + fy(1) y (2) ,
where w (A) denotes the width of the interval A,

1 The function F (T, Y ) is continuous at (T , Y ) if for every  > 0 there is a positive number δ = δ() such that d (F (T, Y ) , F (T , Y )) < 
0 0   0  0
whenever d(T, T0 ) < δ and d(Y, Y0 ) < δ. Here, d denotes the interval metric defined by d (X1 , X2 ) = max |X 1 − X 2 | , X 1 − X 2  , where
   
X1 = X 1 , X 1 and X2 = X 2 , X 2 are two intervals.
20 A. Marciniak, B. Szyszka

• the function Ψ (T, Y ) is defined for all T ⊂ ∆t and a procedure which in interval floating-point arithmetic cal-
Y ⊂ ∆y , culates the number η = tmax for any interval Runge-Kutta
• the function Ψ (T, Y ) is monotonic with respect to in- method (explicit or implicit).
clusion.
For t0 = 0 and y0 ∈ Y0 , where the interval Y0 is given, Finally, we divide the interval [0, η] into n parts
the interval version of Runge-Kutta method is of the form by the points tk = kh (k = 0, 1, ..., n), whereas the inter-
[17, 29] vals Tk , which appear in the method (10)–(11), are selected
in such a way that
h
Yk+1 =Yk + (K1k + 2K2k + 2K3k + K4k ) tk = kh ∈ Tk ⊂ [0, η] .
6
+ (Ψ (Tk , Yk ) + [−α, α]) h , 5 (10) For the method (10)–(11) we have the following

k = 0, 1, . . . , n − 1, Theorem 1 For the exact solution y(t) of the initial value

problem (1) we have y(tk ) ∈ Yk (k = 0, 1, ..., n), where Yk
where are obtained from (10)–(11).

K1k = F (Tk , Yk ) , The proof of this theorem can be found in [17, 29], where
  one can also find a theorem regarding estimations of w(Yk ).
h h
K2k = F Tk + , Yk + K1k ,
2 2
  Theorem 2 If Yk (k = 1, 2, ..., n) are obtained on the basis
h h
K3k = F Tk + , Yk + K2k , of the method (10)–(11), then
2 2 (11)
K4k = F (Tk + h, Yk + hK3k ) , w (Yk ) ≤ Qh4 + Rw (Y0 ) + S max w (Tl ) , (14)
l=0,1,...,k−1
 (6) 
r
 k+1 (θh)  where Q, R and S denote some nonnegative constants.

α = M h0 ,   ≤ M,
 720  The estimation (14) is true for any explicit interval
0 < θ < 1, 0 < h ≤ h0 , Runge-Kutta method of fourth order (not only for the method
(10)–(11)). In the case of the formulas (10)–(11) the con-
and where h0 denotes a given number (an initial value of step stants Q, R and S are as follows (see the proof of the theo-
size). rem in [17, 29]):
According to the theory of interval Runge-Kutta methods 
 S
[17, 29], the step size h of the method (10)–(11) is given by Q = w Ψ (∆t , ∆y ) + 2α ,
γΛ
(15)
η R =exp(γηΛ),
h= , (12)
n S =R − 1,
where where
4
X i−1
X j
η = min {η0 , η2 , η3 , η4 } , (13) γ= wi µij (h0 Λ) ,
i=1 j=0
and where for Y0 ⊂ ∆y and y0 ∈ Y0 the numbers ηi > 0 µi0 = 1, i = 1, 2, 3, 4,
(i = 2, 3, 4) are such that 1
µi1 = , i = 2, 3, 4,
2
Y0 + ηi ci F (∆t , ∆y ) ⊂ ∆y , i = 2, 3, 4, 1 1 1
µ32 = , µ42 = , µ43 = ,
4 2 4
and the number η0 > 0 fulfills the condition
and where w1 = w4 = 1/6, w2 = w3 = 1/3, and Λ is a con-
4
X stant occurring in (9).
Y0 + η0 wi F (∆t , ∆y )
i=1
V. CHANGING THE STEP SIZE
+ (Ψ (∆t , ∆y ) + [−α, α]) h40 ⊂ ∆y .
IN THE INTERVAL RUNGE-KUTTA METHOD

In the above relations we have (according to the Butcher In conventional Runge-Kutta methods the step size
table presented in Sec. III.) c2 = c3 = 1/2, c4 = 1, w1 = is changed to decrease errors of methods. In case of inter-
w4 = 1/6, w2 = w3 = 1/3. In [17] we have described val methods these errors are included in interval enclosures
 Interval Runge-Kutta Methods with Variable Step Sizes 21

of solutions. The only reason to change a given step size Doubling h we should always keep h ≤ h0 in mind.
for such methods is to decrease the widths of interval enclo- If we halve h, we should take into account the inequalities
sures. Since all calculations are performed in floating-point (16). Adding these inequalities by sides we have
arithmetic which produces rounding errors, for each problem   
considered there is a step size decreasing of which will not h
w (Y (h)) + w Y 2 × ≤
give intervals with smaller widths. A procedure for chang- 2
ing the step size for the method (10)–(11), described below, Q(R + 5) 4
is based on halving or doubling the size of a given step. h + R(R + 1)w (Y0 ) + Sw (T0 )
4
Let h ≤ h0 be a step calculated from (12), and let us de-

+S(R + 1)w T1/2 ≤
note by Y (h) an interval obtained from (10) for k = 0, i.e.,
at T1 (t1 = h ∈ T1 ). Let Y 2 × h2 be an interval obtained

Q(R + 5) 4

h + R(R + 1)w(Y0 ) + S(R + 2)w T1/2 ,
at T1 by an application of (10) two times. According to (14) 4
we have
from which it follows that
w (Y (h)) ≤ Qh4 + Rw (Y0 ) + Sw (T0 ) ,
w (Y (h)) + w Y 2 × h2



    4   
h h h
w Y 2× ≤Q


+ Rw Y 4 − R(R + 1)w (Y0 ) − S(R + 2)w T1/2
2 2 2 h ≥4 . (19)
+ Sw T1/2

Q(R + 5)
 4  4
h h Since in
a number of examples we have w (Y0 ) ≈
0 and
≤Q +R Q
2 2 w T1/2 ≈ 0, assuming w (Y0 ) = 0 and w T1/2 = 0,
! (16) from (19) we have
+ Rw (Y0 ) + Sw (T0 ) s
√ h



4 w (Y (h)) + w Y 2 × 2

h≥ 2· . (20)
+ Sw T1/2 Q(R + 5)
 4
h
≤ Q(R + 1) + R2 w (Y0 )
2
Applying the presented procedure, after a number
+ S(R + 1)w T1/2 , of steps we observe that halving and doubling the step size h
do not decrease the width of interval solution at some T1 ,
where we take advantage of w (T0 ) ≤ w T1/2 , t1/2 = h2 ∈


or h > h0 , or the inequalities (20) does not hold. In all
T1/2 . these case we accept the interval obtained with the last cor-
We should consider two cases: rect h. This procedure (presented in part A of block dia-
  
h gram in Fig. 1) gives us only the information  that the cor-
w (Y (h)) ≤ w Y 2 × (17)
rect step size h is within the interval h2 , 2h , but there can

2
and exist a step size h such that h2 ≤ h ≤ h or h ≤ h ≤ 2h,
for which the width of Y is smaller. With a machine ac-
curacy eps we are able to find such an h in the interval
  
h
w Y 2× < w (Y (h)) . (18) [h, 2h] for which the width of Y is w1 (part B in Fig. 1)
2
and such an h0 in the interval h2 , h with the width of Y
In case of (17) we see that by halving step size we do not equals w10 (part 0

C in Fig. 1). If w1 6= w1 , we substitute
obtain a better interval (i.e., with a smaller width) at T1 and h := h + h0 /2 and repeat part A with this new step size.
we can try to double h. In this case we substitute Otherwise (w1 = w10 ), the obtained h is the optimal step size
  (according to the machine accuracy and
because of h ≥ h0 ).
h
Y := Y (h), h := 2h, It means that the width of interval Y h is the smallest.
2 The above procedure may be applied for the next inte-
and calculate Y 2 × h2 and Y (h) for this new value of h.

gration steps, i.e., for k > 0. Since w (Yk+1 ) ≥ w (Yk ),
If (18) occurs, then the step size h should be halved. We sub- which follows immediately from (10), and taking into ac-
stitute count the expense of the described algorithm, for k > 0
we propose to use the constant step size obtained in the first
 
h h integration step, i.e., for k = 0. Although we have con-
Y (h) := Y , h := , sidered only the fourth order interval Runge-Kutta method,
2 2
the same algorithm can be applied to other interval meth-
h ods of Runge-Kutta type (explicit and implicit – see, e.g.,


and find Y 2 × 2 for the new h.
[6, 7, 15–17, 23, 24, 29]).
22 A. Marciniak, B. Szyszka

input: D t , D y, h, M, a, t max , h 0 , eps, T, Y

calculate Q and R from (15)

first_time := true
h OK := false 2 () ( )
calculate Y(h), Y h and Y 2 h
2 A
w11 := width( D y ) w1 := width(Y(h))
w1_less_w2 := false
w22 := w11
w2 _less_w1 := false (( )
w2 := width Y 2 h
2

no yes
w1 < w2

first_time := false no yes

w2_less_w1 := true first_time and (2h > h 0 )

first_time := false
no yes w1_less_w2 := true hOK := true
not w1_less_w2

h no yes no yes
h := w2 # w22 not w2_less_w1
2
hOK := true
h no yes
hOK := true h := hOK := true w1 # w11
2
h
h :=
no yes 2 h := 2h
not (20) hOK := true
no yes
w11 := w1 h := 2h h $ h0
w22 := w2 hOK := true
w11 := w1 h := h0
w22 := w2 hOK := true

no yes
not h OK

h 1/2 := 2
h
C
h N := h h 1/2 + h N
w 2N := w1 h :=
2
h := h calculate Y(h)
h 2:= 2h w1N := width(Y(h))
w2 := w1

no yes
w1N < wN2

h :=
h 2+ h B h 1/2 := h h N := h
2
calculate Y(h) w1N := w2N w2N := w1N
w1 := width(Y(h))

no yes no yes
w1 # w 2 | h N! h 1/2| # eps

h N := h
h 2:= h no yes
h := h w1 = w1N
w1 := w2
w2 := w1
h + hN
h := h := h
2
no yes
| h ! h 2 | # eps
output: h

Fig. 1. A block diagram for step size changing

In fact, we have presented our algorithm only for one- VI. NUMERICAL EXAMPLES
dimensional initial value problems, but it can be easily ex-
tended to systems of ordinary differential equations (interval In the examples presented below we have used our own
methods of Runge-Kutta type for a system of equations have implementation of floating-point interval arithmetic in Del-
been presented, among others, in [7], [17] and [23]). In this phi Pascal. This implementation has been written in the form
case the inequalities (17) and (18) should be replaced by of a unit called IntervalArithmetic32and64 (the current
version of this unit is presented in [18]). This unit takes ad-
N N    vantage of the Delphi Pascal floating-point Extended type.
X
X h
wi Y(i) (h) ≤ wi Y(i) 2 × All programs written in Delphi Pascal for the examples pre-
i=1 i=1
2 sented can be found in [19]. The first two examples concern
a commonly used test problem with a known exact solu-
and tion. For the third example the exact solution is unknown,
N    X N but we compare our results with those obtained by a method
X h

based on variable order Taylor series.
wi Y(i) 2 × < wi Y(i) (h) ,
i=1
2 i=1
VI. 1. Example 1
respectively, where N denotes the number of equations At first, let us consider the following simple initial value
in the system. problem:
 Interval Runge-Kutta Methods with Variable Step Sizes 23

y 0 = 0.5y, y(0) = 1, (21) Thus, let us consider the same initial value problem as in Ex-
ample 1 (with the same additional data), but now let us apply
with the exact solution y = exp(0.5t). Let us assume
the presented algorithm in each integration step. The changes
∆t = {t ∈ R : 0 ≤ t ≤ 10} , of step size in the interval [0, tmax ] are showed in Fig. 2, and
the interval enclosures obtained are presented in Tab. 6.
∆y = {y ∈ R : 0.9 ≤ y ≤ 149} ,

where x denotes the largest machine number less or equal

to x (similarly, by x we will denote the smallest machine
number greater or equal to x). We can find that M = 0.003,
and on the basis of (13) for h0 = 0.05 we have tmax = η ≈
1.9866. It appears that for h < h0 and for the first integra-
tion step, the algorithm presented in Fig. 1 gives the opti-
mal step size h = 7.66261590758908911E-4. This step size
is obtained regardless of whether the starting step size h is
slightly or much smaller/greater than the final step size h
(see Tab. 1–4).
Using optimal step size we obtain interval solutions pre-
sented in Tab. 5. Note that the integration step 2592 is
the last one for which t ∈ [0, tmax ]. Taking h = 0.0005
and h = 0.001 at the last integration for which t ∈ [0, tmax ]
we get

Y ([1.9864999999999997E+0000,
1.9865000000000001E+0000])
= [2.6999952128764781E+0000,
2.6999952128764800E+0000] Fig. 2. Changes of step size during calculations

and From Fig. 2 it follows that there is no regularity in step

size changes. We can only observe that the values of step
Y ([1.9859999999999998E+0000, size fluctuate more for greater t, i.e., for further integration
1.9860000000000001E+0000]) steps. It may be interesting that the mean value of step size
= [2.6993202984410782E+0000, on the whole integration interval [0, tmax ] equals 6.28E-4,
approximately, which differs insignificantly from the op-
2.6993202984410801E+0000]
timal step size for the first integration step (≈ 7.66E-4).
On the other hand, comparing the results in the last lines
with the widths approximately equal to 1.79E-15 and
in Tab. 5 and 6, we see that applying the algorithm for step
1.88E-15, respectively. Both these widths are greater than
size changing in each integration step brings in small profit
the widths of intervals obtained with the optimal step size.
(taking into account the widths of intervals obtained). Ob-
Obviously, in each case the exact solution is within the inter-
viously, the calculation time is much larger in this case and
val enclosure obtained.
to achieve tmax (in 3162 integration steps) we need 19021:00
All calculations have been carried out on a computer
min, approximately, i.e., more than 13 days (compare with
with Intel R
CoreTM i7–5500U CPU@ 2.40 GHz proces-
1:28 + 27:57 min in Example 1). 
sor. To find the optimal step size this computer needed about
1:28 min (1 minute and 28 seconds), while all other calcula-
In Examples 1 and 2 a lot of derivatives ftlp yq (l =
tions (to achieve tmax with this optimal step size, i.e., to exe-
p+q), which are needed to find Ψ (Tk , Yk ), are equal to zero
cute 2592 integration steps) lasted 27:57 min, approximately. (1)
 (only fy equals 0.5). In Example 3 we consider a problem
in which all these derivatives are different from zero.
VI. 2. Example 2 VI. 3. Example 3
In Example 1 we have found the optimal step size only For the initial value problem (the problem A5 from [5],
for the first integration step, and then we have used this step p. 23)
size as a constant one for further steps. Although the proce-
dure for step size changing is very expensive, it may be in- y−t
teresting to see how the step size changes in further steps. y0 = , y(0) = 4, (22)
y+t
24 A. Marciniak, B. Szyszka

Tab. 1. Finding the optimal step size for initial h = 0.0006

Step size considered

Step In interval
(part of algorithm)

1 h(A) [3.00000000000000000E-0004, 1.20000000000000000E-0003]

2 h(A) [6.00000000000000000E-0004, 2.40000000000000000E-0003]
initial h = 6.0000000000000000E-0004
in [h/2, 2h] = [3.00000000000000000E-0004, 1.20000000000000000E-0003],
checking h
in [≥ h, ≤ 2h] = [6.00000000000000000E-0004, 1.20000000000000000E-0003]
3 h(B) [9.00000000000000000E-0004, 1.20000000000000000E-0003]
4 h(B) [9.00000000000000000E-0004, 1.05000000000000000E-0003]
··· ··· ···
62 h(B) [9.54555771447376606E-0004, 9.54555771447376607E-0004]
checking h0
in [≥ h/2, ≤ h] = [3.00000000000000000E-0004, 9.54555771447376606E-0004]
63 h0 (C) [3.00000000000000000E-0004, 6.27277885723688303E-0004]
64 h0 (C) [4.63638942861844152E-0004, 6.27277885723688303E-0004]
··· ··· ···
122 h0 (C) [6.27277885723688303E-0004, 6.27277885723688303E-0004]


new h = h + h0 /2 = 7.90916828585532455E-0004
123 h(A) [3.95458414292766227E-0004, 1.58183365717106491E-0003]
initial h = 3.95458414292766227E-0004
in [h/2, 2h] = [1.97729207146383114E-0004, 7.90916828585532455E-0004]
checking h
in [≥ h, ≤ 2h] = [3.95458414292766227E-0004, 7.90916828585532455E-0004]
124 h(B) [5.93187621439149341E-0004, 7.90916828585532455E-0004]
125 h(B) [6.92052225012340898E-0004, 7.90916828585532455E-0004]
··· ··· ···
182 h(B) [7.66261590758908911E-0004, 7.66261590758908912E-0004]
checking h0
in [≥ h/2, ≤ h] = [1.97729207146383114E-0004, 7.66261590758908911E-0004]
183 h0 (C) [1.97729207146383114E-0004, 4.81995398952646012E-0004]
184 h0 (C) [3.39862303049514563E-0004, 4.81995398952646012E-0004]
··· ··· ···
242 h0 (C) [4.81995398952646012E-0004, 4.81995398952646012E-0004]
w1 = w10 ⇒ h := h = 7.66261590758908911E-0004

let us take presented in Tab. 7. Note that the intervals obtained are very
tiny.
∆t = {t ∈ R : 0 ≤ t ≤ 4} , Taking h = 1.46 − t1786 , from the last result we can find

∆y = y ∈ R : 4 ≤ y ≤ 6.3 , an interval enclosure at t = 1.46. We have
h0 = 0.01, M = 0.0537. [5.0849553259401612E+0000,
5.0849553259401642E+0000]
For these data and the method (10)–(11) we have found
with the width equals 2.87E-15, approximately. On the other
tmax ≈ 1.46. Using our algorithm we get h =
hand, the well-known VNODE-LP package [27] (based
8.17462272838888630E-0004 as the optimal step size for
on variable order Taylor series) produces
the first integration step. Using this value as a constant step
size for next integration steps we obtain interval enclosures 5.0849553259401[565, 699].2

2 Original output from VNODE-LP.

 Interval Runge-Kutta Methods with Variable Step Sizes 25

Tab. 2. Finding the optimal step size for initial h = 0.0008

Step size considered

Step In interval
(part of algorithm)

1 h(A) [4.00000000000000000E-0004, 1.60000000000000000E-0003]

2 h(A) [8.00000000000000000E-0004, 3.20000000000000000E-0003]
initial h = 8.0000000000000000E-0004
in [h/2, 2h] = [4.00000000000000000E-0004, 1.60000000000000000E-0003],
checking h
in [≥ h, ≤ 2h] = [8.00000000000000000E-0004, 1.60000000000000000E-0003]
3 h(B) [1.20000000000000000E-0003, 1.60000000000000000E-0003]
4 h(B) [1.20000000000000000E-0003, 1.40000000000000000E-0003]
··· ··· ···
61 h(B) [1.21444277883821840E-0003, 1.21444277883821841E-0003]
checking h0
in [≥ h/2, ≤ h] = [4.00000000000000000E-0004, 1.21444277883821840E-0003]
62 h0 (C) [4.00000000000000000E-0004, 8.07221389419109202E-0004]
63 h0 (C) [4.00000000000000000E-0004, 6.03610694709554601E-0004]
··· ··· ···
121 h0 (C) [6.03610694709554600E-0004, 6.03610694709554601E-0004]


new h = h + h0 /2 = 9.09026736773886502E-0004
122 h(A) [4.54513368386943251E-0004, 1.81805347354777300E-0003]
123 h(A) [2.27256684193471626E-0004, 9.09026736773886502E-0004]
initial h = 4.54513368386943251E-0004
in [h/2, 2h] = [2.27256684193471626E-0004, 9.09026736773886502E-0004]
checking h
in [≥ h, ≤ 2h] = [4.54513368386943251E-0004, 9.09026736773886502E-0004]
124 h(B) [6.81770052580414877E-0004, 9.09026736773886502E-0004]
125 h(B) [6.81770052580414877E-0004, 7.95398394677150689E-0004]
··· ··· ···
182 h(B) [7.66261590758908911E-0004, 7.66261590758908912E-0004]
checking h0
in [≥ h/2, ≤ h] = [2.27256684193471626E-0004, 7.66261590758908911E-0004]
183 h0 (C) [4.96759137476190268E-0004, 7.66261590758908911E-0004]
184 h0 (C) [6.31510364117549590E-0004, 7.66261590758908911E-0004]
··· ··· ···
241 h0 (C) [7.66261590758908910E-0004, 7.66261590758908911E-0004]
w1 = w10 ⇒ h := h = 7.66261590758908911E-0004

The width of this interval is equal to 1.34E-14. Thus, our Kutta methods it follows that one can use any h which fulfils
method gives better enclosure, but – on the other hand – the inequalities 0 < h ≤ h0 , we usually search for such an h
the VNODE-LP package gives a result very quickly.  which guarantees that the intervals obtained have the small-
est widths. So far, the only way to find such an h has been
VII. CONCLUSIONS based on a number of trials to apply the considered interval
methods for different values of h.
Until now, interval Runge-Kutta methods (explicit and In this paper we have presented an algorithm giving
implicit) have been considered only with a constant time the optimal step size for the next (for the first at the be-
step. This step size (h) has been chosen in such a way that ginning of calculations) integration step in interval Runge-
0 < h ≤ h0 , where h0 denotes some given step size for Kutta methods. The optimal step size means the step size
which the error of interval method is determined (see (10)– which gives a resulting interval with the smallest width. Al-
(11) for details). Although from the theory of interval Runge- though our algorithm can be applied in each integration step,
26 A. Marciniak, B. Szyszka

Tab. 3. Finding the optimal step size for initial h = 0.00001

Step size considered

Step In interval
(part of algorithm)

1 h(A) [5.00000000000000000E-0006, 2.00000000000000000E-0005]

2 h(A) [1.00000000000000000E-0005, 4.00000000000000000E-0005]
··· ··· ···
8 h(A) [6.40000000000000000E-0004, 2.56000000000000000E-0003]
initial h = 6.4000000000000000E-0004
in [h/2, 2h] = [3.20000000000000000E-0004, 1.28000000000000000E-0003],
checking h
in [≥ h, ≤ 2h] = [6.40000000000000000E-0004, 1.28000000000000000E-0003]
9 h(B) [9.60000000000000000E-0004, 1.28500000000000000E-0003]
10 h(B) [9.60000000000000000E-0004, 1.12000000000000000E-0003]
··· ··· ···
67 h(B) [1.05723384520852133E-0004, 1.05723384520852134E-0004]
checking h0
in [≥ h/2, ≤ h] = [4.00000000000000000E-0004, 1.21444277883821840E-0003]
68 h0 (C) [3.20000000000000000E-0004, 6.88616922604260667E-0004]
69 h0 (C) [5.04308461302130334E-0004, 6.88616922604260667E-0004]
··· ··· ···
127 h0 (C) [6.88616922604260667E-0004, 6.88616922604260667E-0004]


new h = h + h0 /2 = 8.72925383906391001E-0004
128 h(A) [4.36462691953195500E-0004, 1.74585076781278200E-0003]
initial h = 4.36462691953195500E-0004
in [h/2, 2h] = [2.18231345976597750E-0004, 8.72925383906391001E-0004]
checking h
in [≥ h, ≤ 2h] = [4.36462691953195500E-0004, 8.72925383906391001E-0004]
129 h(B) [6.54694037929793251E-0004, 8.72925383906391001E-0004]
130 h(B) [7.63809710918092126E-0004, 8.72925383906391001E-0004]
··· ··· ···
187 h(B) [7.66261590758908911E-0004, 7.66261590758908912E-0004]
checking h0
in [≥ h/2, ≤ h] = [2.18231345976597750E-0004, 7.66261590758908911E-0004]
188 h0 (C) [2.18231345976597750E-0004, 4.92246468367753331E-0004]
189 h0 (C) [3.55238907172175540E-0004, 4.92246468367753331E-0004]
··· ··· ···
246 h0 (C) [4.92246468367753330E-0004, 4.92246468367753331E-0004]
w1 = w10 ⇒ h := h = 7.66261590758908911E-0004

due to its complexity (see Example 2) we recommend using ries are commonly considered as most universal, sometimes
the algorithm only for the first integration step and then tak- other interval methods (not only of Runge-Kutta type) give
ing the obtained h as a constant step size for each succeeding better enclosures of the exact solutions (see examples pre-
steps. sented in [20–22, 24] and Example 3). This is the main rea-
In the examples presented in this paper we have ap- son to consider also possibilities of applying such methods
plied our algorithm to an interval version of classical Runge- for solving various initial value problems.
Kutta method of fourth order. However, there is no restric-
tion to apply the presented procedure to other kinds of in- Acknowledgment
terval methods of Runge-Kutta type (explicit and implicit – The paper was supported by the Poznan University
see, e.g., [16, 17, 23, 24, 29]). It should also be noted that of Technology (Poland) as part of Grant No. 09/91/DSPB/
although interval methods based on high-order Taylor se- 0628.
 Interval Runge-Kutta Methods with Variable Step Sizes 27

Tab. 4. Finding the optimal step size for initial h = 0.01

Step size considered

Step In interval
(part of algorithm)

1 h(A) [5.00000000000000000E-0003, 2.00000000000000000E-0002]

2 h(A) [2.50000000000000000E-0003, 1.00000000000000000E-0002]
··· ··· ···
5 h(A) [3.12500000000000000E-0004, 1.25000000000000000E-0003]
initial h = 6.25000000000000000E-0004
in [h/2, 2h] = [3.12500000000000000E-0004, 1.25000000000000000E-0003],
checking h
in [≥ h, ≤ 2h] = [6.25000000000000000E-0004, 1.25000000000000000E-0003]
6 h(B) [6.25000000000000000E-0004, 9.37500000000000000E-0004]
7 h(B) [6.25000000000000000E-0004, 7.81250000000000000E-0004]
··· ··· ···
65 h(B) [7.66261590758908911E-0004, 7.66261590758908912E-0004]
checking h0
in [≥ h/2, ≤ h] = [3.12500000000000000E-0004, 7.66261590758908911E-0004]
66 h0 (C) [5.39380795379454456E-0004, 7.66261590758908911E-0004]
67 h0 (C) [6.52821193069181683E-0004, 7.66261590758908911E-0004]
··· ··· ···
124 h0 (C) [7.66261590758908910E-0004, 7.66261590758908911E-0004]
w1 = w10 ⇒ h := h = 7.66261590758908911E-0004

Tab. 5. The interval enclosures to the solution of (21) obtained with the optimal step size for the first integration step

Step Y Width Exact solution

500 [1.2111440365720224E+0000, ≈ 1.78E-16 1.2111440365720225E+0000

1.2111440365720227E+0000]
1000 [1.4668698773229724E+0000, ≈ 3.96E-16 1.4668698773239726E+0000
1.4668698773239729E+0000]
1500 [1.7765907043480633E+0000, ≈ 6.59E-16 1.7765907043480636E+0000
1.7765907043480640E+0000]
2000 [2.1517072370004459E+0000, ≈ 9.98E-16 2.1517072370004464E+0000
2.1517072370004470E+0000]
2500 [2.6060273885419534E+0000, ≈ 1.45E-15 2.6060273885419541E+0000
2.6060273885419549E+0000]
2592 [2.6995228134287445E+0000, ≈ 1.57E-15 2.6995228134287453E+0000
2.6995228134287462E+0000]

References [3] J.C. Butcher, The Numerical Analysis of Ordinary Differen-

tial Equations: Runge-Kutta and General Linear Methods,
[1] M. Berz, G. Hoffstätter, Computation and Application of Tay- John Wiley & Sons, Chichester (1987).
lor Polynomials with Interval Remainder Bounds, Reliable [4] G.F. Corliss, R. Rihm, Validating an A Priori Enclosure
Computing 4(1), 83–97 (1998). Using High-Order Taylor Series, [In:] Scientific Comput-
[2] M. Berz, K. Makino, Performance of Taylor Model Methods ing, Computer Arithmetic, and Validated Numerics, 228–
for Validated Integration of ODEs, [In:] J. Dongarra, K. Mad- 238, Akademie Verlag (1996).
sen, J. Wasniewski (eds.) Applied Parallel Computing. State [5] W.H. Enright, J.D. Pryce, Two Fortran Packages for Assess-
of the Art in Scientific Computing, Lecture Notes in Computer ing Initial Value Methods, ACM Transactions on Mathemat-
Science 3732, 65–73 (2005). ical Software 13(1), 1–27 (1987).
28 A. Marciniak, B. Szyszka

Tab. 6. The interval enclosures to the solution of (21) obtained with a step size changing in each integration step

Step Y Width Exact solution

500 [1.2018073081170456E+0000, ≈ 1.08E-16 1.2018073081170458E+0000

1.2018073081170458E+0000]
1000 [1.4179877484171219E+0000, ≈ 2.17E-16 1.4179877484171221E+0000
1.4179877484171222E+0000]
1500 [1.6487894962346351E+0000, ≈ 3.26E-16 1.6487894962346353E+0000
1.6487894962346355E+0000]
2000 [1.8723210819363658E+0000, ≈ 4.49E-16 1.8723210819363661E+0000
1.8723210819363663E+0000]
2500 [2.1773586900207728E+0000, ≈ 6.39E-15 2.1773586900207732E+0000
2.1773586900207735E+0000]
3000 [2.5485293132332228E+0000, ≈ 9.12E-15 2.5485293132332234E+0000
2.5485293132332239E+0000]
3162 [2.6990947890180960E+0000, ≈ 1.01E-15 2.6990947890180967E+0000
2.6990947890180971E+0000]

Tab. 7. The interval enclosures to the solution of the problem (22) obtained with the optimal step size for the first integration step

Step T Y Width

500 [4.0873113641944430E-0001, [4.3717586653031176E+0000, ≈ 6.70E-16

4.0873113641944432E-0001] 4.3717586653031183E+0000]
1000 [8.1746227283888861E+0000, [4.6836807485176570E+0000, ≈ 1.46E-15
8.1746227283888865E+0000] 4.6836807485176585E+0000]
1500 [1.2261934092583329E+0000, [4.9498209108613228E+0000, ≈ 2.34E-15
1.2261934092583330E+0000] 4.9498209108613252E+0000]
1786 [1.4599876192902550E+0000, [5.0849484688085775E+0000, ≈ 2.87E-15
1.4599876192902552E+0000] 5.0849484688085804E+0000]

[6] K. Gajda, M. Jankowska, A. Marciniak, B. Szyszka, A Survey [11] K.R. Jackson, N.S. Nedialkov, Some Recent Advances in
of Interval Runge-Kutta and Multistep Methods for Solving Validated Methods for IVPs for ODEs, Applied Numerical
the Initial Value Problem, [In:] R. Wyrzykowski, J. Don- Mathematics 42, 269–284 (2002).
garra, K. Karczewski, J. Wasniewski (eds.) Parallel Process- [12] M. Jankowska, A. Marciniak, Implicit Interval Methods for
ing and Applied Mathematics, Lecture Notes in Computer Solving the Initial Value Problem, Computational Methods
Science 4967, 1361–1371. Springer-Verlag, Berlin (2007). in Science and Technology 8(1), 17–30 (2002).
[7] K. Gajda, A. Marciniak, B. Szyszka, Three-and Four-Stage [13] M. Jankowska, A. Marciniak, On Explicit Interval Methods
Implicit Interval Methods of Runge-Kutta Type, Computa- of Adams-Bashforth Type, Computational Methods in Sci-
tional Methods in Science and Technology 6, 41–59 (2000). ence and Technology 8(2), 46–57 (2002).
[8] E. Hairer, S.P. Norsett, G. Wanner, Solving Ordinary Dif- [14] M. Jankowska, A. Marciniak, On Two Families of Implicit
ferential Equations I – Nonstiff Problems, Springer–Verlag, Interval Methods of Adams-Moulton Type, Computational
Berlin (1987). Methods in Science and Technology 12(2), 109–113 (2006).
[9] R. Hammer, M. Hocks, U. Kulisch, D. Ratz, Numerical Tool- [15] S.A. Kalmykov, J.I. Shokin, Z.H. Juldashev, Solving Ordi-
box for Verified Computing I. Basic Numerical Problems, nary Differential Equations by Interval Methods [in Rus-
Theory, Algorithms, and Pascal-XSC Programs, Springer- sian], Doklady AN SSSR 230(6) (1976).
Verlag, Berlin (1993). [16] A. Marciniak, Implicit Interval Methods for Solving the Ini-
[10] E.R. Hansen, Topics in Interval Analysis, Oxford University tial Value Problem, Numerical Algorithms 37, 241–251
Press, London (1969). (2004).
 Interval Runge-Kutta Methods with Variable Step Sizes 29

[17] A. Marciniak, Selected Interval Methods for Solving the Ini- [23] A. Marciniak, B. Szyszka, One-and Two-Stage Implicit Inter-
tial Value Problem, Publishing House of Poznan University val Methods of Runge-Kutta Type, Computational Methods
of Technology, Poznan (2009). http://www.cs.put.poznan.pl/ in Science and Technology 5, 53–65 (1999).
amarciniak/IMforIVP-book/IMforIVP.pdf [24] A. Marciniak, B. Szyszka, T. Hoffmann, An Interval Version
[18] A. Marciniak, Interval Arithmetic Unit (2016). http://www.cs of Kuntzmann-Butcher Method for Solving the Initial Value
.put.poznan.pl/amarciniak/IAUnits/IntervalArithmetic32and6 Problem (in review, available from the authors).
4.pas [25] R.E. Moore, Interval Analysis, Prentice-Hall, Englewood
Cliffs (1966).
[19] A. Marciniak, Delphi Pascal Programs for Step Size Control [26] R.E. Moore, Methods and Applications of Interval Analy-
in Interval Runge-Kutta Methods (2017). http://www.cs.put.p sis, SIAM Society for Industrial & Applied Mathematics,
oznan.pl/amarciniak/VSSIRKM-Examples Philadelphia (1979).
[20] A. Marciniak, M.A. Jankowska, Interval Versions for Special [27] N.S. Nedialkov, VNODE-LP – a Validated Solver for Initial
Kinds of Explicit Linear Multistep Methods (in review, avail- Value Problems in Ordinary Differential Equations, Tech.
able from the authors). Rep. CAS 06-06-NN, Department of Computing and Soft-
ware, McMaster University, Hamilton (2006).
[21] A. Marciniak, M.A. Jankowska, Interval Versions of Milne’s [28] N.S. Nedialkov, K.R. Jackson, G.F. Corliss, Validated So-
Multistep Methods, Numerical Algorithms 79(1), 87–105 lutions of Initial Value Problems for Ordinary Differential
(2018). Equations, Applied Mathematics and Computation 105(1),
[22] A. Marciniak, M.A. Jankowska, T. Hoffmann, On Interval 21–68 (1999).
Predictor-Corrector Methods, Numerical Algorithms 77(3), [29] Y.I. Shokin, Interval Analysis [in Russian], Nauka, Novosi-
777–808 (2017). birsk (1981).

Andrzej Marciniak was born in Poznań (Poland) in 1953. He received the M.Sc. degree in mathe-
matics in 1977, the M.Sc. degree in astronomy in 1979 and Ph.D. degree in mathematics in 1981, all
from the Adam Mickiewicz University in Poznań. In 1993 he received the Dr.Habil. degree in physics
from the Nicolaus Copernicus Univeristy in Toruń (Poland) and in 2010 he received the Professor Ti-
tle from the President of Poland. From 1977 to 1987 and from 2000 to 2011 he held a research
position at the Faculty of Mathematics and Computer Science of the Adam Mickiewicz Univeristy,
and since 1987 he has been an an assistant professor in Institute of Mathematics and then a professor
of computer science at the Faculty of Computing Science of the Poznań University of Technology.
From 2005 to 2008 he held the office of the President of Polish Information Processing Society.
His research interests include computer programming and numerical methods, especially for solving
ordinary and partial defferential equations with applications to dynamical problems. In these fields
he wrote three monographs, more than 20 textbooks and a number of scientific articles.

Barbara Szyszka received the M.Sc. degree in Applied Mathematics in 1994 and Ph.D. degree
in Computer Science in 2003, both from the Poznan University of Technology. From 1994 she works
in the Poznan University of Technology in the Institute of Mathematics. Her fields of research interest
include programming and numerical methods, in particular difference methods for solving ordinary
and partial differential equations.

CMST 25(1) 17–29 (2019) DOI:10.12921/cmst.2019.0000006

View publication stats