Optimal Parameters for Third Order Runge–Kutta Exponential Integrators for Convection–Diffusion Problems

Abstract

We re-visit the semi-Lagrangian Runge–Kutta exponential integrators presented in Celledoni et al. (J Sci Comput 41:139–164, 2009) for Convection-dominated convection-diffusion problems. Third order accurate methods of this class consist of subsets of coefficients of commutator-free exponential integrators including some free parameters. We consider an optimal choice of parameters for the third order accurate methods based on the least-squares methods applied on fourth order commutator-free conditions.

Appendices

Appendices

1.1 A Deriving the Order Conditions

We start with a pair of partitioned Runge–Kutta (PRK) methods (3) of overall temporal order p. We seek for coefficients and , $j,k = 1,\ldots , s,\, l = 1,\ldots ,J$ of commutator-free (CF) exponential integrators [7] defined such that (4) hold. For the sake of convenience we shall interpret the PRK-CF methods ((5) and (6)) such that at each step the numerical solution $y_{n+1}$ is expressed as an extra stage value $Y_{s+1}$. We have (Tables 1, 2, 3, 4)

Table 1 Order conditions for orders $p = 1, 2$, showing coefficients of essential elementary differentials compared

Table 2 Order conditions for order $p = 3$, showing coefficients of essential elementary differentials compared

where

The order conditions are obtained by requiring from the Taylor expansions of (2) and (5) that

(16)

where $p = 1,2,3,\ldots$ is the desired order of convergence of the method.

Taking the $q^\text {th}$ derivative with respect to h of the exact solution to (2) and of the stage values of the numerical solution we obtain the recursive formulas

(17)

(18)

where $\displaystyle \sigma _{k,r}^q {:}{=} (q-k){q\atopwithdelims ()k}{q-k-1\atopwithdelims ()r}$.

We will often simplify higher order conditions using conditions of lower order whenever necessary. The computation of the derivatives of $Y_i$ involves computation of derivatives of $\varphi _i$ and $\varphi _j^{-1}$.

Now let us consider the matrix-valued functions,

We denote by either of or for $l = 0,1,\ldots ,J-1$, and define

Depending on the choice of or , we have $\psi _i = \varphi _i$ or $\varphi _i^{-1}$ respectively.

Recall the definitions of the Ad and dexp operators as follows:

where $u, w, \psi$ are matrices of compatible dimensions, and $\text {ad}_w(u) = [w,u] {:}{=} wu - uw$ denotes the commutator of w and u.

Then

So we can write $\dot{\psi }_i = S_i(h)\psi _i$, where

and as a direct consequence we have

(19)

Analogously for $\psi _i^l$ we get that $\dot{\psi }_i^l = S_i^l(h)\psi _i^l$, where

Now we have the following proposition for finding the derivatives of $S_i(h)$:

Proposition 1

Given that $Z^0 = Z^0(h)$ and $W=W(h)$ matrix-valued differentiable functions of h with W(h) invertible for all h sufficiently small, then

(20)

with

$$\begin{aligned} Z^r = \left[ W^{-1}\dot{W}, Z^{r-1} \right] + \dot{Z}^{r-1}. \end{aligned}$$

(21)

The proof is by induction.

Differentiation of (21) gives

$$\begin{aligned} \dot{Z}^r = \left[ -(W^{-1}\dot{W})^2 + W^{-1}\ddot{W},Z^{r-1} \right] +\left[ W^{-1}\dot{W}, \dot{Z}^{r-1} \right] + \ddot{Z}^{r-1}. \end{aligned}$$

(22)

Using (21) and (22) and assuming $W(0) = I$ we obtain the following:

(23)

Further assuming that $Z^0 = \mathrm {dexp}_{-B}(B)$ for some differentiable matrix-valued function $B = B(h)$ we obtain the derivatives at $h=0$ as follows:

$$\begin{aligned} Z^0(0)&= \dot{B}(0), \nonumber \\ \dot{Z}^0(0)&= \ddot{B}(0), \nonumber \\ \ddot{Z}^0(0)&= \dddot{B}(0) - \frac{1}{2}\left[ \dot{B}(0),\ddot{B}(0) \right] , \nonumber \\ \dddot{Z}^0(0)&= B^{iv}(0) - \left[ \dot{B}(0),\dddot{B}(0) \right] + \frac{1}{2}\left[ \dot{B}(0),\left[ \dot{B}(0),\ddot{B}(0) \right] \right] , \nonumber \\&\vdots \end{aligned}$$

(24)

For clarity of notation we define the following matrix derivatives:

Definition 1

Suppose that is a k-times differentiable matrix function of , where $k = 1,2,3,\ldots$ Let . Then $C'(y)(u)$, $C''(y)(u,v)$, $C'''(y)(u,v,w)$ etc, are defined as

For the sake of simplicity we shall henceforth write $C(y_n),C'(y_n),C''(y_n),...$ as $C,C',C'',\ldots$ to denote the convection matrix and its derivatives at $y_n$. In our calculations, we exploit the fact that each one of these operators is linear in all its arguments.

1.2 A.1 Derivatives of $S_i$

The derivatives of $S_i$ require that we first compute the derivatives of or at $h=0$. Direct computations yields

We obtain derivatives of by changing the signs in the expressions above and replacing the by .

We can now obtain derivatives of $S_i$ at $h=0$ and subsequently the derivatives $\psi _i$, by applying in a recursive manner the steps in (20)-(24) with $W = \psi _i^l$ and . We shall exploit also the fact that , for any $l,r = 0,\ldots ,J-1$. We obtain

(25)

At this point forward, we shall write , simply as without emphasizing evaluation at $h=0$ (unless where that give rise ambiguities).

The third derivative of $S_i$ at $h=0$ contains many more terms compared to the first two derivatives:

(26)

1.3 A.2 Derivatives of y(t)

We note that

We also assume that $y(t_n) = y_n$. Starting with $\dot{y} = (C(y)+A)y$, at $t=t_n$ we get

$$\begin{aligned} \dot{y}_n = (C + A)y_n. \end{aligned}$$

(27)

The second derivative:

$$\begin{aligned} \ddot{y} = C'(y)(\dot{y})y + (C(y)+A)\dot{y} \end{aligned}$$

at $t = t_n$ gives

$$\begin{aligned} \ddot{y}_n = C'((C+A)y_n)y_n + (C+A)^2y_n. \end{aligned}$$

(28)

The third derivative:

$$\begin{aligned} \dddot{y} = C''(y)(\dot{y},\dot{y})y + C'(y)(\ddot{y})y + 2C'(y)(\dot{y})\dot{y} + (C(y)+A)\ddot{y}, \end{aligned}$$

at $t = t_n$ gives

$$\begin{aligned} \dddot{y}_n&= C''((C+A)y_n,(C+A)y_n)y_n + C'(C'((C+A)y_n)y_n)y_n \nonumber \\&\quad + C'((C+A)^2y_n)y_n + 2C'((C+A)y_n)(C+A)y_n \nonumber \\&\quad + (C+A)C'((C+A)y_n)y_n + (C+A)^3y_n \end{aligned}$$

(29)

The fourth derivative:

$$\begin{aligned} y^{iv}&= C'''(y)(\dot{y},\dot{y},\dot{y})y + C''(y)(\ddot{y},\dot{y})y + 2C''(y)(\dot{y},\ddot{y})y\\&\quad + 3C''(y)(\dot{y},\dot{y})\dot{y} + C'(y)(\dddot{y})y + 3C'(y)(\ddot{y})\dot{y} + 3C'(y)(\dot{y})\ddot{y}\\&\quad + (C(y)+A)\dddot{y} \end{aligned}$$

at $t = t_n$ gives

$$\begin{aligned} y_n^{iv}&= C'''(\dot{y}_n,\dot{y}_n,\dot{y}_n)y_n + C''(\ddot{y}_n,\dot{y}_n)y_n + 2C''(\dot{y}_n,\ddot{y}_n)y_n \nonumber \\&\quad + 3C''(\dot{y}_n,\dot{y}_n)\dot{y}_n + C'(\dddot{y}_n)y_n + 3C'(\ddot{y}_n)\dot{y}_n + 3C'(\dot{y}_n)\ddot{y}_n + (C+A)\dddot{y}_n \end{aligned}$$

(30)

1.4 A.3 Derivatives of and $Y_i$

We use the formula (19) together with the derivatives of $S_i$ to obtain derivatives of $\psi _i$ ($= \varphi _i$ or $\varphi _i^{-1}$) at $h=0$. These are then substituted in (18) to obtain derivatives of $Y_i$ at $h=0$.

Clearly from the definition, we see that $\psi _i(0)$, so that

$$\begin{aligned} \varphi _i(0) = \varphi _i^{-1}(0) = I \quad \text { and } Y_i(0) = y_n \end{aligned}$$

where I is the identity matrix.

1.4.1 A.A.3.1. First Derivatives

Therefore

From (18)

Substituting for $\dot{\varphi }_i$, $\varphi _i$ and $\varphi _i^{-1}$ we get

(31)

Comparing the coefficients of all the elementary differentials in (31) (for $i=s+1$) with (27), we obtain the first order conditions:

$$\begin{aligned} \sum _k\hat{b}_k = 1\quad \text { and } \sum _jb_j = 1 \end{aligned}$$

(32)

Thefore we can require that both RK methods and $\mathcal {A}$ are consistent, so that

(33)

where $\hat{c}_{s+1} = c_{s+1} = 1$.

1.4.2 A.A.3.2. Second Derivatives

So

From (18)

and we get

(34)

Comparing (34) and (28) for $i=s+1$ we obtain the second order conditions:

$$\begin{aligned} 2\sum _k\hat{b}_k\hat{c}_k&= 1,&2\sum _k\hat{b}_kc_k = 1, \nonumber \\ 2\sum _jb_j\hat{c}_j&= 1,&2\sum _jb_jc_j = 1. \end{aligned}$$

(35)

1.4.3 A.A.3.3. Third Derivatives

Table 3 Derivatives of $\varphi _i$ and $\varphi _i^{-1}$ at $h=0$

The complete expressions for $\dddot{\varphi }_i(0)$ and $\dddot{\varphi _i^{-1}}(0)$ are reported in Table 3. From (18)

We obtain

(36)

Comparing (36) with (29) for $i=s+1$ we obtain the third order conditions reported in Table 2.

1.4.4 A.A.3.4. Fourth Derivatives

Calculations leading to the fourth derivative of $Y_i$ at $h=0$ though much more tedious, follow similar procedures as for the lower order derivatives. Therefore we are going to spare the reader some of the details.

$$\begin{aligned} \psi _i^{iv}(0)&= \dddot{S}_i(0)\psi _i(0) + 3\ddot{S}_i(0)\dot{\psi }_i(0) + 3\dot{S}_i(0)\ddot{\psi }_i(0) + S_i(0)\dddot{\psi }_i(0)\\&= \dddot{S}_i(0) + 3\ddot{S}_i(0)S_i(0) + 3\dot{S}_i(0)^2 + 3\dot{S}_i(0)S_i(0)^2\\&\quad + S_i(0)\ddot{S}_i(0) + 2S_i(0)\dot{S}_i(0)S_i(0) + S_i(0)^2\dot{S}(0) + S_i(0)^4. \end{aligned}$$

From (18) we get

The resulting fourth order conditions are list in Table 4.

Table 4 Order conditions for order $p = 4$, showing coefficients of essential elementary differentials compared