Fourth Order Partitioned Methods Designed for the Time Integration of Atmospheric Models

Göran Starius

doi:10.16993/tellusa.342

1. Introduction

In a typical 3-dimensional discrete version of an atmospheric model the meshlength Δz in the vertical direction is considerably smaller than the horizontal ones, here both called Δx. The ratio Δx/Δz can, e.g., exceed 100. One way to avoid severe restrictions on the time step, Δt, is to use a horizontally explicit and vertically implicit (HEVI) method for the time integration. This is a considerable simplification, compared to a fully implicit method, since the implicit/vertical part is 1-dimensional and thus a related Jacobian matrix consists of only a narrow dense band.

We mention that in a nonhydrostatic model, see for instance (Giraldo, Kelly and Constantinescu, 2012) and (Staniforth and Wood, 2008), vertically propagating acoustic waves will appear and make the vertical Courant number even greater.

In (Mengaldo et al., 2018), current and emerging time integration strategies for atmospheric models are discussed.

In this paper we develop a family of fourth order two step Runge-Kutta methods with four stages, for an initial value problem for an additively partitioned system of ordinary differential equations y′ = s(y) + f(y), y ∈ R^d, d > 1, where s is the slow and f the fast part of the system. The free parameters in the family are defined by optimising the HEVI-stability, or shorter H-stability, properties related to the test equation y′ = –ik_xy – ik_zy, where k_x and k_z are real numbers and i² = –1, see for instance (Lock, Wood and Weller, 2014). The result, H-stability if |Δtk_x| ≤ 2 with no restriction on k_z, seems to be quite satisfactory. Our methods will not use stage values from previous steps as is done for instance in (Zharovski, Sandu and Zhang, 2015), where there is no emphasis on H-stability but rather on approximation orders and A(α)-stability.

Since it is beyond the scope of this paper to apply our time integration method to realistic atmospheric problems, we will only consider some very simple numerical examples in order to illustrate some properties of our method.

The present paper is organized as follows. In Section 2, we rewrite our additively partitioned system as a partitioned system by replacing y with two new variables p and q such that y = p + q. For this system we define two step Runge-Kutta methods and then by just addition arrive at the corresponding discretisation for y′ = s(y) + f(y). Appropriate notation for the local discretisation errors is also introduced. At the end of the section we define H-stability in detail.

In Section 3, the family of fourth order two step Runge-Kutta methods with four stages are derived in a fairly detailed manner. The free parameters of the family are determined by an attempt to optimise the H-stability properties of the resulting method. This method is given in a direct and simple way, that is without a new kind of Butcher tableau. We note that it is possible to derive fourth order methods in this context with only three stages, but it seems impossible to get good H-stability properties as well.

In Section 4, we compare our two step Runge-Kutta method tsRK4(4,4,4) with the one step method ARS3(4,4,3), by (Ascher, Ruuth and Spiteri, 1997). Both use four stages but have different approximation orders, namely four and three, respectively. According to a convention by (Pareschi and Russo, 2005) the numbers k(ρ, σ, κ) mean: k is the order of the explicit part, ρ is the number of implicit stages, σ is the number of explicit stages, and κ is the order of the whole method. To our knowledge there is no known one step method characterized by 4(4,4,4), which we could use for comparison purposes, instead of ARS3 (4,4,3).

In the beginning of Section 4 the H-stability properties for the two methods are compared. For ARS3(4,4,3) these are significantly less satisfactory than for tsRK4(4,4,4). Further, a comparison is made between the two methods by applying them to two simple test problems. The first problem contains only one time scale and does not require a partitioned method, but is suitable for comparing accuracies and checking approximation orders. The second test problem contains two quite different time scales, of which only the slowest is of interest to us. The fast scale is not resolved and is not allowed to interfere harmfully with the slow part of the numerical solution.

2. Some preliminaries for two step Runge-Kutta methods for partitioned systems of ordinary differential equations

Let us consider an initial value problem for additively partitioned systems of ordinary differential equations

1

y^{'} = s (y) + f (y), for t \geq 0 with y (0) = y_{0} .

In the right hand side of the equation, s and f are given functions from R^d to R^d, where d is an integer with d > 1 and R denotes the set of real numbers. In (1) s(y) can, e.g., represent the horizontal part of a discretisation of an atmospheric model and f(y) the corresponding vertical part. We will treat f(y) implicitly in order to avoid very small time steps and s(y) explicitly to make reasonable computing times possible.

Following (Ascher, Ruuth and Spiteri, 1997) we rewrite the system (1) as a partitioned system by introducing p and q so that p′ = s(y) and q′ = f(y). It follows that p′ + q′ = y′ and we choose p + q = y. Fortunately, p or q do not have to be known, they are only intermediate variables used for determining the coefficients in the Runge-Kutta methods for (1). We write down the partitioned system

2

{\begin{array}{l} p^{'} = s (p + q) \\ q^{'} = f (p + q) \end{array}, where p + q = y,

for later references.

Two step Runge-Kutta methods for (2), using data from the time-levels t_n = nΔt and t_n–₁ = t_n – Δt to evaluate values corresponding to t_n₊₁ = t_n + Δt, can be written as

3

{\begin{array}{l} P_{0} = p_{n - 1}, P_{1} = p_{n} and Q_{0} = q_{n - 1}, Q_{1} = q_{n}, \\ P_{i} = d_{i} p_{n - 1} + (1 - d_{i}) p_{n} + Δ t \sum_{j = 0}^{i - 1} a_{i j} s (P_{j} + Q_{j}), \\ Q_{i} = d_{i} q_{n - 1} + (1 - d_{i}) q_{n} + Δ t \sum_{j = 0}^{i} b_{i j} f (P_{j} + Q_{j}), \end{array}

where i = 2, …, r and r – 1 is the number of stages in the method. Further, p_n₊₁ = P_r and q_n₊₁ = Q_r are the numerical approximations to p(t_n₊₁) and q(t_n₊₁). The second and third rows of (3) are called the explicit and implicit part, respectively.

By addition in (3) and setting Y_i = P_i + Q_i, y_n_–1 = p_n–₁ + q_n–₁ and y_n = p_n + q_n we obtain Runge-Kutta methods for (1)

4

{\begin{array}{l} Y_{0} = y_{n - 1}, Y_{1} = y_{n}, \\ Y_{i} = d_{i} y_{n - 1} + (1 - d_{i}) y_{n} + Δ t \sum_{j = 0}^{i - 1} a_{i j} s (Y_{j}) + Δ t \sum_{j = 0}^{i} b_{i j} f (Y_{j}), \end{array}

where i = 2, …, r and in accordance with the above y_n₊₁ = Y_r is the numerical approximation to y(t_n₊₁). In this paper all stages in (4) will be implicit.

We now consider the local discretisation errors for (3) by first defining P_i and Q_i as

5

{\begin{array}{l} P_{0} = p (t - Δ t), P_{1} = p (t) and Q_{0} = q (t - Δ t), Q_{1} = q (t), \\ P_{i} = d_{i} p (t - Δ t) + (1 - d_{i}) p (t) + Δ t \sum_{j = 0}^{i - 1} a_{i j} s (P_{j} + Q_{j}), \\ Q_{i} = d_{i} q (t - Δ t) + (1 - d_{i}) q (t) + Δ t \sum_{j = 0}^{i} b_{i j} f (P_{j} + Q_{j}), \end{array}

where i = 2, …, r. We note that the corresponding numerical quantities P_i and Q_i only appear before (5) and therefore cannot be confused with P_i and Q_i. The second terms below are the local discretisation errors for each stage

6

{\begin{array}{l} P_{i} = p (t + c_{i} Δ t) + (P_{i} - p (t + c_{i} Δ t)), \\ Q_{i} = q (t + c_{i} Δ t) + (Q_{i} - q (t + c_{i} Δ t)), \end{array}

where i = 0, …, r. The coefficients c_i, i = 2, …, r – 1 are still free but, according to the above, we must set c₀ = –1, c₁ = 0 and c_r = 1. For i = r in (6) we get the local discretisation errors for the method (3), which we want to be of order Δt^κ+1 for some reasonable value of κ.

We now briefly recall the HEVI- or H-stability concept for the Runge-Kutta method (4). Apply the method to the scalar test equation

7

y^{'} = - i k_{x} y - i k_{z} y, where k_{x} and k_{z} are real numbers .

If for a pair (Δtk_x, Δtk_z), all numerical solutions y_n, n = 1, 2, 3, … are bounded, we say that (4) is H-stable for that pair. The set of such pairs is called the H-stability region of (4).

3. Fourth order two step Runge-Kutta methods for partitioned systems of ordinary differential equations

We will now derive fourth order methods with r – 1 = 4 stages, good H-stability properties and finally with b_ii = b, for i = 2, …, 5. For a method (3) to be of order four the local discretisation errors for the method must be one order higher, that is

8

{\begin{array}{l} L_{5}^{p} = P_{5} - p (t + Δ t)) = O (Δ t^{5}), \\ L_{5}^{q} = Q_{5} - q (t + Δ t)) = O (Δ t^{5}) . \end{array}

We investigate the local discretisation errors (6) by using Taylor expansion, first for the case i = 2

9

{\begin{array}{l} \begin{array}{l} P_{2} - p (t + c_{2} Δ t) = Δ t (a_{20} + a_{21} - d_{2} - c_{2}) p^{'} (t) \\ - \frac{1}{2} Δ t^{2} (2 a_{20} - d_{2} + c_{2}^{2}) p^{″} (t) + \frac{1}{6} Δ t^{3} p_{23} p^{‴} (t) \end{array} \\ + O (Δ t^{4}), where p_{23} = 3 a_{20} - d_{2} - c_{2}^{3}, \\ \begin{array}{l} Q_{2} - q (t + c_{2} Δ t) = Δ t (b_{20} + b_{21} + b - d_{2} - c_{2}) q^{'} (t) \\ - \frac{1}{2} Δ t^{2} (2 b_{20} - 2 c_{2} b - d_{2} + c_{2}^{2}) q^{″} (t) + \frac{1}{6} Δ t^{3} q_{23} q^{‴} (t) \end{array} \\ + O (Δ t^{4}), where q_{23} = 3 b_{20} + 3 c_{2}^{2} b - d_{2} - c_{2}^{3} . \end{array}

By requiring that the first two terms in each right hand side of (9) vanish, we obtain

10

\begin{array}{l} P_{2} = p (t + c_{2} Δ t) + \frac{1}{6} Δ t^{3} p_{23} p^{‴} + O (Δ t^{4}) and \\ Q_{2} = q (t + c_{2} Δ t) + \frac{1}{6} Δ t^{3} q_{23} q^{‴} + O (Δ t^{4}) . \end{array}

The conditions for (10) to hold can also be expressed as

11

\begin{array}{l} a_{21} = c_{2} + c_{2}^{2} + a_{20}, d_{2} = c_{2}^{2} + 2 a_{20} and \\ b_{20} = c_{2} b + (d_{2} - c_{2}^{2}) / 2, b_{21} = c_{2} - b_{20} - b + d_{2}, \end{array}

where we have chosen to let c₂, a₂₀ and b(= b₂₂) be free parameters. Later, when values are assigned to them, p₂₃ and q₂₃ in (9) are also determined.

In analogy to (10) we require that

12

\begin{array}{l} P_{3} = p (t + c_{3} Δ t) + \frac{1}{6} Δ t^{3} p_{33} p^{‴} + O (Δ t^{4}) and \\ Q_{3} = q (t + c_{3} Δ t) + \frac{1}{6} Δ t^{3} q_{33} q^{‴} + O (Δ t^{4}), \\ P_{4} = p (t + c_{4} Δ t) + \frac{1}{6} Δ t^{3} p_{43} p^{‴} + O (Δ t^{4}) and \\ Q_{4} = q (t + c_{4} Δ t) + \frac{1}{6} Δ t^{3} q_{43} q^{‴} + O (Δ t^{4}), \end{array}

where p₃₃, q₃₃, p₄₃ and q₄₃ shall be determined such that (8) holds.

We shall now expand the errors in (8) by inserting the stages 2, 3 and 4 in stage 5 and make use of equalities like

\begin{array}{l} s (P_{3} + Q_{3}) = s (p (t + c_{3} Δ t) + q (t + c_{3} Δ t)) \\ + \frac{1}{6} Δ t^{3} \frac{\partial s}{\partial y} (p_{33} p^{‴} + q_{33} q^{‴}) + O (Δ t^{4}), \end{array}

and s(p(t + c₃Δt) + q(t + c₃Δt)) = p′(t + c₃Δt) according to (2). Taylor expansion and requiring that the coefficients of p′, p″, p‴, p⁽⁴⁾ and q′, q″, q‴, q⁽⁴⁾ vanish lead to the two systems of linear equations

13

\begin{array}{l} M (\begin{array}{l} a_{51} \\ a_{52} \\ a_{53} \\ a_{54} \end{array}) = (\begin{array}{l} 1 + d_{5} \\ 1 - d_{5} \\ 1 + d_{5} \\ 1 - d_{5} \end{array}), M (\begin{array}{l} b_{51} \\ b_{52} \\ b_{53} \\ b_{54} \end{array}) = h_{1}, \\ where M = (\begin{array}{l} 1 & 1 & 1 & 1 \\ 0 & 2 c_{2} & 2 c_{3} & 2 c_{4} \\ 0 & 3 c_{2}^{2} & 3 c_{3}^{2} & 3 c_{4}^{2} \\ 0 & 4 c_{2}^{3} & 4 c_{3}^{3} & 4 c_{4}^{3} \end{array}) \end{array}

and h₁ = (1 + d₅ – b, 1 – d₅ – 2b, 1 + d₅ – 3b, 1 – d₅ – 4b)^T. Provided that the two systems in (13) hold true the local discretisation errors can be written as

14

\begin{array}{l} L_{5}^{p} = \frac{1}{6} Δ t^{4} \frac{\partial s}{\partial y} [a_{52} (p_{23} p^{‴} + q_{23} q^{‴}) + a_{53} (p_{33} p^{‴} + q_{33} q^{‴}) \\ + a_{54} (p_{43} p^{‴} + q_{43} q^{‴})] + O (Δ t^{5}), \\ L_{5}^{q} = \frac{1}{6} Δ t^{4} \frac{\partial f}{\partial y} [b_{52} (p_{23} p^{‴} + q_{23} q^{‴}) + b_{53} (p_{33} p^{‴} + q_{33} q^{‴}) \\ + b_{54} (p_{43} p^{‴} + q_{43} q^{‴})] + O (Δ t^{5}) . \end{array}

Since we want the global errors of (3) to be of order four, all terms in (14) of the same order must vanish, which lead to the following equations

15

\begin{array}{l} C (\begin{array}{l} p_{33} \\ p_{43} \end{array}) = - p_{23} (\begin{array}{l} a_{52} \\ b_{52} \end{array}), C (\begin{array}{l} q_{33} \\ q_{43} \end{array}) = - q_{23} (\begin{array}{l} a_{52} \\ b_{52} \end{array}), \\ where C = (\begin{matrix} a_{53} & a_{54} \\ b_{53} & b_{54} \end{matrix}) . \end{array}

We will now give equations for determination of the coefficients in the stages 3 and 4 by Taylor expansion and by requiring (12) to hold. For the explicit part of stage 3 we obtain

16

(\begin{matrix} 1 & 1 & - 1 \\ 0 & 2 c_{2} & 1 \\ 0 & 3 c_{2}^{2} & - 1 \end{matrix}) (\begin{matrix} a_{31} \\ a_{32} \\ d_{3} \end{matrix}) = (\begin{matrix} c_{3} - a_{30} \\ c_{3}^{2} + 2 a_{30} \\ c_{3}^{3} - 3 a_{30} + p_{33} \end{matrix})

and for the implicit part

17

(\begin{matrix} 1 & 1 & 1 \\ - 2 & 0 & 2 c_{2} \\ 3 & 0 & 3 c_{2}^{2} \end{matrix}) (\begin{matrix} b_{30} \\ b_{31} \\ b_{32} \end{matrix}) = (\begin{matrix} c_{3} - b + d_{3} \\ c_{3}^{2} - 2 c_{3} b - d_{3} \\ c_{3}^{3} - 3 c_{3}^{2} b + d_{3} + q_{33} \end{matrix}) .

For stage 4, with a₄₀ = b₄₀ = 0, we obtain in a similar way

18

\begin{array}{l} E (\begin{matrix} a_{41} \\ a_{42} \\ a_{43} \end{matrix}) = (\begin{matrix} c_{4} + d_{4} \\ c_{4}^{2} - d_{4} \\ c_{4}^{3} + d_{4} + p_{43} \end{matrix}), E (\begin{matrix} b_{41} \\ b_{42} \\ b_{43} \end{matrix}) = h_{2}, \\ where E = (\begin{matrix} 1 & 1 & 1 \\ 0 & 2 c_{2} & 2 c_{3} \\ 0 & 3 c_{2}^{2} & 3 c_{3}^{2} \end{matrix}) \end{array}

and h₂ = (c₄ + d₄ – b, c₄² – d₄ – 2c₄b, c₄³ + d₄ – 3c₄²b + q₄₃)^T.

To avoid unnecessarily complex methods we will restrict ourselves to the subclass characterized by

a_{20} = a_{30} = a_{40} = a_{50} = 0, b_{40} = b_{50} = 0, and d_{4} = d_{5} = 0.

The remaining free parameters are c₂, c₃, c₄ and b, which we have used in an attempt to optimise the H-stability properties of our method for (1). Each method occurring in the optimisational process is obtained by using (11), (13) and (15–18). Simple tabulation and subtabulation has given us the values

c_{2} = 0.4, c_{3} = 1.2, c_{4} = 0.5 and b = 0.6,

which correspond to a H-stability region containing the set

19

{(Δ t k_{x}, Δ t k_{z}) : - 2 \leq Δ t k_{x} \leq 2.1 and k_{z} an arbitrary real number} .

For the convenience of potential users, we write down the corresponding method in detail

20

{\begin{array}{l} \begin{array}{l} Y_{2} = d_{2} y_{n - 1} + (1 - d_{2}) y_{n} + Δ t [a_{21} s (y_{n}) + b_{20} f (y_{n - 1}) \\ + b_{21} f (y_{n}) + b_{22} f (Y_{2})], \end{array} \\ Y_{3} = d_{3} y_{n - 1} + (1 - d_{3}) y_{n} + Δ t [a_{31} s (y_{n}) + a_{32} s (Y_{2}) \\ + b_{30} f (y_{n - 1}) + b_{31} f (y_{n})] + Δ t [b_{32} f (Y_{2}) + b_{33} f (Y_{3})], \\ \begin{array}{l} Y_{4} = y_{n} + Δ t [a_{41} s (y_{n}) + \sum_{j = 2}^{3} a_{4 j} s (Y_{j}) + b_{41} f (y_{n}) \\ + \sum_{j = 2}^{4} b_{4 j} f (Y_{j})], \end{array} \\ \begin{array}{l} Y_{5} = y_{n} + Δ t [a_{51} s (y_{n}) + \sum_{j = 2}^{4} a_{5 j} s (Y_{j}) + b_{51} f (y_{n}) \\ + \sum_{j = 2}^{5} b_{5 j} f (Y_{j})], \end{array} \end{array}

where b₂₂ = b₃₃ = b₄₄ = b₅₅ = b = 3/5, and further

21

\begin{array}{l} Y_{2} : {\begin{array}{l} d_{2} = 4 / 25, a_{21} = 14 / 25, \\ b_{20} = 6 / 25, b_{21} = - 7 / 25, \end{array} \\ Y_{3} : {\begin{array}{l} d_{3} = 11 / 25, a_{31} = 39 / 100, a_{32} = 5 / 4, \\ b_{30} = 222 / 175, b_{31} = - 57 / 20, b_{32} = 367 / 140, \end{array} \\ Y_{4} : {\begin{array}{l} a_{41} = 49 / 288, a_{42} = 65 / 192, a_{43} = - 5 / 576, \\ b_{41} = 371 / 1440, b_{42} = - 61 / 192, b_{43} = - 23 / 576, \end{array} \\ Y_{5} : {\begin{array}{l} a_{51} = 5 / 24, a_{52} = - 25 / 48, a_{53} = 25 / 336, a_{54} = 26 / 21, \\ b_{51} = 7 / 120, b_{52} = 65 / 48, b_{53} = - 65 / 336, b_{54} = - 86 / 105. \end{array} \end{array}

We recall that y_n₊₁ = Y₅ is the numerical approximation to y(t_n₊₁).

If the functions s or f in (1) explicitly depend on the time t, then s(t+c_i Δt, Y_i) or f(t+c_i Δt, Y_i) shall be used in (4) instead of s(Y_i) or f(Y_i).

4. Numerical experiments

The main purpose of this section is to compare our method tsRK4(4,4,4) with the method ARS3 (4,4,3) by (Ascher, Ruuth and Spiteri, 1997). Both methods require about the same amount of computational work per time step, since they use the same number of stages.

For tsRK4(4,4,4) the first step, from 0 to Δt, will be made by making two steps by ARS3(4,4,3), with time step Δt/2. The error at t = Δt will be of order four, and so will the global error.

The H-stability properties of our method are given in (19) and can be contracted to the simple form |Δtk_x| ≤ 2. For ARS3 (4,4,3) with k_x < 0 the situation is more complicated. Let λ be the amplification factor, then we can summarize the stability properties for ARS3(4,4,3) as

0 \leq Δ t k_{x} \leq 1.5, with | λ | \leq 1 or - 1.3 \leq Δ t k_{x} < 0, with | λ | \leq 1.003.

This is definitively less satisfactory than for the method tsRK4(4,4,4). Note that if the point (Δtk_x, Δtk_z) is H-stable so is (–Δtk_x, –Δtk_z), provided the coefficients in the numerical method are all real.

We will now use the two numerical methods on two simple test problems. The first one contains only one time scale and thus does not require a partitioned method, but is suitable for comparing accuracies and checking approximation orders. Our first problem, with oscillating solution, is

22

{(\begin{matrix} u \\ v \end{matrix})}^{'} = (\begin{matrix} 0 & - a (t) \\ a (t) & 0 \end{matrix}) (\begin{matrix} u \\ v \end{matrix}), (\begin{matrix} u (0) \\ v (0) \end{matrix}) = (\begin{matrix} 1 \\ 0 \end{matrix}),

where a(t) = 1–1/(1 + t)². Denote the right hand side of the system by r_s, which we split like r_s = α ∙ r_s+ (1 – α) ∙ r_s and treat the first term explicitly and the second implicitly. We will choose α =2(1 – α), which gives α = 2/3. By defining y = u+iv, u and v real, then (22) can be rewritten as y′ = ia(t) y, y(0) = 1 and we easily get the solution $y (t) = \cos (t \frac{t}{1 + t}) + i \sin (t \frac{t}{1 + t})$ . We mention that the real and complex form of the problem give the same numerical solution.

We note that the solution is approximately 2π periodic for t ≫ 1. Let us define the time steps as Δt = 2π/m and let the integration intervals be [0, T], where T = 2πN and N is the number of approximate periods we integrate over. Finally, we let m = 5, 10, 20, 40 and N = 5, 10, 20 and compute the errors |y_n – y(2πN)| with n = mN, which since |y(t)| = 1 are also the relative errors, see Table 1 below.

Table 1

Comparison between the methods tsRK4(4,4,4) and ARS3(4,4,3) for the problem (22), with α = 2/3. Time steps Δt = 2π/m and integration intervals [0, 2πN].

m	N	error tsRK4(4,4,4)	error ARS3(4,4,3)
5	5	8.7501e-02	6.6770e-01
10	5	6.4467e-03	1.2622e-01
20	5	4.2897e-04	1.6895e-02
40	5	2.7854e-05	2.1340e-03
5	10	1.8045e-01	9.1760e-01
10	10	1.3314e-02	2.4161e-01
20	10	8.7283e-04	3.4335e-02
40	10	5.5842e-05	4.3733e-03
5	20	3.5877e-01	1.0068e+00
10	20	2.7080e-02	4.2989e-01
20	20	1.7635e-03	6.8352e-02
40	20	1.1197e-04	8.8442e-03

The difference in accuracy is obvious even for small resolutions, m = 5 or 10 say. By for instance using the last two rows in Table 1 we have 1.7635/0.11197 = 2^{3.9772 …} ≈ 2⁴ and 6.8352/0.88442 = 2^{2.9501 …} ≈ 2³. The quotients for tsRK4(4,4,4) and ARS3(4,4,3) clearly indicate approximation orders four and three, respectively.

We now consider our second test problem, which contains two quite different time scales. The fastest of these is of no interest to us, and therefore we will not resolve it and rather not allow it to influence the slow part of the numerical solution in a harmful way.

Let us consider the initial value problem

23

u^{″} - i (ω + 1) u^{'} - ω u = 0, u (0) = 1, u^{'} (0) = i (1 + ϵ),

where ω ≫ 1 and ε > 0 are parameters. The solution is given by

u (t) = (1 - ϵ / (ω - 1)) e^{i t} + ϵ / (ω - 1) e^{i ω t} .

By introducing v = u′ and the vector y^T = (u v) we obtain the system

24

\begin{array}{l} y^{'} = (\begin{matrix} 0 & 1 \\ ω & i (ω + 1) \end{matrix}) y = (\begin{matrix} 0 & 1 \\ 0 & i \end{matrix}) y \\ + (\begin{matrix} 0 & 0 \\ ω & i ω \end{matrix}) y, y (0) = (\begin{matrix} 1 \\ i (1 + ϵ) \end{matrix}), \end{array}

where the first term in the splitted right hand side will be treated explicitly and the second implicitly. The complex system (24) can of course be rewritten as a real 4 × 4 system of the form (1).

Since we do not resolve the fastest time scale, the magnitude of the last term in u(t) = (1 – ε/(ω – 1))e^it + ε/(ω – 1)eⁱω^t, that is

| ϵ / (ω - 1) e^{i ω t} | \leq 5.1 \cdot 10^{- 4},

will be a limit for the maximal reliable accuracy, which seems to be attainable for tsRK4(4,4,4) but not for ARS3(4,4,3), see Table 2. No surprise that the accuracy obtained by tsRK4(4,4,4) increases much faster for decreasing Δt, than is the case for ARS3(4,4,3).

Table 2

Comparison between the methods tsRK4(4,4,4) and ARS3(4,4,3) for the problem (24), with ω = 100, ε = 0.05, time steps Δt =2π/m, integration intervals [0, 2πN] and error = |u_n – u(2πN)|, with n = mN.

m	N	error tsRK4(4,4,4)	error ARS3(4,4,3)
10	10	2.2533e-01	6.7569e-01
20	10	1.5140e-02	1.1932e-01
40	10	1.0841e-03	1.5515e-02
80	10	4.7040e-04	2.2383e-03
160	10	3.3149e-04	8.3100e-04
320	10	5.6479e-04	8.8426e-04
10	20	4.1622e-01	9.3054e-01
20	20	3.0132e-02	2.2622e-01
40	20	2.0105e-03	3.1081e-02
80	20	4.7033e-04	4.1364e-03
160	20	3.3283e-04	1.0762e-03
320	20	5.6482e-04	9.1561e-04

5. Summary and conclusions

Our two step Runge-Kutta method tsRK4(4,4,4) and the one step method ARS3(4,4,3), by (Ascher, Ruuth and Spiteri, 1997), use the same number of stages and therefore require about the same computational work per time step. The results for the problem (22) in Table 1 with only one time scale, and the problem (24) in Table 2 with two quite different scales both clearly show that tsRK4(4,4,4) is much more accurate than ARS3 (4,4,3), for different resolutions and integration intervals.

The H-stability, related to the test equation y′ = –ik_xy – ik_zy, can be summarized for tsRK4(4,4,4) by |Δtk_x| ≤ 2 and no condition on Δk_z. H-stability for ARS3(4,4,3) holds for 0 ≤ Δtk_x ≤ 1.5 and no condition on k_z, but for the case k_x < 0 the situation is unsatisfactory, for details see the beginning of Section 4.

The method tsRK4(4,4,4) is likely to be suitable for spatial discretisations of order four, e.g., for finite volume methods on cubed sphere or icosahedral grids according to for instance (Ullrich, Jablonowski and Leer, 2010) and Pudykiewics (2011) and for centred finite differences on an Equator-Pole grid system according to (Starius, 2018; Starius, 2020). By using an Equator-Pole grid system the spatial discretisations can easily and with low extra cost attain the orders 6, 8 or 10, say. To find time integrators of the same orders κ, with κ or perhaps κ + 1 stages and with H-stability properties, such that the time step Δt can be chosen mainly by accuracy considerartions, that is by balancing spatial and temporal discretisation errors, is probably not possible for κ much greater than four. Thus, it is reasonable to consider methods with spatial order greater than the temporal one. Let us, as a thought experiment, consider the method we get by decreasing the spatial order to the temporal one and assume that it is well balanced. If we now increase the spatial order to its original value, balancing will lead to an increase of the spatial steps Δx and the quotient Δt/Δx will decrease. This means that spatial order greater than the temporal one may imply reduced H-stability requirements on the time integrator, which can be adventageous.

Acknowledgements

I would like to thank the reviewer of my paper Starius (2020), who drew my attention to HEVI methods for 3-dimensional atmospheric models, which influenced me to start working on this paper.

Competing Interests

The author has no competing interests to declare.