About this Journal Submit a Manuscript Table of Contents 
International Journal of Differential Equations
Volume 2011 (2011), Article ID 520840, 16 pages
doi:10.1155/2011/520840
Research Article

Solving Advection Equations by Applying the Crank-Nicolson Scheme Combined with the Richardson Extrapolation

Zahari Zlatev,1 Ivan Dimov,2 István Faragó,3 Krassimir Georgiev,2 Ágnes Havasi,4 and Tzvetan Ostromsky2

1National Environmental Research Institute, Aarhus University, DK-4000 Roskilde, Denmark
2Institute of Information and Communication Technologies, Bulgarian Academy of Sciences, BG-1113 Sofia, Bulgaria
3Department of Applied Analysis and Computational Mathematics, Eötvös Loránd University, H-1117 Budapest, Hungary
4Department of Meteorology, Eötvös Loránd University, H-1117 Budapest, Hungary

Received 23 May 2011; Revised 3 August 2011; Accepted 15 September 2011

Academic Editor: Roderick Melnik

Copyright © 2011 Zahari Zlatev et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Advection equations appear often in large-scale mathematical models arising in many fields of science and engineering. The Crank-Nicolson scheme can successfully be used in the numerical treatment of such equations. The accuracy of the numerical solution can sometimes be increased substantially by applying the Richardson Extrapolation. Two theorems related to the accuracy of the calculations will be formulated and proved in this paper. The usefulness of the combination consisting of the Crank-Nicolson scheme and the Richardson Extrapolation will be illustrated by numerical examples.

1. Statement of the Problem

Consider the following advection equation: 𝜕 𝑐 𝜕 𝑡 = − 𝑢 𝜕 𝑐 𝑎 𝜕 𝑥 , 𝑥 ∈ 1 , 𝑏 1 [ ] ⊂ ( − ∞ , ∞ ) , 𝑡 ∈ 𝑎 , 𝑏 ⊂ ( − ∞ , ∞ ) . ( 1 . 1 ) It will be assumed that 𝑢 = 𝑢 ( 𝑥 , 𝑡 ) is a given function. The physical meaning of this function depends on the area in which (1.1) arises. For example, in fluid dynamics applications it can be interpreted as fluid velocity, while in atmospheric modelling 𝑢 ( 𝑥 , 𝑡 ) represents the wind velocity.

Equation (1.1) must always be considered together with appropriate initial and boundary conditions.

The well-known Crank-Nicolson scheme (see, e.g., [1, page 63]) can be used in the numerical treatment of (1.1). The computations are carried out by the following formula: 𝜎 𝑖 , 𝑛 + 0 . 5 𝑐 𝑖 + 1 , 𝑛 + 1 + 𝑐 𝑖 , 𝑛 + 1 − 𝜎 𝑖 , 𝑛 + 0 . 5 𝑐 𝑖 − 1 , 𝑛 + 1 + 𝜎 𝑖 , 𝑛 + 0 . 5 𝑐 𝑖 + 1 , 𝑛 − 𝑐 𝑖 , 𝑛 − 𝜎 𝑖 , 𝑛 + 0 . 5 𝑐 𝑖 − 1 , 𝑛 = 0 , ( 1 . 2 ) when the Crank-Nicolson scheme is used. The quantity 𝜎 𝑖 , 𝑛 + 0 . 5 is defined by 𝜎 𝑖 , 𝑛 + 0 . 5 = 𝑘 𝑢 𝑥 4 ℎ 𝑖 , 𝑡 𝑛 + 0 . 5 , ( 1 . 3 ) where 𝑡 𝑛 + 0 . 5 𝑘 = 𝑡 𝑛 + 0 . 5 𝑘 and the increments ℎ and 𝑘 of the spatial and time variables are introduced by using two equidistant grids: 𝐺 𝑥 = 𝑥 𝑖 , 𝑖 = 0 , 1 , … , 𝑁 𝑥 ∣ 𝑥 0 = 𝑎 1 , 𝑥 𝑖 = 𝑥 𝑖 − 1 + ℎ , 𝑖 = 1 , 2 , … , 𝑁 𝑥 𝑏 , ℎ = 1 − 𝑎 1 𝑁 𝑥 , 𝑥 𝑁 𝑥 = 𝑏 1 𝐺 , ( 1 . 4 ) 𝑡 = 𝑡 𝑛 , 𝑛 = 0 , 1 , … , 𝑁 𝑡 ∣ 𝑡 0 = 𝑎 , 𝑡 𝑛 = 𝑡 𝑛 − 1 + 𝑘 , 𝑛 = 1 , 2 , … , 𝑁 𝑡 , 𝑘 = 𝑏 − 𝑎 𝑁 𝑡 , 𝑡 𝑁 𝑡 = 𝑏 . ( 1 . 5 ) It should be possible to vary the increments ℎ and 𝑘 (e.g., ℎ → 0 and 𝑘 → 0 must be allowed when the convergent rates are to be studied). It will be assumed that the ratio ℎ / 𝑘 remains constant when the increments are varied. This implies a requirement that if, for example, ℎ is halved then 𝑘 is also halved. More precisely, it will be assumed that if an arbitrary pair of increments ( ℎ 1 , 𝑘 1 ) is replaced with another pair ( ℎ 2 , 𝑘 2 ) , then the following relationship holds: 𝑘 1 ℎ 1 = 𝑘 2 ℎ 2 = 𝛾 , ( 1 . 6 ) where 𝛾 is a constant which does not depend on increments. The requirement to keep ℎ / 𝑘 constant is very reasonable.

Assume that 𝑐 ( 𝑥 𝑖 , 𝑡 𝑛 ) is the exact solution of (1.1) at an arbitrary grid-point belonging to the set defined by the two grids (1.4) and (1.5). Then the values 𝑐 𝑖 , 𝑛 ( 𝑖 = 0 , 1 , … , 𝑁 𝑥 and 𝑛 = 1 , 2 , … , 𝑁 𝑡 ) calculated by (1.2) are approximations of the exact solution at the grid-points ( 𝑥 𝑖 , 𝑡 𝑛 ) . Our major task in the following part of this paper will be to show how the accuracy of the approximations 𝑐 𝑖 , 𝑛 can be improved by using additionally the Richardson Extrapolation.

The application of the Richardson Extrapolation, when an arbitrary advection equation (not only the particular equation which was introduced in the beginning of this section) is treated by any numerical method, will be described in Section 2. The error constants in the leading terms of the numerical error for the Crank-Nicolson scheme will be calculated in Section 3. The order of accuracy of the combination of the Crank-Nicolson scheme and the Richardson Extrapolation will be established in Section 4. Numerical results will be presented in Section 5 to demonstrate the applicability of the theorems proved in Sections 3 and 4. Several concluding remarks will be given and discussed in Section 6.

2. Application of the Richardson Extrapolation

Assume that a one-dimensional hyperbolic equation similar to (1.1) is treated by an arbitrary numerical method, which is of order 𝑝 ≥ 1 with regard to the two independent variables 𝑥 and 𝑡 . Let { 𝑧 𝑖 , 𝑛 + 1 } 𝑁 𝑥 𝑖 = 0 be the set of approximations of the solution of (1.1) calculated for 𝑡 = 𝑡 𝑛 + 1 ∈ 𝐺 𝑡 at all grid-points 𝑥 𝑖 , 𝑖 = 0 , 1 , … , 𝑁 𝑥 , of (1.4) by using the numerical method chosen and consider the corresponding approximations { 𝑧 𝑖 , 𝑛 } 𝑁 𝑥 𝑖 = 0 calculated at the previous time-step, that is, for 𝑡 = 𝑡 𝑛 ∈ 𝐺 𝑡 . Introduce vectors 𝑐 ( 𝑡 𝑛 + 1 ) , 𝑧 𝑛 , and 𝑧 𝑛 + 1 , the components of which are { 𝑐 ( 𝑥 𝑖 , 𝑡 𝑛 + 1 ) } 𝑁 𝑥 𝑖 = 0 , { 𝑧 𝑖 , 𝑛 } 𝑁 𝑥 𝑖 = 0 , and { 𝑧 𝑖 , 𝑛 + 1 } 𝑁 𝑥 𝑖 = 0 , respectively. Since the order of the numerical method is assumed to be 𝑝 with regard both to 𝑥 and to 𝑡 , we can write 𝑐 𝑡 𝑛 + 1 = 𝑧 𝑛 + 1 + ℎ 𝑝 𝐾 1 + 𝑘 𝑝 𝐾 2 𝑘 + 𝑂 𝑝 + 1 , ( 2 . 1 ) where 𝐾 1 and 𝐾 2 are some quantities, which do not depend on ℎ and on 𝑘 . It is convenient to rewrite the last equality in the following form: 𝑐 𝑡 𝑛 + 1 = 𝑧 𝑛 + 1 + 𝑘 𝑝 𝑘 𝐾 + 𝑂 𝑝 + 1 , ( 2 . 2 ) where ℎ 𝐾 = 𝑘 𝑝 𝐾 1 + 𝐾 2 . ( 2 . 3 ) If ℎ and 𝑘 are sufficiently small, then the sum ℎ 𝑝 𝐾 1 + 𝑘 𝑝 𝐾 2 will be a good approximation of the errors in the calculated values of the numerical solution { 𝑧 𝑖 , 𝑛 + 1 } 𝑁 𝑥 𝑖 = 0 . If the relationship (1.6) is satisfied and if ℎ and 𝑘 are sufficiently small, then 𝑘 𝑝 𝐾 will also be a good approximation of the error of the numerical solution 𝑧 𝑛 + 1 . This means that if we succeed to eliminate the term 𝑘 𝑝 𝐾 in (2.2), then we shall obtain approximations of order 𝑝 + 1 . The Richardson Extrapolation can be applied in the attempt to achieve such an improvement of the accuracy. In order to apply the Richardson Extrapolation, it is necessary to introduce an additional grid: 𝐺 2 𝑥 = 𝑥 𝑖 , 𝑖 = 0 , 1 , … , 2 𝑁 𝑥 ∣ 𝑥 0 = 𝑎 1 , 𝑥 𝑖 = 𝑥 𝑖 − 1 + ℎ 2 , 𝑖 = 1 , 2 , … , 2 𝑁 𝑥 𝑏 , ℎ = 1 − 𝑎 1 𝑁 𝑥 , 𝑥 2 𝑁 𝑥 = 𝑏 1 . ( 2 . 4 ) Assume that approximations { 𝑤 𝑖 , 𝑛 } 2 𝑁 𝑥 𝑖 = 0 (calculated at the grid-points of 𝐺 2 𝑥 for 𝑡 = 𝑡 𝑛 ∈ 𝐺 𝑡 ) are available and perform two small steps with a stepsize 𝑘 / 2 to compute { 𝑤 𝑖 , 𝑛 + 1 } 2 𝑁 𝑥 𝑖 = 0 . Use only the components with even indices 𝑖 , 𝑖 = 0 , 2 , 4 , … , 2 𝑁 𝑥 , to form vector 𝑤 𝑛 + 1 . The following equality holds for this vector: 𝑐 𝑡 𝑛 + 1 = 𝑤 𝑛 + 1 + 𝑘 2 𝑝 𝑘 𝐾 + 𝑂 𝑝 + 1 , ( 2 . 5 ) where the quantity 𝐾 is defined as in (2.3).

Now, it is possible to eliminate the quantity 𝐾 from (2.2) and (2.5). This can successfully be done in the following way: ( a ) remove the last two terms in (2.2) and (2.5), ( b ) multiply (2.5) by 2 𝑝 , ( c ) subtract (2.2) from the modified equality (2.5).

The result is 𝑐 𝑡 𝑛 + 1 = 2 𝑝 𝑤 𝑛 + 1 − 𝑧 𝑛 + 1 2 𝑝 𝑘 − 1 + 𝑂 𝑝 + 1 . ( 2 . 6 ) Denote 𝑐 𝑛 + 1 = 2 𝑝 𝑤 𝑛 + 1 − 𝑧 𝑛 + 1 2 𝑝 − 1 . ( 2 . 7 ) It is clear that the approximation 𝑐 𝑛 + 1 , being of order 𝑝 + 1 , will in general be more accurate than both 𝑧 𝑛 + 1 and 𝑤 𝑚 + 1 when ℎ and 𝑘 are sufficiently small. The device used to construct the approximation 𝑐 𝑛 + 1 is often called Richardson Extrapolation (it was introduced in [2]; see also [3–7]). Assume that the partial derivatives of the unknown function 𝑐 ( 𝑥 , 𝑡 ) up to order 𝑝 + 1 exist and are continuous. Then one should expect (2.7) to produce more accurate results than those obtained by the underlying numerical method.

Remark 2.1. The rest terms in the formulae given in this section will in general depend on both ℎ and 𝑘 . However, the application of the relationship (1.6) gives ℎ = 𝑘 / 𝛾 , and this justifies the use only of 𝑘 in all rest terms in this section. This approach will also be used in the remaining part of this paper.

Remark 2.2. No specific assumptions about the particular partial differential equation or about the numerical method used to solve it were made in this section. This was done in order to demonstrate how general the idea, on which the Richardson Extrapolation is based, is. It must be emphasized, however, that in the following part of this paper it will be assumed that (a) equation (1.1) is solved under the assumptions made in the previous section and (b) the underlying numerical algorithm applied to handle it numerically will always be the second-order Crank-Nicolson scheme.

3. Error Constants of the Leading Terms of the Numerical Error for the Crank-Nicolson Scheme

Consider formula (1.2). Following [8], we shall replace the approximations 𝑐 𝑖 , 𝑛 with the corresponding exact values 𝑐 ( 𝑥 𝑖 , 𝑡 𝑛 ) of the solution of (1.1). The result is 𝐿 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 ; ℎ , 𝑘 = 𝜎 𝑖 , 𝑛 + 0 . 5 𝑐 𝑥 𝑖 + 1 , 𝑡 𝑛 + 1 𝑥 + 𝑐 𝑖 , 𝑡 𝑛 + 1 − 𝜎 𝑖 , 𝑛 + 0 . 5 𝑐 𝑥 𝑖 − 1 , 𝑡 𝑛 + 1 + 𝜎 𝑖 , 𝑛 + 0 . 5 𝑐 𝑥 𝑖 + 1 , 𝑡 𝑛 𝑥 − 𝑐 𝑖 , 𝑡 𝑛 − 𝜎 𝑖 , 𝑛 + 0 . 5 𝑐 𝑥 𝑖 − 1 , 𝑡 𝑛 = 𝜎 𝑖 , 𝑛 + 0 . 5 𝑐 𝑥 𝑖 + 1 , 𝑡 𝑛 + 1 𝑥 − 𝑐 𝑖 − 1 , 𝑡 𝑛 + 1 + 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 1 𝑥 − 𝑐 𝑖 , 𝑡 𝑛 + 𝜎 𝑖 , 𝑛 + 0 . 5 𝑐 𝑥 𝑖 + 1 , 𝑡 𝑛 𝑥 − 𝑐 𝑖 − 1 , 𝑡 𝑛 . ( 3 . 1 ) The following theorem can be proved by using this notation in the following theorem.

Theorem 3.1. The quantity 𝐿 [ 𝑐 ( 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 ) ; ℎ , 𝑘 ] can be written (assuming that all involved derivatives exist and are continuous) as 𝐿 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 = ; ℎ , 𝑘 𝛾 ℎ 3 6 𝑢 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 𝑥 3 + 𝑘 3 𝜕 2 4 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 𝑡 3 + 𝑘 3 8 𝑢 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 𝑥 𝜕 𝑡 2 𝑘 + 𝑂 5 . ( 3 . 2 )

Proof. Use Taylor expansions around the point ( 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 ) , where 𝑡 𝑛 + 0 . 5 = 𝑡 𝑛 + 0 . 5 𝑘 as in Section 1, of the functions in two variables involved in (3.1) and keep the terms containing 𝑘 𝑠 , 𝑠 = 0 , 1 , 2 , 3 , 4 . After some rather long but quite straight-forward transformations, (3.2) will be obtained.

Theorem 3.1 ensures that the Crank-Nicolson scheme is a second-order numerical method, which is, of course, well known. It is much more important that (a) it provides the leading terms of the error of this method (which are needed in the proof of Theorem 4.1) and (b) it shows that there are no fourth-order terms in the expression for the numerical error.

4. Order of Accuracy of the Combination Consisting of the Crank-Nicolson Scheme and the Richardson Extrapolation

After presenting some preliminary results connected to (a) the problem solved, (b) the Crank-Nicolson scheme, and (c) the Richardson Extrapolation, everything is now prepared for the proof that the use of the combination of the Crank-Nicolson scheme and the Richardson Extrapolation leads to a fourth-order numerical method when the problem (1.1) is solved. More precisely, the following theorem holds.

Theorem 4.1. The combination of the Crank-Nicolson scheme and the Richardson Extrapolation is a fourth-order numerical method when all derivatives of the solution up to order four exist and are continuous.

Proof. Assume that all approximations { 𝑐 𝑖 , 𝑛 } 𝑁 𝑥 𝑖 = 0 at time-step 𝑛 have already been found. Then the Richardson Extrapolation is carried out by using the Crank-Nicolson scheme to perform one large time-step with stepsize 𝑘 and two small time-steps with stepsize 0 . 5 𝑘 . The major part of the computations during the large time-step is carried out by using the formula: 𝜎 𝑖 , 𝑛 + 0 . 5 𝑧 𝑖 + 1 , 𝑛 + 1 + 𝑧 𝑖 , 𝑛 + 1 − 𝜎 𝑖 , 𝑛 + 0 . 5 𝑧 𝑖 − 1 , 𝑛 + 1 + 𝜎 𝑖 , 𝑛 + 0 . 5 𝑐 𝑖 + 1 , 𝑛 − 𝑐 𝑖 , 𝑛 − 𝜎 𝑖 , 𝑛 + 0 . 5 𝑐 𝑖 − 1 , 𝑛 = 0 . ( 4 . 1 ) The major part of the computations during the two small time-steps is based on the use of the following two formulas: 𝜎 𝑖 , 𝑛 + 0 . 2 5 𝑤 𝑖 + 0 . 5 , 𝑛 + 0 . 5 + 𝑤 𝑖 , 𝑛 + 0 . 5 − 𝜎 𝑖 , 𝑛 + 0 . 2 5 𝑤 𝑖 − 0 . 5 , 𝑛 + 0 . 5 + 𝜎 𝑖 , 𝑛 + 0 . 2 5 𝑐 𝑖 + 0 . 5 , 𝑛 − 𝑐 𝑖 , 𝑛 − 𝜎 𝑖 , 𝑛 + 0 . 2 5 𝑐 𝑖 − 0 . 5 , 𝑛 𝜎 = 0 , 𝑖 , 𝑛 + 0 . 7 5 𝑤 𝑖 + 0 . 5 , 𝑛 + 1 + 𝑤 𝑖 , 𝑛 + 1 − 𝜎 𝑖 , 𝑛 + 0 . 7 5 𝑤 𝑖 − 0 . 5 , 𝑛 + 1 + 𝜎 𝑖 , 𝑛 + 0 . 7 5 𝑤 𝑖 + 0 . 5 , 𝑛 + 0 . 5 − 𝑤 𝑖 , 𝑛 + 0 . 5 − 𝜎 𝑖 , 𝑛 + 0 . 7 5 𝑤 𝑖 − 0 . 5 , 𝑛 + 0 . 5 = 0 . ( 4 . 2 ) Let us start with (4.1). Equality (3.1) will be obtained when all approximate values in (4.1) are replaced with the corresponding values of the exact solution. This means that the assertion of Theorem 3.1, equality (3.2), holds for the large time-step.
The treatment of the two small time-steps is more complicated. Combining (4.2) leads to the formula: 𝜎 𝑖 , 𝑛 + 0 . 7 5 𝑤 𝑖 + 0 . 5 , 𝑛 + 1 + 𝑤 𝑖 , 𝑛 + 1 − 𝜎 𝑖 , 𝑛 + 0 . 7 5 𝑤 𝑖 − 0 . 5 , 𝑛 + 1 + 𝜎 𝑖 , 𝑛 + 0 . 7 5 𝑤 𝑖 + 0 . 5 , 𝑛 + 0 . 5 − 𝑤 𝑖 , 𝑛 + 0 . 5 − 𝜎 𝑖 , 𝑛 + 0 . 7 5 𝑤 𝑖 − 0 . 5 , 𝑛 + 0 . 5 + 𝜎 𝑖 , 𝑛 + 0 . 2 5 𝑤 𝑖 + 0 . 5 , 𝑛 + 0 . 5 + 𝑤 𝑖 , 𝑛 + 0 . 5 − 𝜎 𝑖 , 𝑛 + 0 . 2 5 𝑤 𝑖 − 0 . 5 , 𝑛 + 0 . 5 + 𝜎 𝑖 , 𝑛 + 0 . 2 5 𝑐 𝑖 + 0 . 5 , 𝑛 − 𝑐 𝑖 , 𝑛 − 𝜎 𝑖 , 𝑛 + 0 . 2 5 𝑐 𝑖 − 0 . 5 , 𝑛 = 0 . ( 4 . 3 ) Replace all approximate values participating in (4.3) with the corresponding exact values of the solution to obtain an expression for the local approximation error 𝐿 in the form: 𝐿 𝐿 = 1 + 𝐿 2 , ( 4 . 4 ) where 𝐿 1 = 𝜎 𝑖 , 𝑛 + 0 . 7 5 𝑐 𝑥 𝑖 + 0 . 5 , 𝑡 𝑛 + 1 𝑥 − 𝑐 𝑖 − 0 . 5 , 𝑡 𝑛 + 1 𝑥 + 𝑐 𝑖 + 0 . 5 , 𝑡 𝑛 + 0 . 5 𝑥 − 𝑐 𝑖 − 0 . 5 , 𝑡 𝑛 + 0 . 5 𝑥 + 𝑐 𝑖 , 𝑡 𝑛 + 1 𝑥 − 𝑐 𝑖 , 𝑡 𝑛 + 0 . 5 , 𝐿 2 = 𝜎 𝑖 , 𝑛 + 0 . 2 5 𝑐 𝑥 𝑖 + 0 . 5 , 𝑡 𝑛 + 0 . 5 𝑥 − 𝑐 𝑖 − 0 . 5 , 𝑡 𝑛 + 0 . 5 𝑥 + 𝑐 𝑖 + 0 . 5 , 𝑡 𝑛 𝑥 − 𝑐 𝑖 − 0 . 5 , 𝑡 𝑛 𝑥 + 𝑐 𝑖 , 𝑡 𝑛 + 0 . 5 𝑥 − 𝑐 𝑖 , 𝑡 𝑛 . ( 4 . 5 ) Our aim is to derive an expression for 𝐿 .
First, we consider the terms participating in 𝐿 1 . We use Taylor expansions of the involved functions around the point ( 𝑥 𝑖 , 𝑡 𝑛 + 0 . 7 5 ) and apply a similar transformation as in the proof of Theorem 3.1. In this way, we can obtain the following relation: 𝐿 1 = 𝛾 ℎ 3 𝑢 𝑥 4 8 𝑖 , 𝑡 𝑛 + 0 . 7 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 7 5 𝜕 𝑥 3 + 𝑘 3 𝜕 1 9 2 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 7 5 𝜕 𝑡 3 + 𝑘 3 𝑢 𝑥 6 4 𝑖 , 𝑡 𝑛 + 0 . 7 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 7 5 𝜕 𝑥 𝜕 𝑡 2 𝑘 + 𝑂 5 . ( 4 . 6 ) We repeat the same kind of transformations also when 𝐿 2 is considered. Now we apply Taylor expansions around the point ( 𝑥 𝑖 , 𝑡 𝑛 + 0 . 2 5 ) of the involved functions. Then we obtain 𝐿 2 = 𝛾 ℎ 3 𝑢 𝑥 4 8 𝑖 , 𝑡 𝑛 + 0 . 2 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 2 5 𝜕 𝑥 3 + 𝑘 3 𝜕 1 9 2 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 2 5 𝜕 𝑡 3 + 𝑘 3 𝑢 𝑥 6 4 𝑖 , 𝑡 𝑛 + 0 . 2 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 2 5 𝜕 𝑥 𝜕 𝑡 2 𝑘 + 𝑂 5 . ( 4 . 7 ) The following result can be found by combining (4.6) and (4.7): 𝐿 = 𝛾 ℎ 3 𝑢 𝑥 4 8 𝑖 , 𝑡 𝑛 + 0 . 7 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 7 5 𝜕 𝑥 3 𝑥 + 𝑢 𝑖 , 𝑡 𝑛 + 0 . 2 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 2 5 𝜕 𝑥 3 + 𝑘 3 𝜕 1 9 2 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 7 5 𝜕 𝑡 3 + 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 2 5 𝜕 𝑡 3 + 𝑘 3 𝑢 𝑥 6 4 𝑖 , 𝑡 𝑛 + 0 . 7 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 7 5 𝜕 𝑥 𝜕 𝑡 2 𝑥 + 𝑢 𝑖 , 𝑡 𝑛 + 0 . 2 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 2 5 𝜕 𝑥 𝜕 𝑡 2 𝑘 + 𝑂 5 . ( 4 . 8 ) Now, by expanding all terms in the right-hand side of (4.8) around the point ( 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝑘 ) and after some very long but straight-forward transformations, we obtain 𝐿 = 𝛾 ℎ 3 𝑢 𝑥 2 4 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 𝑥 3 + 𝑘 3 𝜕 9 6 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 𝑡 3 + 𝑘 3 𝑢 𝑥 3 2 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 𝑥 𝜕 𝑡 3 𝑘 + 𝑂 5 . ( 4 . 9 ) It should be noted here that a detailed derivation of the important relationship (4.9) can be found in [9].Since the order of the Crank-Nicolson scheme is 𝑝 = 2 , it is clear that the improved approximate solution by the Richardson Extrapolation at time-step 𝑛 + 1 is obtained by ( a ) multiplying the result obtained at the end of the second small time-step by 4 / 3 , ( b ) multiplying the result obtained at the end of the large time-step by 1 / 3 , ( c ) subtracting the two results obtained in (a) and (b).Performing operations (a)–(c) will give 4 3 1 𝐿 − 3 𝐿 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 = ; ℎ , 𝑘 𝛾 ℎ 3 𝑢 𝑥 1 8 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 𝑥 3 + 𝑘 3 𝜕 7 2 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 𝑡 3 + 𝑘 3 𝑢 𝑥 2 4 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 𝑥 𝜕 𝑡 3 − 𝛾 ℎ 3 𝑢 𝑥 1 8 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 𝑥 3 − 𝑘 3 𝜕 7 2 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 𝑡 3 − 𝑘 3 𝑢 𝑥 2 4 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 3 𝑐 𝑥 𝑖 , 𝑡 𝑛 + 0 . 5 𝜕 𝑥 𝜕 𝑡 2 𝑘 + 𝑂 5 . ( 4 . 1 0 ) It is immediately seen that the first six terms in the right-hand side of (4.10) vanish. Therefore, the order of accuracy of the combined numerical method (the Crank-Nicolson scheme + the Richardson Extrapolation) is four, which completes the proof of the theorem.

It should once again be emphasized that a full proof of Theorem 4.1, containing all needed details, can be found in [9].

5. Numerical Experiments

In this section, it will be shown that the following two statements are true:(a)if the solution is continuously differentiable up to order two, then the direct application of the Crank-Nicolson scheme gives second-order accuracy;(b)if the solution is continuously differentiable up to order four, then the combined method consisting of the Crank-Nicolson scheme and the Richardson Extrapolation behaves as a fourth-order numerical algorithm.

Furthermore, we shall also demonstrate the fact that if the above requirements are not satisfied, then neither the direct use of the Crank-Nicolson scheme leads to second-order accuracy nor the combined method based on the combination of the Crank-Nicolson scheme with the Richardson Extrapolation behaves as a fourth-order numerical algorithm.

5.1. Organization of the Computations

It is convenient for the purposes in this paper, but not necessary, to divide the time-interval [a,b] into 24 equal subintervals and to call each of these subintervals “hour”. By this convention, the length of the time-interval becomes 24 hours in all three examples given in this section; and we shall study the size of the numerical errors at the end of every hour.

In each experiment, the first run is performed by using 𝑁 𝑡 = 1 6 8 and 𝑁 𝑥 = 1 6 0 . Ten additional runs are performed after the first one. When a run is finished, both 𝑘 and ℎ are halved (this means that 𝑁 𝑡 and 𝑁 𝑥 are doubled) and a new run is started. Thus, in the eleventh run we have 𝑁 𝑡 = 1 7 2 0 3 2 and 𝑁 𝑥 = 1 6 3 8 4 0 which means that 1 7 2 0 3 2 systems of linear algebraic equations, each of them containing 1 6 3 8 4 0 equations, are to be solved.

Note too, that the ratio ℎ / 𝑘 is kept constant, that is, the requirement (1.6) is satisfied, when the two increments 𝑘 and ℎ are varied in this way.

We are mainly interested in the behavior of the numerical error. As mentioned above, this error is evaluated at the end of every hour (i.e., 24 times in each run) and at the grid-points of the coarsest spatial grid, in the following way. Assume that run number 𝑟 , 𝑟 = 1 , 2 , … , 1 1 , is to be carried out and let 𝑅 = 2 𝑟 − 1 . Then the error made at the end of hour 𝑚 is calculated by using the following formula: E R R 𝑚 = m a x 𝑗 = 0 , 1 , … , 1 6 0 ⎛ ⎜ ⎜ ⎜ ⎝ | | | 𝑐 ̃ 𝑖 , ̃ 𝑛 − 𝑐 e x a c t ̃ 𝑖 , | | | ̃ 𝑛 | | | 𝑐 m a x e x a c t ̃ 𝑖 , | | | ⎞ ⎟ ⎟ ⎟ ⎠ ̃ ̃ 𝑛 , 1 . 0 , 𝑚 = 1 , 2 , … , 2 4 , 𝑖 = 𝑗 𝑅 , ̃ 𝑛 = 7 𝑚 𝑅 , ( 5 . 1 ) where 𝑐 ̃ 𝑖 , ̃ 𝑛 and 𝑐 e x a c t ̃ 𝑖 , ̃ 𝑛 are the calculated value and the exact solution at the end of hour 𝑚 and at the grid-points of the coarsest grid (i.e., for 𝑁 𝑥 = 1 6 0 ). In the three experiments presented in this section, the exact solution is known.

The global error made during the computations is estimated by using the following formula: E R R = m a x 𝑚 = 1 , 2 , … , 2 4 E R R 𝑚 . ( 5 . 2 ) It is necessary to point out here that the numerical values of the unknown function, which are improved by the Richardson Extrapolation, that is, by applying (2.7) with 𝑝 = 2 , are available only on the coarser spatial grid (1.4). It is necessary to get appropriate approximations for all values on the finer spatial grid (2.4). Several devices for obtaining such approximations have been suggested and tested in [10]. It was shown there that the application of third-order interpolation polynomials gives best results. This device has been used also in the present work.

5.2. Construction of a Test Problem with Steep Gradients of the Unknown Function

Assume that the spatial and the time intervals are given by [ ] [ ] 𝑥 ∈ 0 , 5 0 0 0 0 0 0 0 , 𝑡 ∈ 4 3 2 0 0 , 1 2 9 6 0 0 , ( 5 . 3 ) and consider a function 𝑢 ( 𝑥 , 𝑡 ) defined by 𝑢 ( 𝑥 , 𝑡 ) = 3 2 0 . ( 5 . 4 ) Let the initial condition be given by 𝑓 ( 𝑥 ) = 𝜉 1 + 9 9 . 0 𝑒 − 𝜔 ( 𝑥 − 1 0 0 0 0 0 0 0 ) 2 , 𝜉 = 1 . 4 6 7 9 ⋅ 1 0 1 2 , 𝜔 = 1 0 − 1 2 . ( 5 . 5 ) The exact solution of the Test Problem, which is defined as above, is 𝑐 ( 𝑥 , 𝑡 ) = 𝑓 ( 𝑥 − 3 2 0 ( 𝑡 − 4 3 2 0 0 ) ) . ( 5 . 6 ) The Test Problem introduced in this subsection was run both by using the Crank-Nicolson scheme directly and by applying the combination of this scheme and the Richardson Extrapolation. Numerical results are presented in Table 1.

tab1
Table 1: Results obtained when the test problem defined by (5.3)–(5.5) is handled directly by the Crank-Nicolson scheme and by using the combination of the Crank-Nicolson scheme and the Richardson Extrapolation. The numerical errors calculated by (5.1) and (5.2) are given in the columns under “Error”. In row 𝑠 , 𝑠 = 2 , 3 , … , 1 1 , the ratios of the errors in this row and in the previous row are given in the columns under “Ratio”.

The following conclusions can be drawn by studying the results presented in Table 1.(i)The direct application of the Crank-Nicolson scheme leads to quadratic convergence of the accuracy of the numerical results (i.e., halving the increments 𝑘 and ℎ leads to a decrease of the error by a factor of four). This behaviour should be expected according to Theorem 3.1.(ii)The combination of the Crank-Nicolson scheme and the Richardson Extrapolation behaves in general as a numerical method of order four (or, in other words, halving the increments ℎ and 𝑘 leads to a decrease of the error by a factor of sixteen). This behaviour should also be expected (according to Theorem 4.1).(iii)At the end of the computations with the combined numerical method (the Crank-Nicolson scheme + the Richardson Extrapolation), the convergence rate deteriorates. Two facts are very important when this happens: (a) the computed solution is already very accurate and (b) the rounding errors start to affect the calculated results.

In Figure 1, plots are given showing (a) the initial values, (b) the solution in the middle of the time interval (i.e., after 12 hours), and (c) the solution at the end of the time interval for the test problem defined by (5.3)–(5.5).

fig1
Figure 1: The solution of the one-dimensional advection equation defined in Section 5.2: (a) at the beginning of the interval (the upper plot), (b) at the end of the twelfth hour (the plot in the middle), and (c) at the end of the time-interval (the lower plot).

Remark 5.1. A similar test-example was used in [11]. It should also be noted that a very similar advection module is a part of the large-scale air pollution model UNI-DEM ([12–14]) and the quantities used in (5.3)–(5.5) are either the same or very similar to the corresponding quantities in this model. Note too that the values of the unknown function are of the same order of magnitude as the ozone concentrations in the atmosphere when these are measured in (number of molecules)/(cubic centimetre).

5.3. Construction of an Oscillatory Test Problem

Define the spatial and time intervals by 𝑎 = 𝑎 1 = 0 , 𝑏 = 𝑏 1 = 2 𝜋 . ( 5 . 7 ) Consider a function 𝑢 ( 𝑥 , 𝑡 ) defined by 𝑢 ( 𝑥 , 𝑡 ) = 0 . 5 . ( 5 . 8 ) Let the initial values be defined by 𝑓 [ ] ( 𝑥 ) = 1 0 0 + 9 9 s i n ( 1 0 𝑥 ) ∗ 1 . 4 6 7 9 ∗ 1 0 1 2 . ( 5 . 9 ) The exact solution of the test problem defined by (5.7)–(5.9) is: 𝑐 ( 𝑥 , 𝑡 ) = 𝑓 ( 𝑥 − 0 . 5 𝑡 ) . ( 5 . 1 0 ) As in Section 5.2, the test problem introduced above was run both by using the Crank-Nicolson scheme directly and by applying the combination of this scheme and the Richardson Extrapolation. Numerical results are presented in Table 2.

tab2
Table 2: Results obtained when the test problem defined by (5.7)–(5.9) is handled directly by the Crank-Nicolson scheme and by using the combination of the Crank-Nicolson scheme and the Richardson Extrapolation. The numerical errors calculated by (5.1) and (5.2) are given in the columns under “Error”. In row 𝑠 , 𝑠 = 2 , 3 , … , 1 1 , the ratios of the errors in this row and in the previous row are given in the columns under “Ratio”.

The conclusions, which can be drawn from the results presented in Table 2, are quite similar to those given in Section 5.2. However, for the oscillatory test problem, the actual convergence rate achieved in the runs is less than four (greater than three in the beginning and after that equal to or less than three). It is not very clear what the reason for this behaviour is. Perhaps, the interpolation rule used to improve the precision of the values on the finer spatial grid (see Section 5.1 and [10]) is not sufficiently accurate for this example when grid-points near the boundary are treated. Nevertheless, it is clearly seen that the achieved accuracy is nearly the same as the accuracy achieved in the solution of the previous test problem (compare Table 1 with Table 2).

In Figure 2, plots are given showing (a) the initial values, (b) the solution in the middle of the time interval (i.e., after 12 hours), and (c) the solution at the end of the time interval for the test problem studied in this subsection.

fig2
Figure 2: The solution of the one-dimensional advection equation defined in Section 5.3 is drawn (a) at the beginning of the interval (the upper plot), (b) at the end of the twelfth hour (the plot in the middle), and (c) at the end of the time interval (the lower plot).
5.4. Construction of a Test Problem with Discontinuities

Assume that the spatial interval, the time-interval, and function 𝑢 ( 𝑥 , 𝑡 ) are defined as in Section 5.2, that is, by (5.3) and (5.4), and introduce initial values by 𝑓 ( 𝑥 ) = 1 . 4 6 7 9 ∗ 1 0 1 2 , 𝑥 ≤ 5 0 0 0 0 0 0 o r 𝑥 ≥ 1 5 0 0 0 0 0 0 , 𝑓 ( 𝑥 ) = 1 + 9 9 ∗ 𝑥 − 5 0 0 0 0 0 0 5 0 0 0 0 0 0 ∗ 1 . 4 6 7 9 ∗ 1 0 1 2 𝑓 , 5 0 0 0 0 0 0 ≤ 𝑥 ≤ 1 0 0 0 0 0 0 0 , ( 𝑥 ) = 1 + 9 9 ∗ 1 5 0 0 0 0 0 0 − 𝑥 5 0 0 0 0 0 0 ∗ 1 . 4 6 7 9 ∗ 1 0 1 2 , 1 0 0 0 0 0 0 0 ≤ 𝑥 ≤ 1 5 0 0 0 0 0 0 . ( 5 . 1 1 ) The exact solution of the test problem defined as above is given by (5.6).

As in the previous two sub-sections, the test problem introduced above was run both by using the Crank-Nicolson scheme directly and by applying the combination of this scheme and the Richardson Extrapolation. Numerical results are presented in Table 3.

tab3
Table 3: Results obtained when the test problem defined in Section 5.4 is handled directly by the Crank-Nicolson scheme and by using the combination of the Crank-Nicolson scheme and the Richardson Extrapolation. The numerical errors calculated by (5.1) and (5.2) are given in the columns under “Error”. In row 𝑠 , 𝑠 = 2 , 3 , … , 1 1 , the ratios of the errors in this row and in the previous row are given in the columns under “Ratio”.

Two major conclusions can be drawn from the results presented in Table 3: (a) neither the direct Crank-Nicolson scheme nor the combination of the Crank-Nicolson scheme with the Richardson Extrapolation gives the prescribed by the theory accuracy (orders two and four, resp.) and (b) also in this case, that is, in the presence of discontinuities, the combination of the Crank-Nicolson scheme and the Richardson Extrapolation is considerably more accurate than the direct application of the Crank-Nicolson scheme.

In Figure 3, plots are given showing (a) the initial values, (b) the solution in the middle of the time interval (i.e., after 12 hours), and (c) the solution at the end of the time interval for the test problem studied in this subsection.

fig3
Figure 3: The solution of the one-dimensional advection equation defined in Section 5.4 is drawn (a) at the beginning of the interval (the upper plot), (b) at the end of the twelfth hour (the plot in the middle), and (c) at the end of the time interval (the lower plot).

6. Concluding Remarks

Several conclusions can be drawn by using the theorems proved in Sections 3 and 4 as well as the numerical results presented in Section 6.

The most important conclusion is related to the accuracy of the computed results. The accuracy can be improved considerably if the Crank-Nicolson scheme is combined with the Richardson Extrapolation. However, this effect will be achieved only when the solution is continuously differentiable up to order four and when the boundary conditions needed in the actual computations are sufficiently accurate.

It is highly desirable to investigate carefully the performance of the combined method in the cases where ( a ) some derivatives of the solution are discontinuous, ( b ) the solution is highly oscillatory, ( c ) the problem is stiff (which will normally cause stability problems).

These important topics will be studied in the near future.

Some assumptions were made in order to prove the two theorems. It is worthwhile to investigate carefully the possibilities for removing these assumptions or for relaxing them.

Acknowledgments

The Centre for Supercomputing at the Technical University of Denmark gave us access to several powerful computers for running the experiments related to this study. The European Union and the European Social Fund have provided financial support to the project under the Grant agreement no. TÁMOP 4.2.1./B-09/1/KMR-2010-0003. The National Science Fund of Bulgaria supported partly the research of I. Dimov and T. Ostromsky under Grant DTK 02/44 /2009/ and the research of K. Georgiev under Grant DO 02-147/08.

References

  1. J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, SIAM, Philadelphia, Pa, USA, 2nd edition, 2004. View at Zentralblatt MATH
  2. L. F. Richardson, “The deferred approach to the limit, I—single lattice,” Philosophical Transactions of the Royal Society A, vol. 226, pp. 299–349, 1927. View at Publisher · View at Google Scholar
  3. Z. Zlatev, I. Faragó, and Á. Havasi, “On some stability properties of the Richardson extrapolation applied together with the θ-method,” in Large Scale Scientific Computing, I. Lirkov, S. Margenov, and J. Wasniewski, Eds., vol. 5910 of Lecture Notes in Computer Science, pp. 54–65, Spinger, Berlin, Germany, 2010. View at Publisher · View at Google Scholar
  4. C. Burg and T. Erwin, “Application of Richardson extrapolation to the numerical solution of partial differential equations,” Numerical Methods for Partial Differential Equations, vol. 25, no. 4, pp. 810–832, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  5. I. Faragó, Á. Havasi, and Z. Zlatev, “Efficient implementation of stable Richardson extrapolation algorithms,” Computers and Mathematics with Applications, vol. 60, no. 8, pp. 2309–2325, 2010. View at Publisher · View at Google Scholar · View at Scopus
  6. Z. Zlatev, I. Faragó, and Á. Havasi, “Stability of the Richardson extrapolation applied together with the θ-method,” Journal of Computational and Applied Mathematics, vol. 235, no. 2, pp. 507–522, 2010. View at Publisher · View at Google Scholar · View at Scopus
  7. S. A. Chin and J. Geiser, “Multi-product operator splitting as a general method of solving autonomous and nonautonomous equations,” IMA Journal of Numerical Analysis. In press. View at Publisher · View at Google Scholar
  8. J. D. Lambert, Numerical Methods for Ordinary Differential Equations: The Initial Value Problem, John Wiley & Sons, New York, NY, USA, 1991.
  9. Z. Zlatev, I. Dimov, I. Faragó, K. Georgiev, Á Havasi, and T. Ostromsky, “Solving advection equations by applying the Crank-Nicolson scheme combined with the Richardson extrapolation (extended version),” 2011, http://nimbus.elte.hu/~hagi/IJDE/.
  10. Z. Zlatev, I. Dimov, I. Faragó, K. Georgiev, Á Havasi, and T. Ostromsky, “Implementation of Richardson extrapolation in the treatment of one-dimensional advection equations,” in Numerical Methods and Applications, I. Dimov, S. Dimova, and N. Kolkovska, Eds., vol. 6046 of Lecture Notes in Computer Science, pp. 198–206, Springer, Berlin, Germany, 2011. View at Publisher · View at Google Scholar
  11. Z. Zlatev, R. Berkowicz, and L. P. Prahm, “Testing subroutines solving advection-diffusion equations in atmospheric environments,” Computers and Fluids, vol. 11, no. 1, pp. 13–38, 1983. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  12. Z. Zlatev, Computer Treatment of Large Air Pollution Models, Kluwer Academic Publishers, London, UK, 1995.
  13. Z. Zlatev and I. Dimov, Computational and Numerical Challenges in Environmental Modelling, Elsevier, Amsterdam, The Netherlands, 2006.
  14. V. N. Alexandrov, W. Owczarz, P. G. Thomson, and Z. Zlatev, “Parallel runs of a large air pollution model on a grid of SUN computers,” Mathematics and Computers in Simulation, vol. 65, no. 6, pp. 557–577, 2004. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus