On the Consistency and Convergence of Classical Richardson Extrapolation as Applied to Explicit One-Step Methods

. The consistency of the classical Richardson extrapolation (CRE), a simple and robust computational device, is analysed for the case where the underlying method is an explicit one-step numerical method for ordinary differential equations with order of consistency one or two. It is shown in the classical framework that the CRE increases the order of consistency by one. The convergence of the method is proved by the assumption that the time-stepping operator of the base method has the Lipschitz property in its second argument

In certain model problems the accuracy of the numerical solution obtained by a chosen numerical method is not sufficient. In this case we can either decrease the time-step size, and perform the integration with the smaller timestep, or choose a higher-order numerical method, which is, as a rule, more costly. In both cases our previous calculations are thrown away. Another approach is the classical Richardson extrapolation (CRE) [11,12,14,15,16], where a suitable linear combination of numerical solutions obtained by two different time-step sizes with the same base method (BM) are calculated. An attractive property of this procedure is that its implementation requires relatively simple modifications in the computer code of a BM whose code is available.
Numerous experiments have shown that the CRE increases the accuracy of the applied BM by one [15], and for special (symmetric) numerical methods [5,6], where the global error has an asymptotic expansion in even powers of the time-step, the order increases by two. In these experiments presented in [15], the maximal error committed during the time integration or the error committed at a given time point is investigated as the time step size is decreased, which serves as an experimental confirmation that the numerical solution is convergent to the exact solution. Convergence in this sense is a major requirement for a numerical time integration method, and the higher the convergence order, the smaller time step can be used for achieving a prescribed accuracy. According to the classical theory of numerical analysis, when a onestep numerical method is applied, the Lipschitz property of the time-stepping operator implies convergence in case of consistency, and the order of the convergence is equal to the order of consistency. Therefore, the knowledge of the consistency and its order is of crucial importance to study the convergence. Obviously, a higher order consistency is favourable, especially when accuracy requirements are higher.
In [7] the order of consistency is shown to increase to p + 1 when a p-th order Runge-Kutta base method is combined with the CRE, by considering the equivalent autonomous problem. In this paper we use a different approach. We do not consider the equivalent autonomous problem, but a non-autonomous problem directly. Moreover, any explicit one-step method can serve as a BM in our case. Our analysis is based on a classical Taylor series expansion. We restrict ourselves to the cases where the BM has either order one or order two, since in most engineering applications second or third order of consistency is required.
The structure of the paper is as follows. In Section 2, the conditions of order p of a general explicit one-step method are derived for a non-autonomous Cauchy problem. In Section 3, we prove that the order of consistency is increased by one when the chosen BM of consistency order p = 1 or p = 2 is combined with the CRE. Then, in Section 4, the Lipschitz property of the time-stepping function of any explicit one-step BM of consistency order p = 1 or p = 2 combined with the CRE is proven, provided the BM satisfies a Lipschitz property. This implies that in this case the combined scheme is convergent, and the order of convergence is equal to the order of consistency. Finally, our theoretical results are illustrated with numerical experiments.

Consistency of explicit one-step methods
Consider the Cauchy problem where f : R 2 → R (scalar case). Later on we always assume that this problem has a unique, sufficiently smooth solution u(t). Assume that (2.1) is to be solved by the general explicit one-step method y n = y n−1 + hΦ(t n−1 , y n−1 , h), where Φ : R 3 → R is a smooth increment function. As known, the truncation error of this method [13] is defined as and the method (2.2) is said to be p-th order consistent when T n (h) = O(h p ). First, we examine the conditions of first-order consistency. Applying Taylor series expansion for the exact solution around t n−1 and for the function Φ around zero, we get Hence, T n (h) = O(h) if and only if u ′ (t n−1 ) = Φ(t n−1 , u(t n−1 ), 0). In view of the differential equation, this implies the requirement Φ(t n−1 , u(t n−1 ), 0) = f (t n−1 , u(t n−1 )). Hence, the first order consistency is achieved if and only if the condition holds for any (t, u) ∈ R + × R. This result can be found in [9], p. 152, however, the order is not investigated there. However, in a similar way as above, the necessary and sufficient conditions of p-th order consistency can be given as Particularly, the method is second order consistent if, in addition to (2.3), the condition holds.
Remark 1. Conditions of p-th order consistency are derived in [1], but for autonomous systems, only. Namely, for the autonomous equation Moreover, the method (2.6) is consistent with (2.5) at order p, if F ∈ C p and 1] are functions of the type F : R d → R d defined recursively, and the formulated conditions are rather complicated to check. Our approach is for the scalar non-autonomous case, which results in a condition that is simple and easy to check.
3) trivially holds, while (2.4) does not, so this method has first order.

Example 2. For the modified Euler method (MEE) the increment function reads
3) trivially holds, and so does (2.4), since and the substitution h = 0 directly yields (2.4). Therefore, the method has (at least) second order.

Order conditions for the CRE
Assume that some BM is applied to solve the Cauchy problem (2.1). We divide [0, T ] into N sub-intervals of length h, and define two meshes on [0, T ]: Note that Ω 0 h ⊂ Ω 1 h . We will only deal with the active Richardson extrapolation, which is defined as follows. Assume that the numerical solution has been calculated at the time level t n−1 of the coarse mesh Ω 0 h , and denote it by y n−1 . Take a step of length h by using the BM, the resulting numerical solution will be denoted as z n .
Take also two steps of step size h/2 with the BM from y n−1 , and denote the obtained numerical solution by w n (belonging also to time t n of the mesh Ω 0 h ). The CRE method is defined as Case p = 1: Assume that we have chosen an explicit one-step base method, i.e., a BM of the form (2.2) with time-stepping function Φ b , which is first order consistent, i.e., (2.3) holds. We show that its combination with CRE yields a second-order consistent method. To this aim, we first derive the increment function Φ b,CRE of the combined method.
2. Take two steps of length h/2: by using the notation t n− 1 2 := t n−1 + h 2 we can write Combine the two numerical solutions: y n = 2w n − z n . So, according to (2.1) we have . Hence, the three-variable function determining the combined method reads as  4). Since the base method Φ b has first order, therefore (2.3) holds due to the relation For the condition (2.4) we have Taking this function at h = 0 yields Assume that the base method has order 2. Then the numerical solution obtained by CRE reads Therefore, the function of the combined method is Assuming that the base method satisfies conditions (2.3)-(2.4), we need to show that the following three conditions are satisfied: (3.8) In the sequel, when no argument is written for f , the functions are to be understood at (t, u). Checking condition (3.6): Checking condition (3.7): Substituting h = 0 yields Checking condition (3.8): a simple, but tedious calculation shows that the function ∂ 33 Φ b,CRE (t, u, h) at h = 0 is equal to , and always assuming sufficient smoothness of Φ b and f for the interchangeability of the differentiation with respect to any two variables, we obtain that Here, we used the equalities

Convergence analysis
We refer the following theorem from [13].
Theorem 1. Consider a one-step method (2.2) with order of consistency p, where Φ satisfies the Lipschitz condition with respect to its second argument uniformly in t, i.e., there exist constants L Φ and h 0 such that for 0 ≤ h ≤ h 0 Let e n := u(t n ) − y n be the global error at the n-th time step. Then, in case e 0 = O(h p ), we have that |e n | ≤ O(h p ), so the method is convergent of order p.
Our aim is to show the Lipschitz property for the time-stepping operator of the combined method. We consider the cases p = 1 and p = 2 separately. In both cases we will assume that the BM satisfies the Lipschitz condition We remark that this has been shown for the ERK methods if the right-hand side function f satisfies the Lipschitz property in its second argument [4].

The Lipschitz property of Φ b,CRE for p = 1
Let us show now that Φ b,CRE is Lipschitz in its second variable if Φ b satisfies (4.1) with Lipschitz constant L b,1 . We need to show that there exists L CRE,1 > 0 constant such that where T is the length of the time interval of the problem. So, (4.2) holds by the Lipschitz constant L CRE,1 := 6 + T 2 L b,1 .

The Lipschitz property of Φ b,CRE for p = 2
We show that by the assumption (4.1) with Lipschitz constant L b,2 there exists L CRE,2 > 0 constant s.t. for the increment function (3.5) A straightforward calculation shows that

Numerical experiment
We considered the scalar problem The exact solution is y(t) = 2 cot −1 e t 2 cot( 1 2 ) . We solved the problem with the first order EE method (Table 1) and the second order MEE method ( Table 2) as well as with their combinations with the CRE. For comparison we included the second-order Trapezoidal (TR) method and the third-order Heun's method (Table 3). We also compared our results with an implicit method, the implicit Trapezoidal (ITR) method, which is known to be a symmetric scheme [5]. The global errors were calculated in absolute value at the end of the time interval [0, 1] for time step sizes, successively  reduced by a factor of 2. The columns under "Order" show estimations of the convergence order by the formula p ≈ ln(E(h)) − ln(E(h/2)) ln 2 , where E(h) stands for the global error obtained with time-step size h. As expected, the EE method behaves as a first order, while the MEE, TR, and ITR as second order convergent methods. After the application of CRE, we get estimated orders around 2 for the EE method, 3 for the MEE and TR methods, and 4 for the ITR method. So the order of convergence was successfully increased by one with the method of classical Richardson extrapolation, and the order of convergence is increased by two when the BM is symmetric. Note that the EE + CRE method outperforms the TR method for all the studied time steps, but Heun's method is more accurate that MEE + CRE for each of the studied the time step lengths. We remark however that as far as efficiency is considered, CRE allows parallelization, which can significantly reduce the computational time.

Conclusions
We gave a consistency analysis for any explicit one-step numerical method of order p = 1 or p = 2, combined with the classical Richardson extrapolation when applied to a non-autonomous Cauchy problem. We showed that the CRE increases the order of consistency by one. Since the Lipschitz property of the time-stepping function of the method is inherited by the time-stepping function of the combined (BM + CRE) method, the same order of convergence is guaranteed if the base method possesses the Lipschitz property.