Explicit Runge–Kutta Methods Combined with Advanced Versions of the Richardson Extrapolation

Zahari Zlatev; Ivan Dimov; István Faragó; Krassimir Georgiev; Ágnes Havasi

doi:10.1515/cmam-2019-0016

Article Open Access

Explicit Runge–Kutta Methods Combined with Advanced Versions of the Richardson Extrapolation

Zahari Zlatev , Ivan Dimov , István Faragó , Krassimir Georgiev and Ágnes Havasi

Published/Copyright: July 16, 2019

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Computational Methods in Applied Mathematics Volume 20 Issue 4

Abstract

Richardson Extrapolation is a very general numerical procedure, which can be applied in the solution of many mathematical problems in an attempt to increase the accuracy of the results. It is assumed that this approach is used to handle non-linear systems of ordinary differential equations (ODEs) which arise often in the mathematical description of scientific and engineering models either directly or after the discretization of the spatial derivatives of partial differential equations (PDEs). The major topic is the analysis of eight advanced implementations of the Richardson Extrapolation. Two important properties are analyzed: (a) the possibility to achieve more accurate results and (b) the possibility to improve the stability properties of eight advanced versions of the Richardson Extrapolation. A two-parameter family of test-examples was constructed and used to check both the accuracy and the absolute stability of the different versions of the Richardson Extrapolation when these versions are applied together with several Explicit Runge–Kutta Methods (ERKMs).

Keywords: Several Times Repeated Richardson Extrapolation; Order of Accuracy; Absolute Stability Properties; Systems of ODEs; Runge–Kutta Methods; Numerical Examples

MSC 2010: 65L05; 65L06; 65L04

1 Introduction

It is well known that the Richardson Extrapolation (introduced in [23], see also [7, 32, 33, 34, 36, 37]) is a very general approach, which can successfully be used during the approximate solution of different classes of mathematical problems in the efforts (a) to increase the accuracy of the selected numerical methods and/or (b) to control the stepsize when the treated mathematical problems are time-dependent. On the other hand, it is also well known that the application of the Richardson Extrapolation in connection with time-depending problems is sometimes causing problems related to the absolute stability during the computational process ([9, 34, 36]). Therefore, the absolute stability properties of the new numerical methods, which are combinations of the applied numerical algorithms for treating different mathematical problems with several versions of the Richardson Extrapolation, must be studied very carefully. In this paper we shall

introduce a unified approach for the application of the Richardson Extrapolation,
derive general formulae representing different highly accurate versions of the Richardson Extrapolation,
show that the application of more advanced versions of the Richardson Extrapolation leads to improved absolute stability properties when Explicit Runge–Kutta Methods, ERKMs, are used,
present a two-parameter family of meaningful examples which can and will be used in the verification of both the accuracy of the calculated results and the absolute stability,
discuss the obtained numerical results as well as the possibilities for achieving further improvements,
finish the presentation with several general conclusions.

We shall first introduce some convenient notations, which will after that be consistently used in this paper.

2 Notation

Assume that the mathematical problem solved is a non-linear initial value problem for systems of ordinary differential equations (ODEs) defined in the following manner:

(2.1) d ⁢ y d ⁢ t = f ⁢ ( t , y ) , t ∈ [ a , b ] , a < b , y ∈ ℝ s , s ≥ 1 , f ∈ D ⊂ ℝ × ℝ s , y ⁢ ( a ) = η .

The exact solution of (2.1) is y ⁢ ( t ) and we shall first consider the case where approximated values of this solution have to be calculated by applying some arbitrary one-step method ([5, 6, 13, 14, 18, 20, 22]) for solving systems of ODEs either directly or in combination with some version of the Richardson Extrapolation. We shall assume that approximations of the exact solution of (2.1) are calculated on a set of equidistant grid-points:

(2.2) t 0 = a , t n = t n - 1 + h ( n = 1 , 2 , … , N ) , t N = b , h = b - a N ,

but this is done mainly in order to facilitate the presentation of the results. Most of the results are also valid when non-equidistant grids are used. Under the above assumptions, the following notations will be consistently used in the remaining part of this paper:

The value of the numerical solution of (2.1) calculated by using directly the selected one-step method at the grid-point t n ( n = 1 , 2 , … , N ) of (2.2) will always be denoted by y n . It must be mentioned here that the fact that a one-step method is used in the computations means that the approximation y n ≈ y ⁢ ( t n ) can be calculated, for any n = 1 , 2 , … , N , by using only the previous approximation y n - 1 ≈ y ⁢ ( t n - 1 ) .
The value of the numerical solution of (2.1), which is found by the Classical Richardson Extrapolation (introduced in [23], see also [36, 37, 40]) and which is calculated at the grid-point t n ( n = 1 , 2 , … , N ) of (2.2), will always be denoted by

y n [ 0 ] ≈ y ⁢ ( t n ) .
The numerical solution found by using the q-times Repeated Richardson Extrapolation ( q = 1 , 2 , … , 7 ) and evaluated at the grid-point t n ( n = 1 , 2 , … , N ) of (2.2) will always be denoted by

y n [ q ] ≈ y ⁢ ( t n ) , q = 1 , 2 , … , 7 .

The q-times repeated Richardson Extrapolation will be defined and discussed in the next section. Moreover, it will be shown there how the terms in the right-hand side of the following relationship can be obtained (also for q = 0 ):

(2.3) y n [ q ] = P n [ q ] S [ q ] , n = 1 , 2 , … , N ⁢ and ⁢ q = 1 , 2 , … , 7 .

3 Description of Advanced Versions of the Richardson Extrapolation

Definition 3.1.

Advanced versions of the Richardson Extrapolation can be introduced in the following manner. Let q ∈ { 0 , 1 , 2 , … } be some fixed integer, and assume that the approximation y n - 1 [ q ] ≈ y ⁢ ( t n - 1 ) is available. Consider the application of an arbitrary one-step method in the numerical solution of (2.1) on the grid-points defined by (2.2) and calculate q + 2 approximations z n [ r ] , r = 0 , 1 , 2 , … , q + 1 , under the following three conditions:

the starting approximation used in the calculation of all approximations z n [ r ] is y n - 1 [ q ] ,
the number of steps needed to calculate z n [ r ] is 2 r ,
the stepsize used in the calculation of z n [ r ] is h 2 r , where h is defined by (2.2).

Use the approximations z n [ r ] , r = 0 , 1 , 2 , … , q + 1 , to calculate y n [ q ] ≈ y ⁢ ( t n ) and apply after that the same approach in the computation of the next approximation y n + 1 [ q ] ≈ y ⁢ ( t n + 1 ) .

Two comments are needed:

Definition 3.1 is applicable for q > 7 , too, but only versions of the Richardson Extrapolation for 0 ≤ q ≤ 7 will be studied here,
it is necessary to explain how y n [ q ] ≈ y ⁢ ( t n ) can be calculated by using z n [ r ] , r = 0 , 1 , 2 , … , q + 1 .

The answer to this important question will be postponed to the next section, where the order of accuracy of y n [ q ] will be studied.

Example 1.

The well-known Classical Richardson Extrapolation, [10], is obtained by setting q = 0 . It is necessary to calculate z n [ 0 ] by carrying out one step with a stepsize h and then z n [ 1 ] by performing two steps with a stepsize h 2 . The starting approximation in the calculation of z n [ 0 ] and z n [ 1 ] is y n - 1 [ 0 ] ≈ y ⁢ ( t n - 1 ) . The approximations z n [ 0 ] and z n [ 1 ] are used to calculate y n [ 0 ] ≈ y ⁢ ( t n ) .

Example 2.

The Repeated Richardson Extrapolation, [32], see also [7] and [34], is obtained by setting q = 1 . It is necessary to calculate first z n [ 0 ] and z n [ 1 ] as in Example 1 and z n [ 2 ] by performing four steps with a stepsize h 4 . The starting approximation in the calculation of z n [ 0 ] , z n [ 1 ] and z n [ 2 ] is y n - 1 [ 1 ] ≈ y ⁢ ( t n - 1 ) . The approximations z n [ 0 ] , z n [ 1 ] and z n [ 2 ] are used to calculate y n [ 1 ] ≈ y ⁢ ( t n ) .

Example 3.

The Two-times Repeated Richardson Extrapolation, [33], is obtained by setting q = 2 . It is necessary to calculate the three approximations z n [ 0 ] , z n [ 1 ] and z n [ 2 ] as in Example 2 and z n [ 3 ] by performing eight steps with a stepsize h 8 . The starting approximation in the calculation of z n [ 0 ] , z n [ 1 ] , z n [ 2 ] and z n [ 3 ] is y n - 1 [ 2 ] ≈ y ⁢ ( t n - 1 ) . The approximations z n [ 0 ] , z n [ 1 ] , z n [ 2 ] and z n [ 3 ] are used to calculate y n [ 2 ] ≈ y ⁢ ( t n ) .

Two important conclusions can be drawn from the above examples: (a) the process can easily be continued in an obvious way to calculate successively the Three-, Four-, Five-, Six- and Seven-times Repeated Richardson Extrapolation and (b) the computational work is increasing very quickly, but it will be proved in the next section that the accuracy is also increased very considerably.

4 Accuracy of the Different Versions of the Richardson Extrapolation

Theorem 4.1.

Consider an arbitrary one-step method for solving the system of ODEs (2.1). Assume that

its order of accuracy is p ≥ 1 ,
it is used together with some version of the Richardson Extrapolation with 0 ≤ q ≤ 7 .

Then the order of accuracy of the approximation calculated by the selected version of the Richardson Extrapolation is at least p + q + 1 when the right-hand-side f ⁢ ( t , y ) of (2.1) is at least p + q + 1 times continuously differentiable.

Proof.

The following q + 2 equalities hold when the conditions of Theorem 4.1 are satisfied:

(4.1) y ⁢ ( t n ) - z n [ i ] = ∑ j = 1 q + 1 h p + j - 1 2 i ⁢ ( p + j - 1 ) ⁢ K j + 𝒪 ⁢ ( h p + q + 1 ) , i = 0 , 1 , … , q + 1 ,

where K j , j = 1 , 2 , … , q + 1 , are constants which do not depend on the stepsize h, and the theorem will be proved if all terms in the right-hand side of (4.1) that contain these constants are eliminated. It will be explained below how the elimination of the constants K j , j = 1 , 2 , … , q + 1 , can be performed for q = 0 , 1 , 2 , … , 7 .

Classical Richardson Extrapolation. The q + 2 relationships (4.1) are reduced to the following two relationships in the special case where q = 0 :

(4.2) y ⁢ ( t n ) - z n [ 0 ] = h p ⁢ K 1 + 𝒪 ⁢ ( h p + 1 ) ,

(4.3) y ⁢ ( t n ) - z n [ 1 ] = h p 2 p ⁢ K 1 + 𝒪 ⁢ ( h p + 1 ) .

Multiply (4.3) by 2 p and subtract (4.2) from the so-modified relationship (4.3). The result is

y ⁢ ( t n ) = 2 p ⁢ z n [ 1 ] - z n [ 0 ] 2 p - 1 + 𝒪 ⁢ ( h p + 1 ) .

Denote

(4.4) y n [ 0 ] = 2 p ⁢ z n [ 1 ] - z n [ 0 ] 2 p - 1 = : P n [ 0 ] S [ 0 ] .

Then

(4.5) y ⁢ ( t n ) - y n [ 0 ] = 𝒪 ⁢ ( h p + 1 ) .

Equality (4.4) shows how the approximation y n [ 0 ] can be calculated by using the already computed approximations z n [ 0 ] and z n [ 1 ] , while equality (4.5) is telling us that the order of accuracy of the Classical Richardson Extrapolation is p + 1 .

Repeated Richardson Extrapolation. The q + 2 relationships (4.1) are reduced to the following three relationships in the special case where q = 1 :

(4.6) y ⁢ ( t n ) - z n [ 0 ] = h p ⁢ K 1 + h p + 1 ⁢ K 2 + 𝒪 ⁢ ( h p + 2 ) ,

(4.7) y ⁢ ( t n ) - z n [ 1 ] = h p 2 p ⁢ K 1 + h p + 1 2 p + 1 ⁢ K 2 + 𝒪 ⁢ ( h p + 2 ) ,

(4.8) y ⁢ ( t n ) - z n [ 2 ] = h p 4 p ⁢ K 1 + h p + 1 4 p + 1 ⁢ K 2 + 𝒪 ⁢ ( h p + 2 ) .

Multiply (4.7) with 2 p and subtract (4.6) from the so-modified relationship (4.7). Multiply after that (4.8) with 4 p and subtract (4.6) from the so-modified relationship (4.8). In this way the first constant K 1 will be eliminated and two equalities which depend only on K 2 will be obtained. A quite similar procedure can be applied to eliminate also the second constant K 2 . The final result is

(4.9) y ⁢ ( t n ) = 2 2 ⁢ p + 1 ⁢ z n [ 2 ] - 3 ⋅ 2 p ⁢ z n [ 1 ] + z n [ 0 ] 2 2 ⁢ p + 1 - 3 ⋅ 2 p + 1 + 𝒪 ⁢ ( h p + 2 ) .

Denote

(4.10) y n [ 1 ] = 2 2 ⁢ p + 1 ⁢ z n [ 2 ] - 3 ⋅ 2 p ⁢ z n [ 1 ] + z n [ 0 ] 2 2 ⁢ p + 1 - 3 ⋅ 2 p + 1 = : P n [ 1 ] S [ 1 ] .

Then

(4.11) y ⁢ ( t n ) - y n [ 1 ] = 𝒪 ⁢ ( h p + 2 ) .

Equality (4.10) shows how the approximation y n [ 1 ] can be calculated by using the already computed approximations z n [ 0 ] , z n [ 1 ] and z n [ 2 ] , while equality (4.11) is telling us that the order of accuracy of the Repeated Richardson Extrapolation is p + 2 .

Two-times Repeated Richardson Extrapolation. The q + 2 relationships (4.1) are reduced to the following four relationships in the special case where q = 2 :

(4.12) y ⁢ ( t n ) - z n [ 0 ] = h p ⁢ K 1 + h p + 1 ⁢ K 2 + h p + 2 ⁢ K 3 + 𝒪 ⁢ ( h p + 3 ) ,

(4.13) y ⁢ ( t n ) - z n [ 1 ] = h p 2 p ⁢ K 1 + h p + 1 2 p + 1 ⁢ K 2 + h p + 2 2 p + 2 ⁢ K 3 + 𝒪 ⁢ ( h p + 3 ) ,

(4.14) y ⁢ ( t n ) - z n [ 2 ] = h p 4 p ⁢ K 1 + h p + 1 4 p + 1 ⁢ K 2 + h p + 2 4 p + 2 ⁢ K 3 + 𝒪 ⁢ ( h p + 3 ) ,

(4.15) y ⁢ ( t n ) - z n [ 3 ] = h p 8 p ⁢ K 1 + h p + 1 8 p + 1 ⁢ K 2 + h p + 2 8 p + 2 ⁢ K 3 + 𝒪 ⁢ ( h p + 3 ) .

Multiply (4.13) with 2 p and subtract (4.12) from the so-modified relationship (4.13). Multiply after that (4.14) with 4 p and subtract (4.12) from the so-modified relationship (4.14). Finally, multiply (4.15) with 8 p and subtract (4.12) from the so-modified relationship (4.15). In this way the first constant K 1 will be eliminated and three equalities which depend only on K 2 and K 3 will be obtained. A quite similar procedure can be applied twice in order to eliminate first the second constant K 2 and then the third constant K 3 . The final result is

y ⁢ ( t n ) = 2 3 ⁢ p + 3 ⁢ z n [ 3 ] - 7 ⋅ 2 2 ⁢ p + 1 ⁢ z n [ 2 ] + 7 ⋅ 2 p ⁢ z n [ 1 ] - z n [ 0 ] 2 3 ⁢ p + 3 - 7 ⋅ 2 2 ⁢ p + 1 + 7 ⋅ 2 p - 1 + 𝒪 ⁢ ( h p + 3 ) .

Denote

(4.16) y n [ 2 ] = 2 3 ⁢ p + 3 ⁢ z n [ 3 ] - 7 ⋅ 2 2 ⁢ p + 1 ⁢ z n [ 2 ] + 7 ⋅ 2 p ⁢ z n [ 1 ] - z n [ 0 ] 2 3 ⁢ p + 3 - 7 ⋅ 2 2 ⁢ p + 1 + 7 ⋅ 2 p - 1 = : P n [ 2 ] S [ 2 ] .

Then

(4.17) y ⁢ ( t n ) - y n [ 2 ] = 𝒪 ⁢ ( h p + 3 ) .

As in the previous two cases, equality (4.16) shows how y n [ 2 ] can be calculated by using the already computed approximations z n [ 0 ] , z n [ 1 ] , z n [ 2 ] and z n [ 3 ] , while equality (4.17) is telling us that the order of accuracy of the Two-times Repeated Richardson Extrapolation is p + 3 .

Five More Advanced Versions of the Repeated Richardson Extrapolation. The number of constants is increased if these versions are used, being four when q = 3 , five when q = 4 , six when q = 5 , seven when q = 6 and eight when q = 7 , but the main idea remains the same: all these constants must be eliminated and formulae for performing Three-, Four-, Five-, Six- and Seven-times Repeated Richardson Extrapolation will be derived at the end of the elimination process for the appropriate value of q. It should be noted here that the main principle is very clear, quite straight-forward and extremely simple, but the computational process becomes more and more complicated when the value of q is increased. The constants in (4.1) were successively eliminated for q = 3 , 4 , 5 , 6 , 7 and it has been shown that (2.3), which was satisfied for q = 0 , 1 , 2 , is also valid for q = 3 , 4 , 5 , 6 , 7 . Moreover, it is obvious that if we know S [ q ] in some formula (2.3), then we can easily construct the numerator P n [ q ] by multiplying successively the terms in the denominator by z n [ q + 1 ] , z n [ q ] , … , z n [ 0 ] , respectively. This is why only the denominators S [ q ] of formulae (2.3) for q = 3 , 4 , 5 , 6 , 7 that are obtained after the elimination of the constants in (4.1) are listed below:

(4.18) S [ 3 ] = 2 4 ⁢ p + 6 - 15 ⋅ 2 3 ⁢ p + 3 + 35 ⋅ 2 2 ⁢ p + 1 + 1 ,

(4.19) S [ 4 ] = 2 5 ⁢ p + 10 - 31 ⋅ 2 4 ⁢ p + 6 + 155 ⋅ 2 3 ⁢ p + 3 - 155 ⋅ 2 2 ⁢ p + 1 + 31 ⋅ 2 p - 1 ,

(4.20) S [ 5 ] = 2 6 ⁢ p + 15 - 63 ⋅ 2 5 ⁢ p + 10 + 651 ⋅ 2 4 ⁢ p + 6 - 1395 ⋅ 2 3 ⁢ p + 3 + 651 ⋅ 2 2 ⁢ p + 1 - 63 ⋅ 2 p + 1 ,

S [ 6 ] = 2 7 ⁢ p + 21 - 127 ⋅ 2 6 ⁢ p + 15 + 2667 ⋅ 2 5 ⁢ p + 10 - 11811 ⋅ 2 4 ⁢ p + 6

(4.21) + 11811 ⋅ 2 3 ⁢ p + 3 - 2667 ⋅ 2 2 ⁢ p + 1 + 127 ⋅ 2 p - 1 ,

S [ 7 ] = 2 8 ⁢ p + 28 - 255 ⋅ 2 7 ⁢ p + 21 + 10795 ⋅ 2 6 ⁢ p + 15 - 97155 ⋅ 2 5 ⁢ p + 10

(4.22) + 200787 ⋅ 2 4 ⁢ p + 6 - 97155 ⋅ 2 3 ⁢ p + 3 + 10795 ⋅ 2 2 ⁢ p + 1 - 255 ⋅ 2 p + 1 .

This completes the proof of Theorem 4.1. ∎

Formulae (4.18)–(4.22) were derived by eliminating directly constants in (4.1). The major idea is, as mentioned above, straight-forward, but the computational process is becoming too long and very difficult when q is large. However, there is also a much easier way for obtaining the above results.

Theorem 4.2.

Consider any version of the Richardson Extrapolation with 0 ≤ q ≤ 7 . Then its denominator can be calculated by using the following generic formula:

(4.23) S [ q ] = ∏ j = 1 q + 1 ( 2 p + j - 1 - 1 ) .

Moreover, the following relationships follow immediately from (4.23):

(4.24) S [ q ] = ( 2 p + q - 1 ) ⁢ S [ q - 1 ] , q = 1 , 2 , … , 7 .

Proof.

It can be seen from (4.4) that S [ 0 ] = 2 p - 1 . The application of formula (4.24) with q = 1 leads to S [ 1 ] = ( 2 p + 1 - 1 ) ⁢ S [ 0 ] . Replace S [ 0 ] with 2 p - 1 to obtain S [ 1 ] = ( 2 p + 1 - 1 ) ⁢ ( 2 p - 1 ) = 2 p + 1 - 3 ⋅ 2 p + 1 . The last expression is precisely the denominator of the Repeated Richardson Extrapolation, see (4.9), which means that (4.24) is true for q = 1 . Continuing in a quite similar manner, one can successively prove that (4.24) holds also for q = 2 , 3 , … , 7 . ∎

Corollary 4.3.

If we know S [ q ] for some q = 0 , 1 , … , 6 , then S [ q + 1 ] (and after that also P n [ q ] and the final formula y n [ q + 1 ] = P n [ q + 1 ] S [ q + 1 ] ) can be calculated by using (4.24).

Conjecture 1.

The assertion of Theorem 4.1 is valid not only for 0 ≤ q ≤ 7 , but also for any q > 7 .

Conjecture 2.

Theorem 4.2 and Corollary 4.3 are true for any q > 7 .

The results presented at the end of this paper and obtained by performing runs with the Eight-times Repeated Richardson Extrapolation indicate that the two conjectures may be correct. The denominator of the Eight-times Repeated Richardson Extrapolation used in these experiments was calculated by using the generic formula S [ 8 ] = ( 2 p + 8 - 1 ) ⁢ S [ 7 ] , and it is given by the expression

S [ 8 ] = 2 9 ⁢ p + 36 - 511 ⋅ 2 8 ⁢ p + 28 + 43435 ⋅ 2 7 ⁢ p + 21 - 788035 ⋅ 2 6 ⁢ p + 15

+ 3309747 ⋅ 2 5 ⁢ p + 10 - 3309747 ⋅ 2 4 ⁢ p + 6 + 788035 ⋅ 2 3 ⁢ p + 3

- 43435 ⋅ 2 2 ⁢ p + 1 + 511 ⋅ 2 p - 1 .

5 Stability Properties of the Advanced Versions of the Richardson Extrapolation

Germund Dahlquist proposed in [9] to use a simple linear and scalar equation:

(5.1) d ⁢ y d ⁢ t = λ ⁢ y , t ∈ [ 0 , ∞ ) , y ∈ ℂ , λ = α + β ⁢ i ∈ ℂ - , α ≤ 0 , y ⁢ ( 0 ) = η ∈ ℂ ,

solved on the equidistant grid

t 0 = 0 , t n = t n - 1 + h = t 0 + n h ( n = 1 , 2 , … ) , h > 0 ,

in order to study the absolute stability properties of numerical methods for solving ODEs. The exact solution of (5.1) is y ⁢ ( t ) = η ⁢ e λ ⁢ t and it is bounded when the condition α ≤ 0 is satisfied. This implies that it is necessary to require that the selected numerical method produces arbitrarily long sequences of approximations { y 1 , y 2 , … , y N , … } which are also bounded when α ≤ 0 . The application of an arbitrary one-step numerical method in the treatment of (5.1) leads to the following recursive relations (see, for example, [22]):

(5.2) y n = R ⁢ ( ν ) , y n - 1 = [ R ⁢ ( ν ) ] 2 , y n - 2 = … = [ R ⁢ ( ν ) ] n ⁢ y 0 , ν = λ ⁢ h , n = 1 , 2 , … .

The function R ⁢ ( ν ) from (5.2) is called stability function. This function is a polynomial when the selected one-step numerical method is explicit. R ⁢ ( ν ) is a rational function (a ratio of two polynomials) when implicit one-step numerical methods are used.

The selected one-step method will produce a bounded sequence of approximations { y 1 , y 2 , … , y n , … } to the solution of (5.1) for the applied value h of the time-stepsize if the important inequality | R ⁢ ( ν ) | ≤ 1 is satisfied for ν = h ⁢ λ ∈ ℂ - , see (5.2). The one-step numerical method is called absolutely stable for this value of parameter ν ∈ ℂ - , while the set of all values ν ∈ ℂ - for which the relationship | R ⁢ ( ν ) | ≤ 1 holds is forming the absolute stability region of the chosen one-step numerical method ([9], see also [11, 13, 14, 20, 22, 26, 27]). It is important now to answer the following question: will it be possible to calculate the stability function R [ q ] ⁢ ( ν ) of some version of the Richardson Extrapolation obtained by q ≥ 0 if we select some underlying one-step method with a known stability function R ⁢ ( ν ) ? The answer to this important question is positive for q ≤ 7 . More precisely, the following theorem holds:

Theorem 5.1.

Assume that an arbitrary one-step numerical method with a stability function R ⁢ ( ν ) is used in the solution of (5.1). Then the stability functions R [ q ] ⁢ ( ν ) of any version of the Repeated Richardson Extrapolation with q = 0 , 1 , … , 7 can be calculated by using the formula

(5.3) R [ q ] ⁢ ( ν ) = P ~ [ q ] ⁢ ( ν ) S [ q ] , q = 0 , 1 , … , 7 ,

where S [ q ] are defined as in the previous section, while the functions P ~ [ q ] ⁢ ( ν ) are listed below for the eight values of parameter q used in (5.3):

(5.4) P ~ [ 0 ] ⁢ ( ν ) = 2 p ⁢ [ R ⁢ ( ν 2 ) ] 2 - R ⁢ ( ν ) ,

(5.5) P ~ [ 1 ] ⁢ ( ν ) = 2 2 ⁢ p + 1 ⁢ [ R ⁢ ( ν 4 ) ] 4 - 3 ⋅ 2 p ⁢ [ R ⁢ ( ν 2 ) ] 2 + R ⁢ ( ν ) ,

(5.6) P ~ [ 2 ] ⁢ ( ν ) = 2 3 ⁢ p + 3 ⁢ [ R ⁢ ( ν 8 ) ] 8 - 7 ⋅ 2 2 ⁢ p + 1 ⁢ [ R ⁢ ( ν 4 ) ] 4 + 7 ⋅ 2 p ⁢ [ R ⁢ ( ν 2 ) ] 2 - R ⁢ ( ν ) ,

(5.7) P ~ [ 3 ] ⁢ ( ν ) = 2 4 ⁢ p + 6 ⁢ [ R ⁢ ( ν 16 ) ] 16 - 15 ⋅ 2 3 ⁢ p + 3 ⁢ [ R ⁢ ( ν 8 ) ] 8 + 35 ⋅ 2 2 ⁢ p + 1 ⁢ [ R ⁢ ( ν 4 ) ] 4 - 15 ⋅ 2 p ⁢ [ R ⁢ ( ν 2 ) ] 2 + R ⁢ ( ν ) ,

P ~ [ 4 ] ⁢ ( ν ) = 2 5 ⁢ p + 10 ⁢ [ R ⁢ ( ν 32 ) ] 32 - 31 ⋅ 2 4 ⁢ p + 6 ⁢ [ R ⁢ ( ν 16 ) ] 16 + 155 ⋅ 2 3 ⁢ p + 3 ⁢ [ R ⁢ ( ν 8 ) ] 8

(5.8) - 155 ⋅ 2 2 ⁢ p + 1 ⁢ [ R ⁢ ( ν 4 ) ] 4 + 35 ⋅ 2 p ⁢ [ R ⁢ ( ν 2 ) ] 2 - R ⁢ ( ν ) ,

P ~ [ 5 ] ⁢ ( ν ) = 2 6 ⁢ p + 15 ⁢ [ R ⁢ ( ν 64 ) ] 64 - 63 ⋅ 2 5 ⁢ p + 10 ⁢ [ R ⁢ ( ν 32 ) ] 32 + 651 ⋅ 2 4 ⁢ p + 6 ⁢ [ R ⁢ ( ν 16 ) ] 16

(5.9) - 1395 ⋅ 2 3 ⁢ p + 3 ⁢ [ R ⁢ ( ν 8 ) ] 8 + 651 ⋅ 2 2 ⁢ p + 1 ⁢ [ R ⁢ ( ν 4 ) ] 4 - 63 ⋅ 2 p ⁢ [ R ⁢ ( ν 2 ) ] 2 + R ⁢ ( ν ) ,

P ~ [ 6 ] ⁢ ( ν ) = 2 7 ⁢ p + 21 ⁢ [ R ⁢ ( ν 128 ) ] 128 - 127 ⋅ 2 6 ⁢ p + 15 ⁢ [ R ⁢ ( ν 64 ) ] 64 + 2667 ⋅ 2 5 ⁢ p + 10 ⁢ [ R ⁢ ( ν 32 ) ] 32

- 11811 ⋅ 2 4 ⁢ p + 6 ⁢ [ R ⁢ ( ν 16 ) ] 16 + 11811 ⋅ 2 3 ⁢ p + 3 ⁢ [ R ⁢ ( ν 8 ) ] 8 - 2667 ⋅ 2 2 ⁢ p + 1 ⁢ [ R ⁢ ( ν 4 ) ] 4

(5.10) + 127 ⋅ 2 p ⁢ [ R ⁢ ( ν 2 ) ] 2 - R ⁢ ( ν ) ,

P ~ [ 7 ] ⁢ ( ν ) = 2 8 ⁢ p + 28 ⁢ [ R ⁢ ( ν 256 ) ] 256 - 255 ⋅ 2 7 ⁢ p + 21 ⁢ [ R ⁢ ( ν 128 ) ] 128 + 10795 ⋅ 2 6 ⁢ p + 15 ⁢ [ R ⁢ ( ν 64 ) ] 64

- 97155 ⋅ 2 5 ⁢ p + 10 ⁢ [ R ⁢ ( ν 32 ) ] 32 + 200787 ⋅ 2 4 ⁢ p + 6 ⁢ [ R ⁢ ( ν 16 ) ] 16 - 97155 ⋅ 2 3 ⁢ p + 3 ⁢ [ R ⁢ ( ν 8 ) ] 8

(5.11) + 10795 ⋅ 2 2 ⁢ p + 1 ⁢ [ R ⁢ ( ν 4 ) ] 4 - 255 ⋅ 2 p ⁢ [ R ⁢ ( ν 2 ) ] 2 + R ⁢ ( ν ) .

Proof.

We shall prove the theorem for the Four-times Repeated Richardson Extrapolation, i.e., when q = 4 . It is clear, see formula (5.2), that the following six relationships hold for the Four-times Repeated Richardson Extrapolation:

(5.12)

z n [ 0 ] = R ⁢ ( ν ) ⁢ y n - 1 [ 4 ] , z n [ 1 ] = [ R ⁢ ( ν 2 ) ] 2 ⁢ y n - 1 [ 4 ] , z n [ 2 ] = [ R ( ν 4 ) ] 4 y n - 1 [ 4 ] ,

z n [ 3 ] = [ R ⁢ ( ν 8 ) ] 8 ⁢ y n - 1 [ 4 ] , z n [ 4 ] = [ R ⁢ ( ν 16 ) ] 16 ⁢ y n - 1 [ 4 ] , z n [ 5 ] = [ R ( ν 32 ) ] 32 y n - 1 [ 4 ] .

Multiply equalities (5.12) with -1, 31 ⋅ 2 p , - 155 ⋅ 2 2 ⁢ p + 1 , 155 ⋅ 2 3 ⁢ p + 3 , - 31 ⋅ 2 4 ⁢ p + 6 and 2 5 ⁢ p + 10 , respectively. Form the sum of the modified relationships and divide both sides of the obtained equality with the expression 2 5 ⁢ p + 10 - 31 ⋅ 2 4 ⁢ p + 6 + 155 ⋅ 2 3 ⁢ p + 3 - 155 ⋅ 2 2 ⁢ p + 1 + 35 ⋅ 2 p - 1 to obtain

y n [ 4 ] = P ~ [ 4 ] ⁢ ( ν ) S [ 4 ] ⁢ y n - 1 [ 4 ] = R [ 4 ] ⁢ ( ν ) ⁢ y n - 1 [ 4 ] .

In this way the theorem is proved for q = 4 .The same approach can obviously be applied to prove the theorem in the other seven cases. ∎

Remark 5.2.

The above proof of Theorem 5.1 for q = 4 shows that the stability polynomial of the Four-times Repeated Richardson Extrapolation is expressed by the stability polynomial of the underlying one-step numerical method, but it is different, which means that the two numerical methods, the underlying one-step numerical method and its combination with the Four-times Repeated Richardson Extrapolation, will in general have different absolute stability properties. This implies that it is necessary to study carefully the absolute stability properties of the Four-times Repeated Richardson Extrapolation. The same conclusion can also be drawn for the other versions of the Richardson Extrapolation when q ≤ 7 .

Conjecture 3.

The assertion of Theorem 5.1 remains true for q > 7 .

Remark 5.3.

The stability results presented in this section were obtained by using, as in [9], the linear and scalar problem (5.1). Many attempts to generalize the original Dahlquist stability theory have been carried out after 1963 (see, for example, [6, 13, 14, 20, 22]). It can be shown that many of these extended stability definitions and results remain valid also for the different versions of the Richardson Extrapolation.

6 Using Explicit Runge–Kutta Methods with the Richardson Extrapolation

The class of the Runge–Kutta Methods (RKMs), introduced more than 100 years ago in [19, 21, 25] (by K. Heun, W. Kutta and C. Runge, respectively), is a sub-class of the one-step methods. This means that all results proved in the previous sections for the one-step methods are also valid for the RKMs.

Explicit Runge–Kutta Methods (ERKMs) will be used in the remaining part of this paper in combination with the advanced versions of the Richardson Extrapolation to handle numerically a two-parameter family of numerical examples.

6.1 Introduction of the ERKMs

The m-stage Explicit Runge–Kutta Methods are one-step numerical methods for solving systems of ODEs (2.1) that are based on the following formula (see more details in [6, 13] and [22]):

y n = y n - 1 + h ⁢ ∑ i = 1 m c i ⁢ k i n .

The coefficients c i are constants, while the stage vectors k i n are defined by

k 1 n = f ⁢ ( t n - 1 , y n - 1 ) , k i n = f ⁢ ( t n - 1 + h ⁢ a i , y n - 1 + h ⁢ ∑ j = 1 i - 1 b i ⁢ j ⁢ k j n ) , a i = ∑ j = 1 i - 1 b i ⁢ j ,

where i = 2 , 3 , … , m and b i ⁢ j are constants depending on the particular ERKM.

6.2 Stability Polynomials of Some ERKMs

We shall assume that the order of accuracy p of the chosen ERKM is equal to the number m of the stage vectors k i n , i = 1 , 2 , … , m . This is possible only if m ≤ 4 (the relationship p < m is always satisfied when m > 4 ). If p = m , then the stability polynomial R ⁢ ( ν ) associated with the selected Explicit Runge–Kutta Method is given (see again [22, p. 202]) by

(6.1) R ⁢ ( ν ) = 1 + ν + ν 2 2 ! + ν 3 3 ! + … + ν p p ! , p = m , m = 1 , 2 , 3 , 4 .

The stability polynomials R [ q ] ⁢ ( ν ) related to the Classical Richardson Extrapolation when q = 0 or to different versions of the Repeated Richardson Extrapolation when 1 ≤ q ≤ 7 can easily be obtained by substituting the expression of R ⁢ ( ν ) from (6.1) on the right-hand sides of formulae (5.4)–(5.11).

6.3 Absolute Stability Regions of ERKMs with p = m

The absolute stability regions of the ERKMs when p = m and m = 1 , 2 , 3 , 4 will be compared with the absolute stability regions of the corresponding new numerical methods which are obtained when any of these ERKMs is combined with the eight advanced versions of the Richardson Extrapolation that were introduced and discussed in the previous sections. A straight-forward algorithm which can efficiently be used for obtaining the absolute stability region of an arbitrary Runge–Kutta Method (RKM), not only of an Explicit Runge–Kutta method (ERKM), is presented below.

Algorithm 1.

Let ν = α ¯ + β ¯ ⁢ i with α ¯ ≤ 0 and suppose that ε is a small increment. Set α ¯ = 0 and calculate the values of R ⁢ ( ν ) of the chosen RKM for α ¯ = 0 and β ¯ = 0 , ε , 2 ⁢ ε , 3 ⁢ ε , … . The computations are carried on as long as both

| R ⁢ ( ν ) | ≤ 1 ,
β ¯ ≤ 10 5

and the last value of β ¯ obtained in this way is denoted by β ¯ 0 . Continue the process by setting α ¯ = - ε and by repeating the computations for the new value of α ¯ and for β ¯ = 0 , ε , 2 ⁢ ε , 3 ⁢ ε , … . Let β ¯ 1 be the largest value of β ¯ for which requirements (a) and (b) are satisfied. Set α ¯ equal to - 2 ⁢ ε , - 3 ⁢ ε , … and for each α ¯ repeat the computations described above. Continue always this process until (a), (b) and

α ¯ ≤ 10 5

hold. A large set of points { ( 0 , β ¯ 0 ) , ( - ε , β ¯ 1 ) , ( - 2 ⁢ ε , β ¯ 2 ) , … } located in the negative part of the complex plane and over the real axis will be calculated in this manner. Moreover, all these points are inside a part of the absolute stability region which is located to the left of the imaginary axis and over the negative part of the real axis. Therefore, we shall obtain, by connecting successively these points, a domain in which the selected RKM will produce stable results for the chosen stepsize when (5.1) is handled.

Remark 6.1.

If the selected RKM is explicit, then there is as a rule no need to include requirements (b) and (c) in the definition of Algorithm 1 and the boundary of the part of the absolute stability region which is located to the left of the imaginary axis and over the negative part of the real axis will be obtained when the increment ε is sufficiently small. If the RKM is implicit (such methods will be studied in a paper which is in preparation), then Algorithm 1 with the application of requirements (b) and (c) may result in a huge domain in which the computational process will be stable when (5.1) is solved (this domain will be a part of the absolute stability region in this case). If the computations in Algorithm 1 are always stopped by requirements (b) and (c), then a huge square with vertices ( 0.0 , 0.0 ) , ( 0.0 , 10 5 ) , ( - 10 5 , 10 5 ) and ( - 10 5 , 0.0 ) which is a part of the absolute stability region of the selected method will be obtained. It should be mentioned that the constant 10 5 in requirements (b) and (c) can be replaced with some other large positive number.

The absolute stability regions are symmetric with regard to the real axis, [22], and there is no need to repeat the computations for negative values of the imaginary part β ¯ of ν = h ⁢ λ = α ¯ + β ¯ ⁢ i .

There is no need to require bounded approximations of the solution if α ¯ > 0 , because for positive values of α ¯ the exact solution of (5.1) is unbounded. Therefore, it is relevant to introduce the requirement | R ⁢ ( ν ) | ≤ 1 and to study the absolute stability properties of the numerical methods only when α ¯ ≤ 0 .

The results obtained for the ERKMs with p = m and m = 1 , 2 , 3 , 4 by using Algorithm 1 are given in the four plots in Figure 1 (the boundaries of the absolute stability regions of the ERKMs are the red curves in these plots).

Figure 1
Absolute stability regions of (a) the One-stage First-order ERKM (the Forward Euler Formula) and its combinations with eight versions of the Richardson Extrapolation (the upper left plot), (b) any particular method from the class of the Two-stages Second-order ERKMs and its combinations with eight versions of the Richardson Extrapolation (upper right plot) (c) any particular method from the class of the Three-stages Third-order ERKMs and its combinations with eight versions of the Richardson Extrapolation (lower left plot) (d) any particular method from the class of the Four-stages Fourth-order ERKMs and its combinations with eight versions of the Richardson Extrapolation (lower right plot).

Figure 1

Absolute stability regions of (a) the One-stage First-order ERKM (the Forward Euler Formula) and its combinations with eight versions of the Richardson Extrapolation (the upper left plot), (b) any particular method from the class of the Two-stages Second-order ERKMs and its combinations with eight versions of the Richardson Extrapolation (upper right plot) (c) any particular method from the class of the Three-stages Third-order ERKMs and its combinations with eight versions of the Richardson Extrapolation (lower left plot) (d) any particular method from the class of the Four-stages Fourth-order ERKMs and its combinations with eight versions of the Richardson Extrapolation (lower right plot).

The same procedure can be used to obtain the absolute stability regions of all versions of the Richardson Extrapolation. It is only necessary to replace the stability polynomial R ⁢ ( ν ) with the stability polynomials R [ q ] ⁢ ( ν ) , q = 0 , 1 , … , 7 . The boundaries of the absolute stability regions of all eight versions of the Richardson Extrapolation are also given in Figure 1.

The results shown in Figure 1 can be considered as a graphical proof of the following theorem:

Theorem 6.2.

The absolute stability regions of the numerical methods that are combinations of Explicit Runge–Kutta Methods (ERKMs) with different versions of the Richardson Extrapolation are always larger than the absolute stability regions of the underlying Explicit Runge–Kutta Methods when the conditions p = m , m = 1 , 2 , 3 , 4 and 0 ≤ q ≤ 7 are satisfied.

Remark 6.3.

The above stability results are valid only when equation (5.1) is handled. If (5.1) is solved, then the computational process remains stable when the stepsize is chosen so that the point ν = h ⁢ λ is located inside the absolute stability region of the selected numerical procedure. The stability results can easily be generalized (see, for example, [22]) for the case where the scalar parameter λ in (5.1) is replaced with a constant matrix A ∈ ℝ s × s all eigenvalues of which are different and with non-positive real parts. It is necessary in this case to introduce the requirement that all points ν i = h ⁢ λ i , i = 1 , 2 , … , s , where λ i , i = 1 , 2 , … , s , are eigenvalues of matrix A, are located inside the absolute stability region of the chosen version of the Richardson Extrapolation. The non-linear case is much more complicated, but one should expect that the computations will be stable for a given h if all eigenvalues λ i ⁢ ( t ) , i = 1 , 2 , … , s , of the Jacobian matrix J = ∂ ⁡ f ∂ ⁡ y are with non-positive real parts and if all ν i ⁢ ( t ) = h ⁢ λ i ⁢ ( t ) , i = 1 , 2 , … , s , are inside the absolute stability region for all t ∈ [ a , b ] . Much more details about this topic can be found in [2, 6, 13, 14, 20] and/or [22].

6.4 Selecting Particular ERKMs

Some particular ERKMs satisfying the conditions p = m and m = 1 , 2 , 3 , 4 are needed for the numerical experiments. Only the Forward Euler Formula satisfies the requirement p = m = 1 , but there exists a large class of ERKMs for each pair { p , m } with p = m and m = 2 , 3 , 4 (depending on one or two free parameters for m = 2 and m = 3 , 4 , respectively; see [22]). All methods from such a class have the same absolute stability region and each of these methods can be used in combination with any version of the Richardson Extrapolation but if the version is fixed (i.e., if q is fixed), then all resulting combinations have the same absolute stability region. Methods for each pair { p , m } with p = m and m = 2 , 3 , 4 are listed below.

6.4.1 The Forward Euler Formula

The ERKM with p = m = 1 can be written in the following form:

(6.2) y n = y n - 1 + h ⁢ f ⁢ ( t n - 1 , y n - 1 ) .

6.4.2 The Improved Euler Formula

The ERKM with p = m = 2 can be written as

(6.3)

k 1 = f ⁢ ( t n - 1 , y n - 1 ) ,

k 2 = f ⁢ ( t n - 1 + h , y n - 1 + h ⁢ k 1 ) ,

y n = y n - 1 + 1 2 ⁢ h ⁢ ( k 1 + k 2 ) .

6.4.3 The Heun’s Method

The ERKM with p = m = 3 can be written as

(6.4)

k 1 = f ⁢ ( t n - 1 , y n - 1 ) ,

k 2 = f ⁢ ( t n - 1 + 1 3 ⁢ h , y n - 1 + 1 3 ⁢ h ⁢ k 1 ) ,

k 3 = f ⁢ ( t n - 1 + 2 3 ⁢ h , y n - 1 + 2 3 ⁢ h ⁢ k 2 ) ,

y n = y n - 1 + 1 4 ⁢ h ⁢ ( k 1 + 3 ⁢ k 3 ) .

6.4.4 The Most Popular Runge–Kutta Method

This method (an ERKM with p = m = 4 ) is given below:

(6.5)

k 1 = f ⁢ ( t n - 1 , y n - 1 ) ,

k 2 = f ⁢ ( t n - 1 + 1 2 ⁢ h , y n - 1 + 1 2 ⁢ h ⁢ k 1 ) ,

k 3 = ( t n - 1 + 1 2 ⁢ h , y n - 1 + 1 2 ⁢ h ⁢ k 2 ) ,

k 4 = f ⁢ ( t n - 1 + h , y n - 1 + h ⁢ k 3 ) ,

y n = y n - 1 + 1 4 ⁢ h ⁢ ( k 1 + 2 ⁢ k 2 + 2 ⁢ k 3 + k 4 ) .

7 Numerical Examples

7.1 Introduction of a Two-Parameter Family for Testing Both the Accuracy and the Absolute Stability Properties

Consider the system of three linear ODEs:

(7.1) d ⁢ y d ⁢ t = A ⁢ y , t ∈ [ 0 , 13.1072 ] , y = ( y 1 , y 2 , y 3 ) T , y ⁢ ( 0 ) = ( 1 , 0 , 2 ) T , A ∈ ℝ 3 × 3

The elements of matrix A from (7.1) are given below:

(7.2) a 11 = - γ - β - 0.6 , a 12 = - γ - 0.3 , a 13 = γ + β + 0.3 ,

(7.3) a 21 = γ - 2 ⁢ β + 0.3 , a 22 = γ - β , a 23 = - γ + β - 0.3 ,

(7.4) a 31 = γ - 3 ⁢ β - 0.3 , a 32 = - γ - β - 0.3 , a 33 = γ + 2 ⁢ β .

The three components of the exact solution of the problem defined by (7.1)–(7.4) are given by

y 1 ⁢ ( t ) = e - 0.3 ⁢ t ⁢ sin ⁡ β ⁢ t + e γ ⁢ t ,

y 2 ⁢ ( t ) = e - 0.3 ⁢ t ⁢ cos ⁡ β ⁢ t + e γ ⁢ t ,

y 3 ⁢ ( t ) = e - 0.3 ⁢ t ⁢ ( sin ⁡ β ⁢ t + cos ⁡ β ⁢ t ) + e γ ⁢ t .

The curves representing the three components of the solution for γ = - 750 and β = 8 are given in Figure 2. It should be mentioned here that the eigenvalues of matrix A from (7.1)–(7.4) are

μ 1 = γ ≤ 0 ,

μ 2 = - 0.3 + β ⁢ i ,

μ 3 = - 0.3 - β ⁢ i .

In order to ensure stable computations in the solution of this example, one should require all three points h ⁢ μ 1 , h ⁢ μ 2 and h ⁢ μ 3 to be inside the absolute stability region of the used method.

$Figure 2 The components of the solution of the system of ODEs represented by (7.1)–(7.4) with γ = - 750 {\gamma=-750} and β = 8 {\beta=8} .$

Figure 2

The components of the solution of the system of ODEs represented by (7.1)–(7.4) with γ = - 750 and β = 8 .

It is clear that by increasing the absolute value of γ and keeping the value of β small, we are increasing the stability requirements. Moreover, the stability interval on the negative part of the real axis is important in this case ( h ⁢ μ 1 should be inside the absolute stability interval on the negative part of the real axis). Increasing the value of β leads to an increase of the oscillations and, thus to an increase of the accuracy requirements. However, large values of β will also cause stability problems related to the height of the absolute stability region (the influence of the complex eigenvalues will become dominant when β is much larger in absolute value than γ). This short analysis shows clearly that it is possible to steer both the accuracy requirements and the stability requirements by selecting appropriate values of each of the two parameters.

7.2 Organization of the Computations

We used the four ERKMs listed in (6.2)–(6.5) and their combinations with all eight versions of the Richardson Extrapolation. Twelve runs with different stepsizes were successively performed for each ERKM and for each version of the Richardson Extrapolation. The first run was always performed with a stepsize h = 0.02048 and we halved the stepsize after each run. This means that the last (and also the smallest) stepsize was h = 10 - 5 . The integration interval [ 0 , 13.1072 ] was divided into 128 equal sub-intervals and the accuracy of the results obtained by any of the selected numerical methods was evaluated at the end of each sub-interval. Assume that

t ¯ j , j ∈ { 1 , 2 , … , 128 } , are the end-point of any of the 128 sub-intervals,
y ¯ i ⁢ j ≈ y i ⁢ ( t ¯ j ) is an approximation calculated by any of the used numerical methods.

Then the following formula was used to evaluate the accuracy achieved by the selected method at this point:

ERROR j = ∑ t = 1 3 ( y i ⁢ ( t ¯ j ) - y ¯ i ⁢ j ) 2 max ⁡ [ ∑ t = 1 3 ( y i ⁢ ( t ¯ j ) ) 2 , 1.0 ] .

The global error is computed by using the formula

ERROR = max j = 1 , 2 , … , 128 ⁡ ( ERROR j ) .

7.3 Moderate Accuracy and Absolute Stability Requirements

The ERKMs defined by (6.2)–(6.5) were used in the solution of the example defined with setting with β = 32 and γ = - 750 in (7.1)–(7.4). The errors are given in Tables 1–4, while the corresponding convergence rates can be found in Tables 5–8. It should be mentioned here that a requirement to achieve at least two correct significant digits in the calculated solution was imposed in all tables given below, which implies that ERROR must be less than 10 - 2 . Four important conclusions can be drawn from Tables 1–8:

The results are becoming more and more accurate when q is increased. The repeated versions of the Richardson Extrapolation are sometimes so accurate that the rounding errors are starting to interfere with the truncation errors in spite of the fact that extended precision (with 32 significant digits) was used in all runs.
The expected convergence rates are nearly always achieved.
The repeated versions of the Repeated Richardson Extrapolation are sometimes stable in runs with large stepsizes for which the underlying ERKMs are unstable.
The expected order of accuracy of the Two-times Repeated Richardson Extrapolation applied in the combination with the Two-stages Second-order ERKM is five. The order which is actually achieved is clearly six (which is seen from the results given in the column under “Two-times REP” in Table 6). This means that the coefficient of the first term in 𝒪 ⁢ ( h 5 ) is equal to zero for this particular example. Similar situation arises when the Six-times Repeated Richardson Extrapolation is run together with the Two-stages Second-order ERKM (the actually achieved order of accuracy for his particular example is approximately ten (not nine as it should be).

We use the following notation in Tables 1–8:

“N.S.” means that the method is not stable.
“N.A.” means that the error is greater than 10 - 2 .
“V.L.” means that the rate is very large (due presumably to the fact that ν is very close to the boundary of the absolute stability region).
The rounding errors are dominant in the grey zone.

Table 1

Accuracy results obtained in solving the test-problem with β = 32 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The First-order One-stage Explicit Runge–Kutta Method (the Explicit Euler Formula) and its combinations with eight versions of the Richardson Extrapolation are run.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	9.4E - 10	2.2E - 12
N.S.	N.S.	N.S.	N.S.	4.2E - 05	3.2E - 07	1.6E - 09	3.7E - 12	4.3E - 15
N.S.	N.S.	N.S.	1.5E - 03	1.3E - 06	5.6E - 09	1.3E - 11	1.4E - 14	8.4E - 18
N.A.	N.A.	1.2E - 03	9.3E - 06	4.1E - 08	8.8E - 11	1.0E - 13	5.6E - 17	1.6E - 20
N.A.	N.A.	1.5E - 04	5.8E - 07	1.3E - 09	1.4E - 12	7.8E - 16	2.2E - 19	3.2E - 23
N.A.	4.6E - 03	1.8E - 05	3.7E - 08	4.0E - 11	2.1E - 14	6.1E - 18	8.5E - 22	6.2E - 26
N.A.	1.2E - 03	2.3E - 06	2.3E - 09	1.2E - 12	3.3E - 16	4.8E - 20	3.3E - 24	1.2E - 28
N.A.	3.0E - 04	2.8E - 07	1.4E - 10	3.9E - 14	5.2E - 18	3.7E - 22	1.3E - 26	An error in the conversion from LaTeX to XML has occurred here. gray5.8E - 31
N.A.	7.3E - 05	3.5E - 08	8.9E - 12	1.2E - 15	8.2E - 20	2.9E - 24	5.1E - 29	An error in the conversion from LaTeX to XML has occurred here. gray7.6E - 31
N.A.	1.8E - 05	4.4E - 09	5.6E - 13	3.8E - 17	1.3E - 21	2.3E - 26	An error in the conversion from LaTeX to XML has occurred here. gray5.5E - 31	An error in the conversion from LaTeX to XML has occurred here. gray1.6E - 30
N.A.	4.5E - 06	5.5E - 10	3.5E - 14	1.2E - 18	2.0E - 23	1.8E - 28	An error in the conversion from LaTeX to XML has occurred here. gray1.1E - 30	An error in the conversion from LaTeX to XML has occurred here. gray2.0E - 30
N.A.	1.1E - 06	6.9E - 11	2.2E - 15	3.7E - 20	3.1E - 25	1.6E - 30	An error in the conversion from LaTeX to XML has occurred here. gray1.5E - 30	An error in the conversion from LaTeX to XML has occurred here. gray1.8E - 30

Table 2

Accuracy results obtained in solving the test-problem with β = 32 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The Second-order Two-stage Explicit Runge–Kutta Method (the Improved Euler Formula) and its combinations with eight versions of the Richardson Extrapolation are run.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	1.6E - 05	1.3E - 14	7.8E - 17
N.S.	N.S.	N.S.	N.S.	4.6E - 08	2.1E - 10	2.3E - 13	1.1E - 17	7.6E - 20
N.S.	6.2E - 03	4.3E - 05	1.7E - 08	7.2E - 10	1.6E - 12	9.0E - 16	1.1E - 20	7.5E - 23
N.A.	7.7E - 04	2.7E - 06	2.6E - 10	1.1E - 11	1.3E - 14	3.5E - 18	1.1E - 23	7.3E - 26
N.A.	9.7E - 05	1.7E - 07	4.1E - 12	1.8E - 13	9.9E - 17	1.4E - 20	1.1E - 26	An error in the conversion from LaTeX to XML has occurred here. gray7.1E - 29
4.6E - 03	1.2E - 05	1.0E - 08	6.5E - 14	2.8E - 15	7.7E - 19	5.4E - 23	An error in the conversion from LaTeX to XML has occurred here. gray1.1E - 29	An error in the conversion from LaTeX to XML has occurred here. gray1.3E - 31
1.2E - 03	1.5E - 06	6.5E - 10	1.0E - 15	4.3E - 17	6.0E - 21	2.1E - 25	An error in the conversion from LaTeX to XML has occurred here. gray3.4E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.0E - 31
3.0E - 04	1.9E - 07	4.1E - 11	1.6E - 17	6.7E - 19	4.7E - 23	8.2E - 28	An error in the conversion from LaTeX to XML has occurred here. gray2.7E - 31	An error in the conversion from LaTeX to XML has occurred here. gray1.9E - 31
7.3E - 05	2.4E - 08	2.6E - 12	2.5E - 19	1.1E - 20	3.7E - 25	3.2E - 30	An error in the conversion from LaTeX to XML has occurred here. gray4.2E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.4E - 31
1.8E - 05	3.0E - 09	1.6E - 13	3.9E - 21	1.6E - 22	2.9E - 27	An error in the conversion from LaTeX to XML has occurred here. gray3.0E - 31	An error in the conversion from LaTeX to XML has occurred here. gray3.1E - 31	An error in the conversion from LaTeX to XML has occurred here. gray6.4E - 31
4.5E - 06	3.7E - 10	1.0E - 14	6.0E - 23	2.6E - 24	2.2E - 29	An error in the conversion from LaTeX to XML has occurred here. gray6.1E - 31	An error in the conversion from LaTeX to XML has occurred here. gray4.4E - 31	An error in the conversion from LaTeX to XML has occurred here. gray1.2E - 30
1.1E - 06	4.6E - 11	6.2E - 16	9.4E - 25	4.0E - 26	An error in the conversion from LaTeX to XML has occurred here. gray7.2E - 31	An error in the conversion from LaTeX to XML has occurred here. gray5.4E - 31	An error in the conversion from LaTeX to XML has occurred here. gray6.7E - 31	An error in the conversion from LaTeX to XML has occurred here. gray9.4E - 31

Table 3

Accuracy results obtained in solving the linear test-problem with β = 32 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The Third-order Three-stage Explicit Runge–Kutta Method (Heun’s Formula) and its combinations with eight versions of the Richardson Extrapolation are run.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	2.1E - 11	5.2E - 17	1.6E - 20
N.S.	N.S.	N.S.	2.5E - 04	7.2E - 11	1.6E - 13	1.5E - 16	5.1E - 20	7.7E - 24
N.S.	7.4E - 03	4.0E - 07	5.1E - 10	5.6E - 13	6.1E - 16	2.9E - 19	5.0E - 23	3.7E - 27
1.6E - 03	7.1E - 06	1.2E - 08	8.0E - 12	4.4E - 15	2.4E - 18	5.7E - 22	4.9E - 26	An error in the conversion from LaTeX to XML has occurred here. gray1.9E - 30
1.9E - 04	4.5E - 07	3.9E - 10	1.3E - 13	3.4E - 17	9.3E - 21	1.1E - 24	An error in the conversion from LaTeX to XML has occurred here. gray4.8E - 29	An error in the conversion from LaTeX to XML has occurred here. gray2.4E - 31
2.4E - 05	2.8E - 08	1.2E - 11	2.0E - 15	2.7E - 19	3.6E - 23	2.2E - 27	An error in the conversion from LaTeX to XML has occurred here. gray2.5E - 31	An error in the conversion from LaTeX to XML has occurred here. gray1.0E - 31
3.0E - 06	1.7E - 09	3.8E - 13	3.1E - 17	2.1E - 21	1.4E - 25	4.2E - 30	An error in the conversion from LaTeX to XML has occurred here. gray1.3E - 31	An error in the conversion from LaTeX to XML has occurred here. gray3.2E - 31
3.8E - 07	1.1E - 10	1.2E - 14	4.8E - 19	1.6E - 23	5.5E - 28	An error in the conversion from LaTeX to XML has occurred here. gray2.0E - 31	An error in the conversion from LaTeX to XML has occurred here. gray4.4E - 31	An error in the conversion from LaTeX to XML has occurred here. gray3.2E - 31
4.3E - 08	6.8E - 12	3.7E - 16	7.5E - 21	1.3E - 25	2.2E - 30	An error in the conversion from LaTeX to XML has occurred here. gray3.3E - 31	An error in the conversion from LaTeX to XML has occurred here. gray1.6E - 31	An error in the conversion from LaTeX to XML has occurred here. gray4.9E - 31
5.9E - 09	4.2E - 13	1.2E - 17	1.2E - 22	9.9E - 28	An error in the conversion from LaTeX to XML has occurred here. gray3.5E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.5E - 31	An error in the conversion from LaTeX to XML has occurred here. gray5.1E - 3	An error in the conversion from LaTeX to XML has occurred here. gray2.4E - 31
7.4E - 10	2.7E - 14	3.6E - 19	1.8E - 24	An error in the conversion from LaTeX to XML has occurred here. gray8.0E - 30	An error in the conversion from LaTeX to XML has occurred here. gray1.2E - 31	An error in the conversion from LaTeX to XML has occurred here. gray5.6E - 31	An error in the conversion from LaTeX to XML has occurred here. gray3.0E - 31	An error in the conversion from LaTeX to XML has occurred here. gray1.1E - 30
9.2E - 11	1.7E - 15	1.1E - 20	2.9E - 26	An error in the conversion from LaTeX to XML has occurred here. gray4.4E - 31	An error in the conversion from LaTeX to XML has occurred here. gray5.5E - 31	An error in the conversion from LaTeX to XML has occurred here. gray4.4E - 3	An error in the conversion from LaTeX to XML has occurred here. gray1.2E - 31	An error in the conversion from LaTeX to XML has occurred here. gray1.3E - 30

Table 4

Accuracy results obtained in solving the test-problem with β = 32 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The Fourth-order Four-stage Explicit Runge–Kutta Method (the most popular Runge–Kutta formula) and its combinations with eight versions of the Richardson Extrapolation are run.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.A.	4.3E - 09	3.1E - 17	1.9E - 21	3.3E - 24
N.S.	N.S.	N.S.	2.7E - 10	1.4E - 13	1.4E - 17	3.0E - 20	9.0E - 24	8.1E - 28
N.S.	1.9E - 06	3.1E - 09	2.1E - 12	5.5E - 16	1.3E - 20	2.0E - 23	4.4E - 27	An error in the conversion from LaTeX to XML has occurred here. gray1.8E - 31
2.5E - 05	5.8E - 08	4.8E - 11	1.6E - 14	2.2E - 18	1.3E - 23	2.9E - 26	An error in the conversion from LaTeX to XML has occurred here. gray2.0E - 30	An error in the conversion from LaTeX to XML has occurred here. gray2.2E - 31
1.6E - 06	1.8E - 09	7.5E - 13	1.3E - 16	8.4E - 21	1.3E - 26	An error in the conversion from LaTeX to XML has occurred here. gray2.8E - 29	An error in the conversion from LaTeX to XML has occurred here. gray1.8E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.6E - 31
9.7E - 08	5.6E - 11	1.2E - 14	9.9E - 19	3.3E - 23	1.2E - 29	An error in the conversion from LaTeX to XML has occurred here. gray1.8E - 31	An error in the conversion from LaTeX to XML has occurred here. gray1.9E - 31	An error in the conversion from LaTeX to XML has occurred here. gray1.9E - 31
6.1E - 09	1.8E - 12	1.8E - 16	7.7E - 21	1.3E - 25	An error in the conversion from LaTeX to XML has occurred here. gray1.7E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.3E - 31	An error in the conversion from LaTeX to XML has occurred here. gray1.9E - 31	An error in the conversion from LaTeX to XML has occurred here. gray1.7E - 31
3.8E - 10	5.5E - 14	2.9E - 18	6.1E - 23	5.0E - 28	An error in the conversion from LaTeX to XML has occurred here. gray2.4E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.1E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.0E - 31	An error in the conversion from LaTeX to XML has occurred here. gray3.5E - 31
2.4E - 11	1.7E - 15	4.5E - 20	4.7E - 25	2.0E - 30	An error in the conversion from LaTeX to XML has occurred here. gray2.5E - 31	An error in the conversion from LaTeX to XML has occurred here. gray1.8E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.8E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.5E - 31
1.5E - 12	5.4E - 17	7.0E - 22	3.7E - 27	An error in the conversion from LaTeX to XML has occurred here. gray2.5E - 31	An error in the conversion from LaTeX to XML has occurred here. gray3.5E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.1E - 31	An error in the conversion from LaTeX to XML has occurred here. gray3.6E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.7E - 31
9.3E - 14	1.7E - 18	1.1E - 23	2.9E - 29	An error in the conversion from LaTeX to XML has occurred here. gray3.2E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.3E - 31	An error in the conversion from LaTeX to XML has occurred here. gray4.2E - 31	An error in the conversion from LaTeX to XML has occurred here. gray4.2E - 31	An error in the conversion from LaTeX to XML has occurred here. gray3.9E - 32
5.8E - 14	5.2E - 20	1.8E - 25	An error in the conversion from LaTeX to XML has occurred here. gray6.6E - 31	An error in the conversion from LaTeX to XML has occurred here. gray2.2E - 31	An error in the conversion from LaTeX to XML has occurred here. gray4.8E - 31	An error in the conversion from LaTeX to XML has occurred here. gray3.1E - 31	An error in the conversion from LaTeX to XML has occurred here. gray6.6E - 31	An error in the conversion from LaTeX to XML has occurred here. gray8.6E - 32

Table 5

Convergence rates obtained in solving the test-problem with β = 32 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The First-order One-stage Explicit Runge–Kutta Method (the Forward Euler Formula) and its combinations with eight versions of the Richardson Extrapolation are run. The perfect convergence rates are 2 , 4 , 8 , 16 , 32 , 64 , 128 , 256 and 512.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	–	–
N.S.	N.S.	N.S.	N.S.	–	–	–	256.01	515.15
N.S.	N.S.	N.S.	–	32.15	63.86	128.13	255.94	513.16
N.A.	N.A.	–	158.63	32.02	63.89	128.39	256.27	512.88
N.A.	N.A.	8.00	15.98	31.97	63.87	128.15	255.96	512.58
N.A.	–	8.01	16.00	32.07	64.02	127.96	256.17	511.25
N.A.	4.00	7.97	16.01	31.94	64.07	128.00	255.56	509.84
N.A.	3.92	8.02	15.94	32.04	63.98	128.03	256.15	An error in the conversion from LaTeX to XML has occurred here. gray209.98
N.A.	4.08	7.99	16.05	31.98	64.05	127.93	254.90	An error in the conversion from LaTeX to XML has occurred here. gray0.76
N.A.	4.01	7.99	16.00	32.10	64.17	128.32	An error in the conversion from LaTeX to XML has occurred here. gray92.39	An error in the conversion from LaTeX to XML has occurred here. gray0.48
N.A.	4.00	8.01	16.01	31.95	63.83	127.68	An error in the conversion from LaTeX to XML has occurred here. gray 0.49	An error in the conversion from LaTeX to XML has occurred here. gray0.79
N.A.	4.01	8.00	15.96	31.98	63.99	110.63	An error in the conversion from LaTeX to XML has occurred here. gray0.77	An error in the conversion from LaTeX to XML has occurred here. gray1.12

Table 6

Convergence rates obtained in solving the test-problem with β = 32 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The selected Second-order Two-stage Explicit Runge–Kutta Method (the Improved Euler Formula) and its combinations with eight versions of the Richardson Extrapolation are run. The perfect convergence rates are 4 , 8 , 16 , 32 , 64 , 128 , 256 , 512 and 1024.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	–	–	–
N.S.	N.S.	N.S.	N.S.	–	–	V.L.	1175.44	1019.63
N.S.	–	–	–	63.89	128.40	255.56	1027.03	1020.03
N.A.	7.98	16.03	64.02	63.72	127.56	255.68	1018.35	1023.22
N.A.	8.00	15.99	63.92	64.20	128.54	256.93	1028.30	1014.03
4.00	7.99	16.06	63.93	63.77	127.98	255.60	1009.52	An error in the conversion from LaTeX to XML has occurred here. gray544.27
3.92	8.01	15.95	63.96	64.04	128.03	255.24	An error in the conversion from LaTeX to XML has occurred here. gray31.25	An error in the conversion from LaTeX to XML has occurred here. gray0.65
4.08	7.99	16.02	63.92	64.04	128.03	256.72	An error in the conversion from LaTeX to XML has occurred here. gray1.24	An error in the conversion from LaTeX to XML has occurred here. gray1.06
4.01	8.01	15.96	64.23	64.10	127.99	256.43	An error in the conversion from LaTeX to XML has occurred here. gray0.65	An error in the conversion from LaTeX to XML has occurred here. gray0.78
4.00	8.00	16.04	63.90	63.02	128.22	An error in the conversion from LaTeX to XML has occurred here. gray10.60	An error in the conversion from LaTeX to XML has occurred here. gray1.34	An error in the conversion from LaTeX to XML has occurred here. gray0.38
4.01	7.99	15.98	63.95	63.81	128.70	An error in the conversion from LaTeX to XML has occurred here. gray0.49	An error in the conversion from LaTeX to XML has occurred here. gray0.71	An error in the conversion from LaTeX to XML has occurred here. gray0.54
4.00	8.00	16.00	64.04	64.09	An error in the conversion from LaTeX to XML has occurred here. gray30.80	An error in the conversion from LaTeX to XML has occurred here. gray 1.14	An error in the conversion from LaTeX to XML has occurred here. gray0.65	An error in the conversion from LaTeX to XML has occurred here. gray1.28

Table 7

Convergence rates obtained in solving the test-problem with β = 32 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The selected Third-order Three-stage Explicit Runge–Kutta Method (the Heun’s Formula) and its combinations with eight versions of the Richardson Extrapolation are run. The perfect convergence rates are 8 , 16 , 32 , 64 , 128 , 256 , 512 , 1024 and 2048.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	–	–	–
N.S.	N.S.	N.S.	–	–	–	V.L.	1013.70	2044.56
N.S.	–	–	V.L.	128.39	255.35	513.61	1019.96	2045.58
–	V.L.	31.85	63.72	128.44	255.04	512.20	1022.45	1963.16
8.03	15.98	32.12	64.16	128.24	255.91	512.50	1027.25	An error in the conversion from LaTeX to XML has occurred here. gray7.92
7.98	16.01	32.17	63.78	127.82	256.20	511.42	An error in the conversion from LaTeX to XML has occurred here. gray193.90	An error in the conversion from LaTeX to XML has occurred here. gray1.24
8.01	15.98	32.10	64.05	128.50	255.63	526.44	An error in the conversion from LaTeX to XML has occurred here. gray1.95	An error in the conversion from LaTeX to XML has occurred here. gray0.61
7.99	15.96	31.95	63.88	127.78	261.73	An error in the conversion from LaTeX to XML has occurred here. gray19.81	An error in the conversion from LaTeX to XML has occurred here. gray0.30	An error in the conversion from LaTeX to XML has occurred here. gray0.98
8.01	16.05	32.07	64.04	127.56	254.13	An error in the conversion from LaTeX to XML has occurred here. gray0.64	An error in the conversion from LaTeX to XML has occurred here. gray2.77	An error in the conversion from LaTeX to XML has occurred here. gray0.65
8.00	16.01	32.00	63.93	128.41	An error in the conversion from LaTeX to XML has occurred here. gray6.19	An error in the conversion from LaTeX to XML has occurred here. gray1.34	An error in the conversion from LaTeX to XML has occurred here. gray0.31	An error in the conversion from LaTeX to XML has occurred here. gray2.05
7.99	16.00	32.03	63.94	An error in the conversion from LaTeX to XML has occurred here. gray124.09	An error in the conversion from LaTeX to XML has occurred here. gray2.32	An error in the conversion from LaTeX to XML has occurred here. gray0.44	An error in the conversion from LaTeX to XML has occurred here. gray1.72	An error in the conversion from LaTeX to XML has occurred here. gray0.23
8.00	15.96	32.05	64.21	An error in the conversion from LaTeX to XML has occurred here. gray18.24	An error in the conversion from LaTeX to XML has occurred here. gray0.22	An error in the conversion from LaTeX to XML has occurred here. gray1.28	An error in the conversion from LaTeX to XML has occurred here. gray2.47	An error in the conversion from LaTeX to XML has occurred here. gray0.80

Table 8

Convergence rates obtained in solving the test-problem with β = 32 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The selected Fourth-order Four-stage Explicit Runge–Kutta Method (the most popular Runge–Kutta formula) and its combinations with eight versions of the Richardson Extrapolation are run. The perfect convergence rates are 16 , 32 , 64 , 128 , 256 , 512 , 1024 , 2048 and 4096.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.A.	–	–	–	–
N.S.	N.S.	N.S.	–	–	V.L.	1026.58	2060.13	4014.69
N.S.	–	–	128.85	259.06	1015.04	1023.81	2054.92	4508.38
–	32.06	63.69	127.61	256.74	1023.08	1020.83	2185.00	An error in the conversion from LaTeX to XML has occurred here. gray0.81
16.06	32.06	64.01	128.35	255.95	1023.62	1021.28	An error in the conversion from LaTeX to XML has occurred here. gray10.99	An error in the conversion from LaTeX to XML has occurred here. gray0.84
15.93	31.97	63.81	128.15	256.10	1032.52	An error in the conversion from LaTeX to XML has occurred here. gray155.80	An error in the conversion from LaTeX to XML has occurred here. gray 0.94	An error in the conversion from LaTeX to XML has occurred here. gray1.41
16.00	31.99	64.13	128.04	256.25	An error in the conversion from LaTeX to XML has occurred here. gray70.69	An error in the conversion from LaTeX to XML has occurred here. gray 0.77	An error in the conversion from LaTeX to XML has occurred here. gray1.03	An error in the conversion from LaTeX to XML has occurred here. gray1.11
16.00	32.00	64.11	127.93	255.49	An error in the conversion from LaTeX to XML has occurred here. gray0.73	An error in the conversion from LaTeX to XML has occurred here. gray1.09	An error in the conversion from LaTeX to XML has occurred here. gray0.95	An error in the conversion from LaTeX to XML has occurred here. gray0.48
15.97	31.98	64.37	127.91	246.80	An error in the conversion from LaTeX to XML has occurred here. gray0.95	An error in the conversion from LaTeX to XML has occurred here. gray1.21	An error in the conversion from LaTeX to XML has occurred here. gray0.71	An error in the conversion from LaTeX to XML has occurred here. gray1.41
15.97	32.03	64.05	128.18	An error in the conversion from LaTeX to XML has occurred here. gray8.25	An error in the conversion from LaTeX to XML has occurred here. gray0.72	An error in the conversion from LaTeX to XML has occurred here. gray0.84	An error in the conversion from LaTeX to XML has occurred here. gray0.76	An error in the conversion from LaTeX to XML has occurred here. gray0.92
16.06	31.96	63.73	127.24	An error in the conversion from LaTeX to XML has occurred here. gray0.77	An error in the conversion from LaTeX to XML has occurred here. gray1.54	An error in the conversion from LaTeX to XML has occurred here. gray0.51	An error in the conversion from LaTeX to XML has occurred here. gray0.86	An error in the conversion from LaTeX to XML has occurred here. gray0.69
16.00	32.06	64.33	An error in the conversion from LaTeX to XML has occurred here. gray44.01	An error in the conversion from LaTeX to XML has occurred here. gray1.43	An error in the conversion from LaTeX to XML has occurred here. gray0.47	An error in the conversion from LaTeX to XML has occurred here. gray1.33	An error in the conversion from LaTeX to XML has occurred here. gray0.64	An error in the conversion from LaTeX to XML has occurred here. gray0.46

7.4 Increasing the Accuracy and Stability Requirements

It is interesting to investigate what will happen if the accuracy requirements are increased, which can be done by simply increasing the value of β. We set β = 8192 instead of β = 32 and kept γ = - 750 in a second series of runs. In this way the oscillations are increased and, thus, the accuracy requirements are also increased. It is necessary to point out, however, that not only are the accuracy requirements increased. The absolute stability requirements are also increased because the absolute value of the complex eigenvalues of matrix A is becoming more than ten times larger than the absolute value of the real eigenvalue when this choice is made. Therefore, both stability and accuracy problems will occur in this case. Results for this choice of the two parameters are given in Tables 9–12 for the computational errors and in Tables 13–16 for the convergence rates. Four additional conclusions can also be drawn by studying the results given in Tables 9–16:

The results are becoming more and more accurate when q is increased. However, this problem cannot be solved by using the largest stepsizes even when the Fourth-order Four-stage ERKM is used in combination with the Seven-times Repeated Richardson Extrapolation. The stability is causing greater problems in the treatment of this test-example.
The expected convergence rates are nearly always achieved when the runs were successful (i.e., neither the accuracy nor the absolute stability is causing problems), but now the results are not as accurate as the results presented in Tables 5–8.
The underlying ERKMs are always unstable or not accurate (excepting the Fourth-order Four-stage ERKM when it is used with the smallest stepsize).
Also in these runs, the convergence rate actually achieved is sometimes greater than the expected convergence rate (see the results presented in Tables 14–16).

We use the following notation in Tables 9–16:

“N.S.” means that the method is not stable.
“N.A.” means that the error is greater than 10 - 2 .
“V.L.” means very large (which is caused presumably by the fact that ν i = h ⁢ μ i for i = 2 , 3 is very close to the boundary of the absolute stability region).
The rounding errors are dominant in the grey zone.

Table 9

Accuracy results obtained in solving the test-problem with β = 8192 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The First-order One-stage Explicit Runge–Kutta Method (the Explicit Euler Formula) and its combinations with eight versions of the Richardson Extrapolation are run with twelve stepsizes.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.A.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	1.0E - 02
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	9.8E - 03	1.6E - 04
N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	6.7E - 03	6.0E - 05	2.7E - 07
N.S.	N.S.	N.S.	N.S.	N.A.	5.7E - 03	5.0E - 05	2.4E - 07	5.1E - 10
N.S.	N.S.	N.S.	N.A.	9.9E - 03	9.0E - 05	3.8E - 07	9.2E - 10	9.3E - 13
N.S.	N.S.	N.A.	N.A.	3.0E - 04	1.4E - 06	2.9E - 09	3.6E - 12	1.9E - 15
N.S.	N.A.	N.A.	2.4E - 03	9.3E - 06	2.2E - 08	2.3E - 11	1.4E - 14	3.7E - 18

Table 10

Accuracy results obtained in solving the test-problem with β = 8192 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The Second-order Two-stage Explicit Runge–Kutta Method (the Improved Euler Formula) and its combinations with eight versions of the Richardson Extrapolation are run.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	N.A.
N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	N.A.	2.5E - 03	2.2E - 05
N.S.	N.S.	N.S.	N.S.	N.A.	N.A.	9.7E - 04	2.5E - 06	1.9E - 08
N.S.	N.S.	N.S.	N.A.	N.A.	8.7E - 04	3.8E - 06	2.8E - 09	1.9E - 11
N.S.	N.S.	N.S.	N.A.	7.2E - 04	6.5E - 06	1.5E - 08	2.9E - 12	1.9E - 14
N.S.	N.S.	N.A.	2.7E - 04	1.2E - 05	4.9E - 08	5.8E - 11	2.9E - 15	1.9E - 17
N.S.	N.A.	N.A.	4.3E - 06	1.8E - 07	3.7E - 10	2.3E - 13	2.8E - 18	1.9E - 20
N.A.	N.A.	6.8E - 04	6.7E - 08	2.9E - 09	2.9E - 12	8.9E - 16	2.8E - 21	1.9E - 23

Table 11

Accuracy results obtained in solving the test-problem with β = 8192 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The Third-order Three-stage Explicit Runge–Kutta Method (the Heun’s Formula) and its combinations with eight versions of the Richardson Extrapolation are run. “N.S.” means that the method is not stable. “N.A.” means that the error is greater than 10 - 2 .

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	N.A.	9.5E - 03	6.9E - 06
N.S.	N.S.	N.S.	N.S.	N.A.	N.A.	1.7E - 03	9.5E - 06	2.0E - 08
N.S.	N.S.	N.S.	N.S.	N.A.	6.8E - 04	4.9E - 06	1.0E - 08	1.5E - 11
N.S.	N.S.	N.A.	N.A.	3.0E - 04	2.3E - 06	9.9E - 09	1.2E - 11	7.8E - 15
N.S.	N.S.	N.A.	4.9E - 04	2.3E - 06	9.6E - 09	1.9E - 11	1.3E - 14	3.7E - 18
N.S.	N.A.	3.0E - 03	8.1E - 06	1.7E - 08	3.9E - 11	3.6E - 14	1.3E - 17	1.8E - 21
N.A.	N.A.	9.1E - 05	1.3E - 07	1.3E - 10	1.5E - 13	6.8E - 17	1.3E - 20	8.6E - 25
N.A.	1.8E - 03	2.8E - 06	2.0E - 09	1.0E - 12	6.1E - 16	1.3E - 19	1.2E - 23	4.2E - 28

Table 12

Accuracy results obtained in solving the test-problem with β = 8192 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The Fourth-order Four-stage Explicit Runge–Kutta Method (the most popular and the most used explicit Runge–Kutta formula) and its combinations with eight versions of the Richardson Extrapolation are run.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	2.1E - 03
N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	6.5E - 03	5.9E - 05	7.4E - 08
N.S.	N.S.	N.S.	N.S.	N.A.	2.8E - 03	1.2E - 05	3.0E - 08	3.2E - 11
N.S.	N.S.	N.S.	N.A.	7.9E - 04	3.2E - 06	1.1E - 08	1.9E - 11	1.1E - 14
N.S.	N.A.	N.A.	1.1E - 03	2.7E - 06	3.1E - 09	7.2E - 12	9.6E - 15	3.2E - 18
N.S.	N.A.	3.1E - 03	8.3E - 06	9.5E - 09	3.3E - 12	7.4E - 15	4.4E - 18	8.1E - 22
N.S.	N.A.	4.9E - 05	6.2E - 08	3.6E - 11	3.4E - 15	7.5E - 18	2.1E - 21	2.0E - 25
N.A.	4.3E - 04	7.8E - 07	4.8E - 10	1.4E - 13	3.3E - 18	7.4E - 21	1.0E - 24	5.1E - 29
6.3E - 03	1.3E - 05	1.2E - 08	3.7E - 12	5.5E - 16	3.3E - 21	7.3E - 24	7.3E - 28	An error in the conversion from LaTeX to XML has occurred here. gray2.8E - 30

Table 13

Convergence rates obtained in solving the test-problem with β = 8192 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The First-order One-stage Explicit Runge–Kutta Method (the Forward Euler Formula) and its combinations with eight versions of the Richardson Extrapolation are run. The perfect convergence rates should be 2 , 4 , 8 , 16 , 32 , 64 , 128 , 256 and 512.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.A.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	–
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	–	63.63
N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	–	164.21	572.99
N.S.	N.S.	N.S.	N.S.	N.A.	–	134.33	254.47	535.16
N.S.	N.S.	N.S.	N.A.	–	63.03	131.84	254.33	551.72
N.S.	N.S.	N.A.	N.A.	32.93	63.69	129.69	255.25	488.42
N.S.	N.A.	N.A.	–	32.33	63.51	129.07	254.93	513.51

Table 14

Convergence rates obtained in solving the test-problem with β = 8192 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The Second-order Two-stage Explicit Runge–Kutta Method (the Forward Euler Formula) and its combinations with eight versions of the Richardson Extrapolation are run. The perfect convergence rates should be 4 , 8 , 16 , 32 , 64 , 128 , 256 , 512 and 1024.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	N.A.
N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	N.A.	–	–
N.S.	N.S.	N.S.	N.S.	N.A.	N.A.	–	1000.00	1235.63
N.S.	N.S.	N.S.	N.A.	N.A.	–	256.23	900.36	920.63
N.S.	N.S.	N.S.	N.A.	–	137.73	254.73	975.69	979.27
N.S.	N.S.	N.A.	–	62.17	132.17	254.30	1006.99	1005.21
N.S.	N.A.	N.A.	63.14	63.19	130.48	256.26	1014.18	1015.87
N.A.	N.A.	–	63.58	63.64	129.41	256.05	1021.74	1016.13

Table 15

Convergence rates obtained in solving the test-problem with β = 8192 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The Third-order Three-stage Explicit Runge–Kutta Method (the Heun’s Formula) and its combinations with eight versions of the Richardson Extrapolation are run. The perfect convergence rates should be 8 , 16 , 32 , 64 , 128 , 256 , 512 , 1024 and 2048.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	N.A.	–	–
N.S.	N.S.	N.S.	N.S.	N.A.	N.A.	–	382.11	352.55
N.S.	N.S.	N.S.	N.S.	N.A.	–	345.68	913.46	1289.47
N.S.	N.S.	N.A.	N.A.	–	302.65	490.41	845.53	1943.73
N.S.	N.S.	N.A.	–	130.43	236.40	527.13	960.94	2090.91
N.S.	N.A.	–	60.20	132.95	247.69	529.58	1000.00	2101.12
N.A.	N.A.	32.57	62.71	132.06	252.29	522.83	1007.87	2077.01
N.A.	–	32.80	63.55	131.00	251.64	518.32	1024.19	2060.10

Table 16

Convergence rates obtained in solving the test-problem with β = 8192 and γ = - 750 for twelve successive runs, starting with a stepsize h = 0.02048 and reducing the stepsize by a factor of two after every run. The Fourth-order Four-stage Explicit Runge–Kutta Method (the most popular and the most used explicit Runge–Kutta formula) and its combinations with eight versions of the Richardson Extrapolation are run. The perfect convergence rates should be 16 , 32 , 64 , 128 , 256 , 512 , 1024 , 2048 and 4096.

Direct	CRE	Rep.	2TRRE	3TRRE	4TRRE	5TRRE	6TRRE	7TRRE
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.
N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	–
N.S.	N.S.	N.S.	N.S.	N.S.	N.A.	–	–	V.L.
N.S.	N.S.	N.S.	N.S.	N.A.	–	541.67	1983.22	2338.61
N.S.	N.S.	N.S.	N.A.	–	8703.70	1132.10	1560.21	2846.85
N.S.	N.A.	N.A.	–	295.15	1035.14	1476.32	2000.00	3512.66
N.S.	N.A.	–	V.L.	281.52	937.13	967.65	2279.24	3886.84
N.S.	N.A.	62.40	132.21	263.71	988.17	994.64	2127.96	4004.93
N.A.	–	63.24	120.27	257.86	1011.98	1009.47	2089.11	3957.11
–	32.50	63.77	129.11	255.94	1015.19	1016.51	1387.36	18.19

8 Major Achievements and Plans for Future Research

The basic properties and the performance of eight versions of the Richardson Extrapolation were studied in the previous sections. The following important topics were mainly discussed.

8.1 Accuracy Properties of the Different Versions of the Richardson Extrapolation

It was shown (Theorem 4.1) that the order of accuracy of the combinations of eight versions of the Richardson Extrapolation with different one-step methods for solving systems of ODEs d ⁢ y d ⁢ t = f ⁢ ( t , y ) is always higher than that of the underlying method when the right-hand-side vector f is sufficiently smooth.

8.2 Stability Functions

It was proved (Theorem 5.1) that if we know the stability function of the selected one-step method for solving systems of ODEs that are represented by formula (2.1), then it will be possible to construct the stability functions of any of the studied versions of the Richardson Extrapolation.

8.3 Absolute Stability Regions Related to the ERKMs

If ERKMs with m = p and m = 1 , 2 , 3 , 4 are used together with some version of the Richardson Extrapolation, then it was shown that the combined numerical methods have always better stability properties (i.e., larger absolute stability regions) than the underlying ERKM (see Figure 1 and Theorem 6.2).

8.4 Development of a Two-Parameter Family of Test-Examples

A two-parameter family, by which both the accuracy and the stability properties of the numerical methods can easily be controlled by varying in an appropriate way the parameters, was developed and a long series of experiments was carried out. All eight versions were very systematically tested in the cases where the underlying methods were Explicit Runge–Kutta Methods (ERKMs).

8.5 Use of More Complicated Examples

We are planning to apply advanced versions of the Richardson Extrapolation in the numerical treatment of an atmospheric chemical scheme. This scheme is described by a non-linear system of ODEs, which is very stiff, badly scaled and extremely ill-conditioned (it was found, by using subroutines from [3, 28, 38] that the condition numbers of the involved Jacobian matrices of the non-linear system of ODEs varied in the interval [ 4.56 ⋅ 10 8 , 9.27 ⋅ 10 12 ] . The atmospheric chemical scheme is used in the Unified Danish Eulerian Model, which has been used in many studies related to the pollution levels and to the increases of some dangerous pollution levels as a results of the climatic changes both in Europe as a whole ([4, 12, 29, 30, 31]) and in different parts of Europe, such as the UK ([1]), the Balkan Peninsula ([39], Bulgaria ([35, 43, 44]), Denmark ([42]), Hungary ([17, 41]), several countries in Eastern and Central Europe [8] and the North Sea [15]. The Unified Danish Eulerian Model and other large-scale European air pollution models were also studied in [16, 24].

8.6 Use of Implicit Runge–Kutta Methods

The Implicit Runge–Kutta Methods deserve some special treatment, because it is not easy to study the absolute stability properties of the numerical methods resulting after the application of the advanced versions of the Richardson Extrapolation together with IRKMs. In fact, that is extremely difficult. The following example explains the difficulties. The investigation of the Seven-times Repeated Richardson Extrapolations together with the simplest IRKMs, the first-order Backward Differentiation Formula, implies study of stability functions, which contain polynomials of degree 256 both in the numerator and in the denominator. Most of the difficulties are now successfully resolved and we are preparing a paper about the use of IRKMs together with the Advanced Versions of the Richardson Extrapolation, which will focus on the absolute stability properties of the resulting methods.

8.7 Linear Multistep Methods

The Richardson Extrapolation is not suitable for application together with linear multistep methods. Assume that a linear k-step method is to be used in the calculation of y n . Then approximations y n - 1 , y n - 2 , … , y n - k of the unknown function are needed at the points x n - 1 = x n - h , x n - 2 = x n - 2 ⁢ h , … , x n - k = x n - k ⁢ h when the linear k-step method is used. The application of the Classical Richardson Extrapolation will require (when computations with a stepsize 0.5 ⁢ h are to be carried out) additional approximations y n - 1 2 , y n - 3 2 , … , y n + 1 2 - k at the points x n - 1 2 = x n - 0.5 ⁢ h , x n - 3 2 = x n - 1.5 ⁢ h , … , x n - k + 1 2 = x n - ( 0.5 ⁢ h ) ⁢ k . This means that many extra evaluations of the unknown function are needed. The situation becomes much more complicated when repeated versions of the Richardson Extrapolation are to be used together with linear k-step methods. However, it should be mentioned that it is possible to use another efficient approach, predictor-corrector schemes, when linear multistep methods are to be handled. This possibility is discussed in Chapter 3 of [34].

8.8 Efficient Implementation

We have not discussed the problems related to parallel implementations of the Richardson Extrapolation in computer codes, because the paper is rather long and we left this issue for a future investigation. If the problem solved is large, then it will be necessary to implement the algorithms in an efficient way and/or to carry out the most time-consuming operations in parallel. Note that the most time-consuming operations used when Explicit Runge–Kutta Methods are applied are Basic Linear Algebra Operations (mainly matrix-vector multiplications). Well-developed kernels for performing Basic Linear Algebra Operations in parallel are available on many computers, see, for example, [3] and can easily be called in order to achieve parallelism and to improve the performance.

Furthermore, the implementation of Explicit Runge–Kutta Methods should be made very carefully when large-scale problems are to be handled. For example, k 1 n has to be used eight times when the Seven-time Repeated Richardson Extrapolation is selected. It will be appropriate to develop a code where this quantity is calculated only once and used every times when this is needed. Some other similar approaches are to be used during the development of an efficient code for the treatment of the different versions of the Richardson Extrapolation described in this paper.

8.9 Plans for Future Research

Accuracy and absolute stability results related to the application of eight versions of the Richardson Extrapolation used together with Explicit Runge–Kutta methods were established and verified in this paper by using appropriate sets of test-examples. The correctness of the proved results was illustrated for the case where explicit numerical methods were used. The next step will be the use of some implicit methods (the paper being in preparation). Furthermore, it will be useful, as mentioned above, to develop a reliable and efficient software, which will allow the users to solve their problems by automatically changing both the stepsize and the used versions of the Richardson Extrapolation according to the accuracy and absolute stability requirements of the solved by them scientific and engineering problems.

8.10 Conjectures

Several conjectures were formulated in this paper. Some experiments with First-order One-stage Explicit Runge–Kutta Method in combination with the Eight-times Repeated Richardson Extrapolation have been carried out. We found that two conditions are satisfied when the Eight-times Repeated Richardson Extrapolation is used:

the coefficients of this version can be calculated by using the generic formula (4.24),
the derived in this way version of the Richardson Extrapolation behaves as a numerical method the order of accuracy of which is ten (as should be expected according to Theorem 4.1 and Conjecture 1).

This indicates that the statement made in this conjecture is true, but an attempt to prove strictly all conjectures made above will be useful and will be made.

Dedicated to Alexander Andreevich Samarskii. We are expressing our great respect to his numerous valuable contributions in different fields of applied mathematics. He greatly influenced the foundation of numerical analysis as well as the development of the applied mathematical schools in Bulgaria and Hungary as well as in many other countries.

Funding statement: This research was partially carried out in the ELTE Institutional Excellence Program (1783-3/2018/FEKUTSRAT) supported by the Hungarian Ministry of Human Capacities, and supported by the Hungarian Scientific Research Fund OTKA, No. K112157 and No. SNN125119, and National Scientific Program ICT-SDM-SES. The research reported in this paper has been supported by the National Research, Development and Innovation Fund (TUDFO/51757/2019-ITM, Thematic Excellence Program).

References

[1] S. Abdalmogith, R. M. Harrison and Z. Zlatev, Intercomparison of inorganic aerosol concentrations in the UK with predictions of the Danish Eulerian Model, J. Atmospheric Chem. 54 (2004), 43–66. 10.1007/s10874-006-9012-3Search in Google Scholar

[2] R. Alexander, Diagonally implicit Runge–Kutta methods for stiff O.D.E.s, SIAM J. Numer. Anal. 14 (1977), no. 6, 1006–1021. 10.1137/0714068Search in Google Scholar

[3] V. A. Barker, L. S. Blackford, J. Dongarra, J. Du Croz, S. Hammarling, M. Marinova, J. Waśniewski and P. Yalamov, LAPACK95 Users’ Guide, Software Environ. Tools 13, SIAM, Philadelphia, 2001. 10.1137/1.9780898718201Search in Google Scholar

[4] A. Bastrup-Birk, J. Brandt, I. Uria and Z. Zlatev, Studying cumulative ozone exposures in Europe during a 7-year period, J. Geophys. Res. 102 (1997), 23917–23035. 10.1029/97JD01966Search in Google Scholar

[5] K. Burrage, Parallel and Sequential Methods for Ordinary Differential Equations, Numer. Math. Sci. Comput., Oxford University, New York, 1995. 10.1093/oso/9780198534327.001.0001Search in Google Scholar

[6] J. C. Butcher, Numerical Methods for Ordinary Differential Equations, 2nd ed., John Wiley & Sons, Chichester, 2003. 10.1002/0470868279Search in Google Scholar

[7] E. Christiansen and H. G. Petersen, Estimation of convergence orders in repeated Richardson extrapolation, BIT 29 (1989), no. 1, 48–59. 10.1007/BF01932705Search in Google Scholar

[8] P. Csomós, R. Cuciureanu, G. Dimitriu, I. Dimov, A. Doroshenko, I. Faragó, K. Georgiev, Á. Havasi, R. Horváth, S. Margenov, T. Ostromsky, V. Prusov, D. Syrakov and Z. Zlatev, Impact of climate changes on pollution levels in Europe, Final report for a NATO Linkage Project (Grant 980505), 2006, http://www.cs.elte.hu/~faragois/NATO.pdf. Search in Google Scholar

[9] G. G. Dahlquist, A special stability problem for linear multistep methods, Nordisk Tidskr. Inform.-Behand. 3 (1963), 27–43. 10.1007/BF01963532Search in Google Scholar

[10] I. Faragó, A. Havasi and Z. Zlatev, Efficient implementation of stable Richardson extrapolation algorithms, Comput. Math. Appl. 60 (2010), no. 8, 2309–2325. 10.1016/j.camwa.2010.08.025Search in Google Scholar

[11] C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, 1971. Search in Google Scholar

[12] G. Geernaert and Z. Zlatev, Studying the influence of the biogenic emissions on the AOT40 levels in Europe, Int. J. Environ. Pollution 23 (2004), no. 12, 29–41. 10.1504/IJEP.2004.005485Search in Google Scholar

[13] E. Hairer, S. P. Nørsett and G. Wanner, Solving Ordinary Differential Equations: I. Nonstiff Problems, Springer, Berlin, 1987. 10.1007/978-3-662-12607-3Search in Google Scholar

[14] E. Hairer and G. Wanner, Solving Ordinary Differential Equations: II. Stiff and Differential-Algebraic Problems, Springer, Berlin, 1991. 10.1007/978-3-662-09947-6Search in Google Scholar

[15] R. M. Harrison, Z. Zlatev and C. J. Ottley, A comparison of the predictions of a Eulerian atmospheric transport chemistry model with experimental measurements over the North Sea, Atmospheric Environ. 28 (1994), 497–516. 10.1016/1352-2310(94)90127-9Search in Google Scholar

[16] H. Hass, M. van Loon, C. Kessler, R. Stern, J. Mathijsen, F. Sauter, Z. Zlatev, J. Langner, V. Foltescu and M. Schaap, Aerosol modelling: Results and intercomparison from European regional-scale modelling systems, GSF–National Research Center for Environment and Health, International Scientific Secretariat (ISS), EUROTRAC-2, Munich, 2004, http://www.trumf.fu-berlin.de/veranstaltungen/events/gloream, http://www.rivm.nl/bibliotheek/digitaaldepot/GLOREAM_PMmodel-comparison.pdf. Search in Google Scholar

[17] Á. Havasi and Z. Zlatev, Trends of Hungarian air pollution levels on a long time-scale, Atmospheric Environ. 36 (2002), 4145–4156. 10.1016/S1352-2310(02)00198-XSearch in Google Scholar

[18] P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, John Wiley & Sons, New York, 1962. Search in Google Scholar

[19] K. Heun, Neue Methode zur approximativen Interation der Differentialgleichungen einer unabhängigen Veränderlichen, Z. Math. Phys. 45 (1900), 23–38. Search in Google Scholar

[20] W. Hundsdorfer and J. Verwer, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, Springer Ser. Comput. Math. 33, Springer, Berlin, 2003. 10.1007/978-3-662-09017-6Search in Google Scholar

[21] W. Kutta, Beitrag zur näherungsweisen Integration totaler Differentialgleichungen, Z. Math. Phys. 46 (1901), 435–453. Search in Google Scholar

[22] J. D. Lambert, Numerical Methods for Ordinary Differential Systems. The Initial Value Problem, John Wiley & Sons, Chichester, 1991. Search in Google Scholar

[23] L. F. Richardson, The deferred approach to the limit, I-single lattice, Philos. Trans. Roy. Soc. Lond. Ser. A 226 (1927), 299–349. 10.1098/rsta.1927.0008Search in Google Scholar

[24] M. Roemer, M. Beekman, R. Bergsröm, G. Boersen, H. Feldmann, F. Flatoy, C. Honore, J. Langner, J. E. Jonson, J. Matthijsen, M. Memmesheimer, D. Simpson, P. Smeets, S. Solberg, D. Stevenson, P. Zandveld and Z. Zlatev, Ozone trends according to ten dispersion models, GSF – National Research Center for Environment and Health, International Scientific Secretariat (ISS), EUROTRAC-2, Munich, 2004, http://www.mep.tno.nl/eurotrac/EUROTRAC-trends.pdf. Search in Google Scholar

[25] C. Runge, Ueber die numerische Auflösung von Differentialgleichungen, Math. Ann. 46 (1895), no. 2, 167–178. 10.1007/BF01446807Search in Google Scholar

[26] L. F. Shampine, Numerical Solution of Ordinary Differential Equations, Chapman & Hall, New York, 1994. Search in Google Scholar

[27] L. F. Shampine and M. K. Gordon, Computer Solution of Ordinary Differential Equations. The Initial Value Problem, W. H. Freeman, San Francisco, 1975. Search in Google Scholar

[28] Z. Zlatev, Computational Methods for General Sparse Matrices, Math. Appl. 65, Kluwer Academic, Dordrecht, 1991. 10.1007/978-94-017-1116-6Search in Google Scholar

[29] Z. Zlatev, Computer Treatment of Large Air Pollution Models, Kluwer Academic, Dordrecht, 1995. 10.1007/978-94-011-0311-4Search in Google Scholar

[30] Z. Zlatev, Impact of future climate changes on high ozone levels in European suburban areas, Climatic Change 101 (2010), 447–483. 10.1007/s10584-009-9699-7Search in Google Scholar

[31] Z. Zlatev and I. Dimov, Computational and Numerical Challenges in Environmental Modelling, Elsevier, Amsterdam, 2006. Search in Google Scholar

[32] Z. Zlatev, I. Dimov, I. Faragó, K. Georgiev and Á. Havasi, Stability of the Richardson extrapolation combined with some implicit Runge–Kutta methods, J. Comput. Appl. Math. 310 (2017), 224–240. 10.1016/j.cam.2016.03.018Search in Google Scholar

[33] Z. Zlatev, I. Dimov, I. Faragó, K. Georgiev and Á. Havasi, Stability properties of the Two-times Repeated Richardson Extrapolation combined with some Explicit Runge–Kutta Methods, Talk presented at the conference in Lozenets, June 2018. 10.1007/978-3-030-11539-5_80Search in Google Scholar

[34] Z. Zlatev, I. Dimov, I. Faragó and Á. Havasi, Richardson Extrapolation: Practical Aspects and Applications, De Gruyter Ser. Appl. Numer. Math. 2, De Gruyter, Berlin, 2018. 10.1515/9783110533002Search in Google Scholar

[35] Z. Zlatev, I. Dimov and K. Georgiev, Relations between climatic changes and high pollution levels in Bulgaria, Open J. Appl. Sci. 6 (2016), 386–401. 10.4236/ojapps.2016.67040Search in Google Scholar

[36] Z. Zlatev, I. Faragó and Á. Havasi, Stability of the Richardson extrapolation applied together with the θ-method, J. Comput. Appl. Math. 235 (2010), no. 2, 507–517. 10.1016/j.cam.2010.05.052Search in Google Scholar

[37] Z. Zlatev, I. Faragó and Á. Havasi, Richardson extrapolation combined with the sequential splitting procedure and the θ-method, Cent. Eur. J. Math. 10 (2012), no. 1, 159–172. 10.2478/s11533-011-0099-7Search in Google Scholar

[38] Z. Zlatev and K. Georgiev, Applying approximate LU-factorizations as preconditioners in eight iterative methods for solving systems of linear algebraic equations, Cent. Eur. J. Math. 11 (2013), no. 8, 1510–1530. 10.2478/s11533-013-0248-2Search in Google Scholar

[39] Z. Zlatev, K. Georgiev and I. Dimov, Influence of climatic changes on pollution levels in the Balkan Peninsula, Comput. Math. Appl. 65 (2013), no. 3, 544–562. 10.1016/j.camwa.2012.07.006Search in Google Scholar

[40] Z. Zlatev, K. Georgiev and I. Dimov, Studying absolute stability properties of the Richardson extrapolation combined with explicit Runge–Kutta methods, Comput. Math. Appl. 67 (2014), no. 12, 2294–2307. 10.1016/j.camwa.2014.02.025Search in Google Scholar

[41] Z. Zlatev, Á. Havasi and I. Faragó, Influence of climatic changes on pollution levels in Hungary and its surrounding countries, Atmosphere 2 (2011), 201–221. 10.3390/atmos2030201Search in Google Scholar

[42] Z. Zlatev and L. Moseholm, Impact of climate changes on pollution levels in Denmark, Ecological Model. 217 (2008), 305–319. 10.1016/j.ecolmodel.2008.06.030Search in Google Scholar

[43] Z. Zlatev and D. Syrakov, A fine resolution modelling study of pollution levels in Bulgaria. Part 1: SOx and NOx pollution, Int. J. Environment Pollution 22 (2004), no. 1–2, 186–202. 10.1504/IJEP.2004.005508Search in Google Scholar

[44] Z. Zlatev and D. Syrakov, A fine resolution modelling study of pollution levels in Bulgaria. Part 2: High ozone levels, Int. J. Environment Pollution 22 (2004), no. 1–2, 203–222. 10.1504/IJEP.2004.005513Search in Google Scholar

Received: 2019-01-27

Revised: 2019-05-20

Accepted: 2019-06-17

Published Online: 2019-07-16

Published in Print: 2020-10-01

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/cmam-2019-0016

Keywords for this article

Several Times Repeated Richardson Extrapolation; Order of Accuracy; Absolute Stability Properties; Systems of ODEs; Runge–Kutta Methods; Numerical Examples

Creative Commons

BY 4.0