Coupled Multirate Infinitesimal GARK Schemes for Stiff Systems with Multiple Time Scales

TL;DR

This paper introduces coupled implicit multirate GARK schemes to improve stability when integrating stiff systems with multiple time scales, extending previous decoupled methods and validating the approach with numerical experiments.

Contribution

It develops coupled implicit MRI-GARK schemes for stiff multiscale systems, enhancing stability over previous decoupled methods.

Findings

Coupled schemes show improved stability for stiff systems.

Order conditions and stability analysis are derived for the new methods.

Numerical experiments validate the theoretical advantages of coupled methods.

Abstract

Traditional time discretization methods use a single timestep for the entire system of interest and can perform poorly when the dynamics of the system exhibits a wide range of time scales. Multirate infinitesimal step (MIS) methods (Knoth and Wolke, 1998) offer an elegant and flexible approach to efficiently integrate such systems. The slow components are discretized by a Runge-Kutta method, and the fast components are resolved by solving modified fast differential equations. Sandu (2018) developed the Multirate Infinitesimal General-structure Additive Runge-Kutta (MRI-GARK) family of methods that includes traditional MIS schemes as a subset. The MRI-GARK framework allowed the construction of the first fourth order MIS schemes. This framework also enabled the introduction of implicit methods, which are decoupled in the sense that any implicitness lies entirely within the fast or slow integrations. It was shown by Sandu that the stability of decoupled implicit MRI-GARK methods has limitations when both the fast and slow components are stiff and interact strongly. This work extends the MRI-GARK framework by introducing coupled implicit methods to solve stiff multiscale systems. The coupled approach has the potential to considerably improve the overall stability of the scheme, at the price of requiring implicit stage calculations over the entire system. Two coupling strategies are considered. The first computes coupled Runge-Kutta stages before solving a single differential equation to refine the fast solution. The second alternates between computing coupled Runge-Kutta stages and solving fast differential equations. We derive order conditions and perform the stability analysis for both strategies. The new coupled methods offer improved stability compared to the decoupled MRI-GARK schemes. The theoretical properties of the new methods are validated with numerical experiments.

Figures

Equations

$$y^{\prime }=f(t,y)=f^{\left\{\mathfrak{f}\right\}}(t,y)+f^{\left\{\mathfrak{s}\right\}}(t,y),y(t_0)=y_0\in \text{mathbbm}{R}^d,$$

$$Y_1=y_n$$

$$\{\begin{array}{rl}{v}_i(0)& =Y_i,\\ T_i& =t_n+\text{mathbase}{c}_i^{\left\{\mathfrak{s}\right\}}H,\\ v_i^{\prime }& ={\Delta \text{mathbase}c}_i^{\left\{\mathfrak{s}\right\}}{f}^{\left\{\mathfrak{f}\right\}}\left(T_i+{\Delta \text{mathbase}c}_i^{\left\{\mathfrak{s}\right\}}\theta ,v_i\right)+\sum _{j=1}^{i+1}{\gamma }_{i,j}\left(\frac{\theta }{H}\right)f^{\left\{\mathfrak{s}\right\}}\left(T_j,Y_j\right),\\ & \text{for }\theta \in [0,H],\\ Y_{i+1}& =v_i(H),i=1,\dots ,s^{\left\{\mathfrak{s}\right\}},\end{array}$$

$$y_{n+1}=Y_{s^{\left\{\mathfrak{s}\right\}}+1},$$

$$\begin{array}[]{c|c}\mathbase{c}^{\left\{\mathfrak{s}\right\}}&\mathbase{A}^{\left\{\mathfrak{s},\mathfrak{s}\right\}}\\ \hline\cr&\mathbase{b}^{\left\{\mathfrak{s}\right\}T}\\ \hline\cr&\widehat{\mathbase{b}}^{\left\{\mathfrak{s}\right\}T}\end{array}$$

$$Y_i=y_n+H\sum _{j=1}^{s^{\left\{\mathfrak{s}\right\}}}\text{mathbase}{a}_{i,j}^{\left\{\mathfrak{s}\right\}}{f}_j,i=1,\dots ,s^{\left\{\mathfrak{s}\right\}},$$

$$\{\begin{array}{rl}v(0)& =y_n,\\ v^{\prime }& =f^{\left\{\mathfrak{f}\right\}}(t_n+\theta ,v)+\sum _{j=1}^{s^{\left\{\mathfrak{s}\right\}}}{\gamma }_j\left(\frac{\theta }{H}\right)f_j^{\left\{\mathfrak{s}\right\}},\text{for }\theta \in [0,H],\\ y_{n+1}& =v(H),\end{array}$$

$${\gamma }_i(t):=\sum _{k\ge 0}{\gamma }_i^k{t}^k,{\tilde{\gamma }}_i(t):={\int }_0^t{\gamma }_i(\tau )\text{differential}\tau =\sum _{k\ge 0}{\gamma }_i^k\frac{t^{k+1}}{k+1},{\overline{\gamma }}_i:={\tilde{\gamma }}_i(1).$$

$$\{\begin{array}{rl}\hat{v}(0)& =y_n,\\ {\hat{v}}^{\prime }& =f^{\left\{\mathfrak{f}\right\}}\left(t_n+\theta ,\hat{v}\right)+\sum _{j=1}^{s^{\left\{\mathfrak{s}\right\}}}{\hat{\gamma }}_j\left(\frac{\theta }{H}\right)f_j^{\left\{\mathfrak{s}\right\}},\text{for }\theta \in [0,H],\\ {\hat{y}}_{n+1}& =\hat{v}(H),\end{array}$$

$$y_{n+1}=y_n+{\int }_0^H\sum _{j=1}^{s^{\left\{\mathfrak{s}\right\}}}{\gamma }_j\left(\frac{\theta }{H}\right)f_j^{\left\{\mathfrak{s}\right\}}\text{differential}\theta =y_n+H\sum _{j=1}^{s^{\left\{\mathfrak{s}\right\}}}{\overline{\gamma }}_j{f}_j^{\left\{\mathfrak{s}\right\}}.$$

$$\text{mathbase}{b}^{\left\{\mathfrak{s}\right\}}=\overline{\gamma }.$$

$${\left[\begin{array}{c}{y}^{\left\{\mathfrak{f}\right\}}\\ y^{\left\{\mathfrak{s}\right\}}\end{array}\right]}^{\prime }=\left[\begin{array}{c}{f}^{\left\{\mathfrak{f}\right\}}\left(t,y^{\left\{\mathfrak{f}\right\}},y^{\left\{\mathfrak{s}\right\}}\right)\\ f^{\left\{\mathfrak{s}\right\}}\left(t,y^{\left\{\mathfrak{f}\right\}},y^{\left\{\mathfrak{s}\right\}}\right)\end{array}\right]=\left[\begin{array}{c}{f}^{\left\{\mathfrak{f}\right\}}\left(t,y^{\left\{\mathfrak{f}\right\}},y^{\left\{\mathfrak{s}\right\}}\right)\\ 0\end{array}\right]+\left[\begin{array}{c}0\\ f^{\left\{\mathfrak{s}\right\}}\left(t,y^{\left\{\mathfrak{f}\right\}},y^{\left\{\mathfrak{s}\right\}}\right)\end{array}\right].$$

$$\left[\begin{array}{c}{Y}_i^{\left\{\mathfrak{f}\right\}}\\ Y_i^{\left\{\mathfrak{s}\right\}}\end{array}\right]=\left[\begin{array}{c}{y}_n^{\left\{\mathfrak{f}\right\}}+H{\sum }_{j=1}^{s^{\left\{\mathfrak{s}\right\}}}\text{mathbase}{a}_{i,j}^{\left\{\mathfrak{s}\right\}}{f}_j^{\left\{\mathfrak{f}\right\}}\\ y_n^{\left\{\mathfrak{s}\right\}}+H{\sum }_{j=1}^{s^{\left\{\mathfrak{s}\right\}}}\text{mathbase}{a}_{i,j}^{\left\{\mathfrak{s}\right\}}{f}_j^{\left\{\mathfrak{s}\right\}}\end{array}\right],$$

$$\{\begin{array}{rl}{v}^{\left\{\mathfrak{f}\right\}}(0)& =y_n^{\left\{\mathfrak{f}\right\}},\\ {v^{\left\{\mathfrak{f}\right\}}}^{\prime }& =f^{\left\{\mathfrak{f}\right\}}\left(t_n+\theta ,v,y_n^{\left\{\mathfrak{s}\right\}}+H\sum _{j=1}^{s^{\left\{\mathfrak{s}\right\}}}{\tilde{\gamma }}_j\left(\frac{\theta }{H}\right)f_j^{\left\{\mathfrak{s}\right\}}\right),\text{for }\theta \in [0,H]\\ \left[\begin{array}{c}{y}_{n+1}^{\left\{\mathfrak{f}\right\}}\\ y_{n+1}^{\left\{\mathfrak{s}\right\}}\end{array}\right]& =\left[\begin{array}{c}{v}^{\left\{\mathfrak{f}\right\}}(H)\\ y_n^{\left\{\mathfrak{s}\right\}}+H{\sum }_{j=1}^{s^{\left\{\mathfrak{s}\right\}}}\text{mathbase}{b}_j^{\left\{\mathfrak{s}\right\}}{f}_j^{\left\{\mathfrak{s}\right\}}\end{array}\right],\end{array}$$

$$Y_i^{\left\{\mathfrak{f},\right\}}$$

$$\displaystyle=y_{n}+H\sum_{j=1}^{s^{\left\{\mathfrak{f}\right\}}}a^{\left\{\mathfrak{f},\mathfrak{f}\right\}}_{i,j}\,f^{\left\{\mathfrak{f},\right\}}_{j}+H\sum_{\ell=1}^{s^{\left\{\mathfrak{s}\right\}}}\left(\sum_{j=1}^{s^{\left\{\mathfrak{f}\right\}}}a^{\left\{\mathfrak{f},\mathfrak{f}\right\}}_{i,j}\,\gamma_{\ell}\big{(}c^{\left\{\mathfrak{f}\right\}}_{j}\big{)}\right)f^{\left\{\mathfrak{s}\right\}}_{\ell},$$

$$=y_n+H\sum _{j=1}^{s^{\left\{\mathfrak{f}\right\}}}{a}_{i,j}^{\left\{\mathfrak{f},\mathfrak{f}\right\}}{f}_j^{\left\{\mathfrak{f},\right\}}+H\sum _{j=1}^{s^{\left\{\mathfrak{s}\right\}}}{\left(\sum _{k\ge 0}{A}^{\left\{\mathfrak{f},\mathfrak{f}\right\}}{c}^{\left\{\mathfrak{f}\right\}\times k}{{\gamma }^k}^T\right)}_{i,j}{f}_j^{\left\{\mathfrak{s}\right\}},$$

$$y_{n+1}$$

$$=y_n+H\sum _{j=1}^{s^{\left\{\mathfrak{f}\right\}}}{b}_j^{\left\{\mathfrak{f}\right\}}{f}_j^{\left\{\mathfrak{f},\right\}}+H\sum _{j=1}^{s^{\left\{\mathfrak{s}\right\}}}\left(\sum _{k\ge 0}\frac{1}{k+1}{\gamma }_j^k\right)f_j^{\left\{\mathfrak{s}\right\}},$$

$$=y_n+H\sum _{j=1}^{s^{\left\{\mathfrak{f}\right\}}}{b}_j^{\left\{\mathfrak{f}\right\}}{f}_j^{\left\{\mathfrak{f},\right\}}+H\sum _{j=1}^{s^{\left\{\mathfrak{s}\right\}}}{\overline{\gamma }}_j{f}_j^{\left\{\mathfrak{s}\right\}},$$

$$=y_n+H\sum _{j=1}^{s^{\left\{\mathfrak{f}\right\}}}{b}_j^{\left\{\mathfrak{f}\right\}}{f}_j^{\left\{\mathfrak{f},\right\}}+H\sum _{j=1}^{s^{\left\{\mathfrak{s}\right\}}}\text{mathbase}{b}_j^{\left\{\mathfrak{s}\right\}}{f}_j^{\left\{\mathfrak{s}\right\}}.$$

$$\begin{array}[]{c|cc|c|c}\mathbase{c}^{\left\{\mathfrak{s}\right\}}&\mathbase{A}^{\left\{\mathfrak{s},\mathfrak{s}\right\}}&0&\mathbase{A}^{\left\{\mathfrak{s},\mathfrak{s}\right\}}&\mathbase{c}^{\left\{\mathfrak{s}\right\}}\\ c^{\left\{\mathfrak{f}\right\}}&0&A^{\left\{\mathfrak{f},\mathfrak{f}\right\}}&\sum_{k\geq 0}A^{\left\{\mathfrak{f},\mathfrak{f}\right\}}\,c^{\left\{\mathfrak{f}\right\}\times k}\,{\gamma^{k}}^{T}&c^{\left\{\mathfrak{f},\mathfrak{s}\right\}}\\ \hline\cr\mathbase{c}^{\left\{\mathfrak{s}\right\}}&\mathbase{A}^{\left\{\mathfrak{s},\mathfrak{s}\right\}}&0&\mathbase{A}^{\left\{\mathfrak{s},\mathfrak{s}\right\}}&\mathbase{c}^{\left\{\mathfrak{s}\right\}}\\ \hline\cr&0&b^{\left\{\mathfrak{f}\right\}T}&\mathbase{b}^{\left\{\mathfrak{s}\right\}T}\end{array},$$

$${\mathbf{c}}^{\left\{\mathfrak{s},\mathfrak{f}\right\}}=\mathbf{c}\text{SS}\equiv \text{mathbase}{c}^{\left\{\mathfrak{s}\right\}}\text{and}{\mathbf{c}}^{\left\{\mathfrak{f},\mathfrak{f}\right\}}={\mathbf{c}}^{\left\{\mathfrak{f},\mathfrak{s}\right\}}$$

$${{\gamma }^0}^T\text{mathbbm}{1}^{\left\{\mathfrak{s}\right\}}=1\text{and}{{\gamma }^k}^T\text{mathbbm}{1}^{\left\{\mathfrak{s}\right\}}=0\forall k\ge 1.$$

$$c^{\left\{\mathfrak{f}\right\}}$$

$${\mathfrak{t}}_k:=[\underset{k\text{ times}}{\underbrace{\tau ,\dots ,\tau }}]\in T,$$

$${\zeta }_k$$

$${\omega }_k$$

$${\xi }_k$$

$$\frac{1}{6}$$

$$\frac{1}{8}$$

$$\frac{1}{12}$$

$$\frac{1}{24}$$

$$\frac{1}{6}={\mathbf{b}}^{\left\{\mathfrak{f}\right\}T}{\mathbf{A}}^{\left\{\mathfrak{f},\mathfrak{s}\right\}}{\mathbf{c}}^{\left\{\mathfrak{s}\right\}}=\sum _{k\ge 0}{\zeta }_k{{\gamma }^k}^T\text{mathbase}{c}^{\left\{\mathfrak{s}\right\}}.$$

$$\frac{1}{6}={\mathbf{b}}^{\left\{\mathfrak{s}\right\}T}{\mathbf{A}}^{\left\{\mathfrak{s},\mathfrak{f}\right\}}{\mathbf{c}}^{\left\{\mathfrak{f}\right\}}=\text{mathbase}{b}^{\left\{\mathfrak{s}\right\}T}\text{mathbase}{A}^{\left\{\mathfrak{s},\mathfrak{s}\right\}}\text{mathbase}{c}^{\left\{\mathfrak{s}\right\}}.$$

$$\frac{1}{8}={\left({\mathbf{b}}^{\left\{\mathfrak{f}\right\}}\times {\mathbf{c}}^{\left\{\mathfrak{f}\right\}}\right)}^T{\mathbf{A}}^{\left\{\mathfrak{f},\mathfrak{s}\right\}}{\mathbf{c}}^{\left\{\mathfrak{s}\right\}}=\sum _{k\ge 0}{\omega }_k{{\gamma }^k}^T\text{mathbase}{c}^{\left\{\mathfrak{s}\right\}}.$$

$$\frac{1}{8}={\left({\mathbf{b}}^{\left\{\mathfrak{s}\right\}}\times {\mathbf{c}}^{\left\{\mathfrak{s}\right\}}\right)}^T{\mathbf{A}}^{\left\{\mathfrak{s},\mathfrak{f}\right\}}{\mathbf{c}}^{\left\{\mathfrak{f}\right\}}={\left(\text{mathbase}{b}^{\left\{\mathfrak{s}\right\}}\times \text{mathbase}{c}^{\left\{\mathfrak{s}\right\}}\right)}^T\text{mathbase}{A}^{\left\{\mathfrak{s},\mathfrak{s}\right\}}\text{mathbase}{c}^{\left\{\mathfrak{s}\right\}}.$$

Full text

{AMS}

65L05, 65L06

1 Introduction

In this paper, we consider the additively partitioned ordinary differential equation (ODE)

[TABLE]

where component $f^{\left\{\mathfrak{f}\right\}}$ represents the fast dynamics of the system, and component $f^{\left\{\mathfrak{s}\right\}}$ the slow dynamics. This structure models a feature appearing in many dynamical systems of practical interest: multiple characteristic time scales.

Multirate time integration methods are designed to efficiently solve Eq. 1 by using different timesteps for the fast and slow components. First explored by Rice [28] and Andrus [2, 3], the multirating strategy has been expanded to numerous types of traditional time integration methods. This includes Runge–Kutta methods [7, 13, 14, 23, 24, 33, 29], linear multistep methods [11, 18, 31], Rosenbrock-W methods [12], extrapolation methods [8, 10], Galerkin discretizations [25], and combined multiscale methodologies [9].

Multirate infinitesimal step (MIS) methods, first proposed by Knoth and Wolke [22], and later extended by others [21, 36, 37, 38, 40], introduce a new multirating philosophy in which the fast method solves a modified ODE that advances the solution between slow stages. While the slow system is solved discretely, the fast system can be solved with arbitrary small steps, hence the naming “infinitesimal step”. In [15], Günther and Sandu cast MIS methods into the General-structure Additive Runge–Kutta (GARK) framework. This framework was subsequently leveraged by Sandu in [30] to create the multirate infinitesimal GARK (MRI-GARK) class of methods. One step of an MRI-GARK method advances the solution from $t_{n}$ to $t_{n}+H$ by

[TABLE]

where $\mathbase{c}^{\left\{\mathfrak{s}\right\}}$ are the slow method abscissae, and the modified fast ODEs $v_{i}^{\prime}=\dots$ advance the solution between the slow stages. In [30], Sandu presents MRI-GARK methods Eq. 2 of orders up to four that are explicit or implicit in the fast and slow systems but are not coupled across partitions. This work also provides new techniques in investigating the stability of partitioned methods that we have adopted in our paper. Recent developments in the field include the work of Sexton and Reynolds [39] where a new structure for fast integration weights is considered and shown to help reduce order conditions; the “relaxed MIS” methods derived retain the same order as traditional MIS and it is possible to pair them for error control and adaptivity purposes.

This work extends the MRI-GARK family [30] to include implicit methods with coupling between slow and fast systems. We construct two new families of schemes designed to offer improved stability compared to decoupled MRI-GARK methods. The first family is step predictor-corrector MRI-GARK methods that perform a coupled discrete prediction over the entire timestep, then use that information to perform an infinitesimal step correction with small timesteps. This is similar to the strategy used by Rice in [28] and Savcenco in [34, 35]. The second family is internal stage predictor-corrector MRI-GARK methods that alternate between coupled discrete predictor stages and infinitesimal step correction ones. Figure 1 shows the differences between the new strategies and previous approaches. Order condition theories and stability analyses are developed for both new families of methods, and schemes up to order four are designed. Numerical experiments are employed to verify the theoretical findings.

The paper is organized as follows. The new family of step predictor corrector MRI-GARK schemes is introduced in Section 2, followed by its order condition theory and the stability analysis. In Section 3, internal stage predictor corrector MRI-GARK schemes are defined and their order conditions and stability are established. Numerical results are reported in Section 4, and concluding remarks are drawn in Section 5. Appendix A presents the lists of coefficients and stability plots for the newly developed methods.

2 Step predictor-corrector MRI-GARK methods

One coupling strategy commonly used in discrete multirate methods is a predictor-corrector approach, where the predictor evolves the entire system, while the corrector is only applied to the fast partition whose solution was “predicted” inaccurately (see [28, 34, 35]). First, a combined Runge–Kutta macro-step is taken which serves as the predictor. The fast parts of the predicted stages are inaccurate and are refined by sub-stepping the fast component only. Approximations of the slow values needed during the micro-steps are obtained from interpolating the slow predicted values. The step predictor-corrector MRI-GARK methods, as depicted in Fig. 1(b), can be viewed as an extreme case of this coupling strategy where the multirate ratio is infinite, i.e., the corrector takes infinitely many steps to refine the fast solution.

2.1 Method definition

We start with a “slow” Runge–Kutta base method

[TABLE]

with $s^{\left\{\mathfrak{s}\right\}}$ stages. Unlike other multirate infinitesimal strategies, the base method is not restricted to be explicit or diagonally implicit.

Definition 2.1** (Step predictor corrector MRI-GARK methods).**

One step of a step predictor-corrector MRI-GARK (SPC-MRI-GARK) scheme applied to Eq. 1 is given by

[TABLE]

*where f_{j}\coloneqq f\big{(}t_{n}+\mathbase{c}^{\left\{\mathfrak{s}\right\}}_{j}\,H,Y_{j}\big{)} and f_{j}^{\left\{\mathfrak{s}\right\}}\coloneqq f^{\left\{\mathfrak{s}\right\}}\big{(}t_{n}+\mathbase{c}^{\left\{\mathfrak{s}\right\}}_{j}\,H,Y_{j}\big{)}. *

Definition 2.2** (Slow tendency coefficients [30, Definition 2.2]).**

The time-dependent coefficients in Eq. 4b are defined as polynomials:

[TABLE]

Remark 2.3** (Embedded method).**

An embedded solution for an SPC-MRI-GARK method can be computed by solving the additional ODE

[TABLE]

*which uses the embedded polynomials $\widehat{\gamma}_{i}$ and produces a solution of a different order. We note that although this additional integration can be expensive, it can be done in parallel with Eq. 4b. *

Consider the trivial partitioning $f^{\left\{\mathfrak{f}\right\}}=0$, $f^{\left\{\mathfrak{s}\right\}}=f$ of Eq. 1. In this case, it is natural to expect an SPC-MRI-GARK method to degenerate into the slow base method. Note that the final solution of Eq. 4b simplifies to

[TABLE]

Thus, we enforce the condition

[TABLE]

An important special case of (1) are component partitioned systems:

[TABLE]

One step of an SPC-MRI-GARK method Eq. 4 applied to Eq. 8 reads:

[TABLE]

where f^{\left\{\mathfrak{f}\right\}}_{j}\coloneqq f^{\left\{\mathfrak{f}\right\}}\big{(}t_{n}+\mathbase{c}^{\left\{\mathfrak{s}\right\}}_{j}\,H,Y^{\left\{\mathfrak{f}\right\}}_{j},Y^{\left\{\mathfrak{s}\right\}}_{j}\big{)} and f^{\left\{\mathfrak{s}\right\}}_{j}\coloneqq f^{\left\{\mathfrak{s}\right\}}\big{(}t_{n}+\mathbase{c}^{\left\{\mathfrak{s}\right\}}_{j}\,H,Y^{\left\{\mathfrak{f}\right\}}_{j},Y^{\left\{\mathfrak{s}\right\}}_{j}\big{)}. With Eqs. 6 and 7 the internal ODE integrates (and corrects) only the fast component, while the slow component is solved with the traditional base Runge–Kutta method (3).

2.2 Order conditions

Following [30], we use an arbitrarily accurate Runge–Kutta method $\left(A^{\left\{\mathfrak{f},\mathfrak{f}\right\}},b^{\left\{\mathfrak{f}\right\}},c^{\left\{\mathfrak{f}\right\}}\right)$ to discretize the continuous ODE appearing in the method formulation, which casts the SPC-MRI-GARK scheme (4) into the GARK framework. The discrete corrector stages, denoted $Y^{\left\{\mathfrak{f},\right\}}_{i}$, are computed as

[TABLE]

where f^{\left\{\mathfrak{f},\right\}}_{j}\coloneqq f^{\left\{\mathfrak{f}\right\}}\big{(}t_{n}+c^{\left\{\mathfrak{f}\right\}}_{j}\,H,Y^{\left\{\mathfrak{f},\right\}}_{j}\big{)} and the superscript $\times k$ denotes the elementwise vector power. Similarly, the final solution reads

[TABLE]

Now, the corresponding GARK tableau for an SPC-MRI-GARK method is

[TABLE]

with $c^{\left\{\mathfrak{f},\mathfrak{s}\right\}}=\sum_{k\geq 0}A^{\left\{\mathfrak{f},\mathfrak{f}\right\}}\,c^{\left\{\mathfrak{f}\right\}\times k}\,{\gamma^{k}}^{T}\,\mathbbm{1}^{\left\{\mathfrak{s}\right\}}$.

2.2.1 Internal consistency

Theorem 2.4** (Internal consistency conditions).**

An SPC-MRI-GARK method Eq. 4 satisfies the “internal consistency” conditions

[TABLE]

for any fast method iff the following conditions hold:

[TABLE]

Proof 2.5**.**

All internal consistency equations are automatically satisfied except for the following one, which needs to be imposed explicitly:

[TABLE]

*It is easy to confirm Eq. 10 is sufficient to satisfy this condition, and thus, internal consistency. Since the equality must hold for all $A^{\left\{\mathfrak{f},\mathfrak{f}\right\}}$, it must hold when all $A^{\left\{\mathfrak{f},\mathfrak{f}\right\}}\,c^{\left\{\mathfrak{f}\right\}\times k}$ are linearly independent. Matching powers of the left- and right-hand sides proves the necessity of Eq. 10. *

If an SPC-MRI-GARK method has a slow base method Eq. 3 of order two, then internal consistency is sufficient to guarantee the method is order two [32].

2.2.2 Fourth order conditions

In this section, we derive order conditions of the SPC-MRI-GARK schemes for up to order four. First, we define a set of useful coefficients.

Definition 2.6** (Some useful coefficients [30, Definition 3.3]).**

Consider the “bushy” Butcher tree [16]

[TABLE]

where $\tau\in T$ is the tree of order one and $[\cdot]$ is the operation of joining subtrees by a root. An arbitrarily accurate Runge–Kutta method $\left(A^{\left\{\mathfrak{f},\mathfrak{f}\right\}},b^{\left\{\mathfrak{f}\right\}},c^{\left\{\mathfrak{f}\right\}}\right)$ satisfies the following equations:

[TABLE]

Theorem 2.7** (Fourth order coupling conditions).**

An internally consistent SPC-MRI-GARK method Eq. 4 satisfying Eq. 7 has order four iff the slow base scheme Eq. 3 has order at least four, and the following coupling conditions hold:

[TABLE]

Proof 2.8**.**

An internally consistent GARK scheme is order four iff the base methods are order four and the 12 coupling conditions up to order four are satisfied [32]. We proceed with checking each coupling condition.

Condition 3a

The first third order condition gives Eq. 13a:

[TABLE]

Condition 3b

The other third order condition is automatically satisfied if the slow base method is order three:

[TABLE]

Condition 4a

The first fourth order condition gives Eq. 13b:

[TABLE]

Condition 4b

This condition is automatically satisfied if the slow base method is order four:

[TABLE]

Condition 4c

This condition proves Eq. 13c:

[TABLE]

Condition 4d

This condition is automatically satisfied if the slow base method is order four:

[TABLE]

Condition 4e

This condition is the redundant since it is the difference of Condition 3a and 4a:

[TABLE]

Condition 4f

This condition proves Eq. 13d:

[TABLE]

Condition 4g

The following condition is identical to Condition 4f:

[TABLE]

Condition 4h

This condition is automatically satisfied if the slow base method is order four:

[TABLE]

Condition 4i

This condition is automatically satisfied if the slow base method is order four:

[TABLE]

Condition 4j

This condition is automatically satisfied if the slow base method is order four:

[TABLE]

Remark 2.9**.**

*In the proof of Theorem 2.7, all coupling order conditions that start with $\mathbf{b}^{\left\{\mathfrak{s}\right\}T}$ collapse onto the order conditions of the slow base method. In the context of two-trees [4, 32], trees containing both slow and fast nodes with a slow root can be recolored into purely slow trees. The purely fast trees are of no concern since the fast base method is arbitrarily accurate. The remaining trees contain slow and fast nodes with a fast root, which correspond to coupling conditions Eq. 13 that must be explicitly enforced through the choice of $\gamma$. *

2.3 Stability analysis

2.3.1 Scalar stability analysis

Consider the partitioned, linear, scalar test problem

[TABLE]

and let $z^{\left\{\mathfrak{f}\right\}}\coloneqq H\,\lambda^{\left\{\mathfrak{f}\right\}}$, $z^{\left\{\mathfrak{s}\right\}}\coloneqq H\,\lambda^{\left\{\mathfrak{s}\right\}}$, and $z\coloneqq z^{\left\{\mathfrak{f}\right\}}+z^{\left\{\mathfrak{s}\right\}}$. Applying the SPC-MRI-GARK method (4) to Eq. 14 yields

[TABLE]

which leads to the stability function

[TABLE]

where, following [30]:

[TABLE]

Of special interest are cases when a partition becomes infinitely stiff. If the base method has bounded internal stability, the stability function Eq. 15 enjoys the following property:

[TABLE]

Although Eq. 16b cannot be zero for all $z^{\left\{\mathfrak{f}\right\}}$ due to the linear independence of $\varphi$ functions, its modulus is bounded for $z^{\left\{\mathfrak{f}\right\}}\in\mathbb{C}^{-}$.

2.3.2 Matrix stability analysis

Following [23, 30], consider the matrix test problem

[TABLE]

The change of variables that produces $\boldsymbol{\Omega}$ [30]

[TABLE]

allows the matrix eigenvalues to be written as linear combinations of the diagonal entries: $\xi\,\lambda^{\left\{\mathfrak{f}\right\}}+(1-\xi)\,\lambda^{\left\{\mathfrak{s}\right\}}$ and $(1-\xi)\,\lambda^{\left\{\mathfrak{f}\right\}}+\xi\,\lambda^{\left\{\mathfrak{s}\right\}}$. The coupling between the fast and slow variables is controlled by $\xi$. Values close to zero indicate the slow system is weakly influenced by the fast one, while values close to one indicate the fast system is weakly influenced by the slow one.

Let

[TABLE]

The component partitioned SPC-MRI-GARK method Eq. 9 applied to the matrix test problem Eq. 17 gives

[TABLE]

We write the solution of the ODE as

[TABLE]

The transfer matrix for the matrix test problem can be written as

[TABLE]

where $\mathfrak{Y}(Z)$ is the internal stability matrix:

[TABLE]

2.4 Construction of practical SPC-MRI-GARK methods

We develop new implicit SPC-MRI-GARK methods of up to order four. Their coefficients are presented in Section A.1. The base methods are chosen to be existing, high-quality schemes that have either singly diagonally implicit (SDIRK) or explicit first stage single diagonally implicit (ESDIRK) structures. These offer a nice balance between stability and computational complexity. We note that explicit and fully implicit base methods can be employed as well. The $\gamma(t)$ coupling coefficients for each method are determined by first enforcing the order conditions, and then using remaining free parameters to optimize for stability. Plots of the scalar and matrix stability regions are provided in Figs. 5 and 7, respectively. These regions are significantly larger than those of the decoupled MRI-GARK counterparts developed in [30].

3 Internal stage predictor-corrector MRI-GARK methods

Traditional multirate infinitesimal methods subdivide the integration interval $[t_{n},t_{n+1}]$ into subintervals $[t_{n}+\mathbase{c}^{\left\{\mathfrak{s}\right\}}_{i}\,H,t_{n}+\mathbase{c}^{\left\{\mathfrak{s}\right\}}_{i+1}\,H]$, and solve a fast ODE over each subinterval. This advances the solution from one abscissa to the next, and then to the final solution. As illustrated in Fig. 1(c), an internal stage predictor-corrector MRI-GARK method follows this strategy, but also incorporates a predictor-corrector strategy similar to that used in SPC-MRI-GARK schemes. On each subinterval, the solution is first predicted with a traditional Runge-Kutta stage calculation. Next, the fast components are refined by solving an ODE which uses previous predictor and corrector stages, as well as the current predictor stage, to implement the slow tendencies.

3.1 Method definition

Again, we start with a slow Runge–Kutta base method Eq. 3, but now enforce that is has a diagonally implicit structure and the abscissae are non-decreasing:

[TABLE]

This ensures that each ODE between stages is not integrated backward in time. We define the abscissa increments:

[TABLE]

The final integration from $\mathbase{c}^{\left\{\mathfrak{s}\right\}}_{s^{\left\{\mathfrak{s}\right\}}}$ to $1$ can introduce special cases that increase the complexity of the notation, order conditions, and stability analysis. We will impose that the base slow method is stiffly accurate [17], which makes the last stage equal to the final solution, and simplifies the subsequent analyses. This comes at no loss of generality since we can always rewrite a Runge–Kutta method into a reducible, but stiffly accurate form. In Butcher tableau notation we have:

[TABLE]

Definition 3.1** (Internal stage predictor corrector MRI-GARK methods).**

One step of an internal stage predictor-corrector MRI-GARK (IPC-MRI-GARK) scheme applied to Eq. 1 is given by

[TABLE]

*with $f^{\left\{\mathfrak{s}\right\}\ast}_{j}\coloneqq f^{\left\{\mathfrak{s}\right\}}\left(T_{j},Y^{\ast}_{j}\right)$ and $f^{\ast}_{j}\coloneqq f\left(T_{j},Y^{\ast}_{j}\right)$. Stages and functions with an asterisk are predictor values, and terms without the asterisk are corrector values. In order to enforce that only previously computed stages appear in the ODE, we require that $\gamma_{i,j}(\tau)=0$ for $j\geq i$ and $\psi_{i,j}(\tau)=0$ for $j>i$. *

Once again, we can take each $\gamma_{i,j}(t)$ and $\psi_{i,j}(t)$ to be polynomial in time. These and their integral terms $\widetilde{\gamma}_{i,j}(t)$, $\overline{\gamma}_{i,j}$, $\widetilde{\psi}_{i,j}(t)$, $\overline{\psi}_{i,j}$ are defined analogously to Eq. 5. The capitalized versions are used to denote the matrices of coefficients.

Remark 3.2** (Embedded method).**

Following the strategy used in [30], an embedded solution can be obtained via the additional integration

[TABLE]

With the trivial partitioning $f^{\left\{\mathfrak{s}\right\}}=f$, $f^{\left\{\mathfrak{f}\right\}}=0$, the corrector stages simplify to

[TABLE]

The IPC-MRI-GARK method becomes a $2s^{\left\{\mathfrak{s}\right\}}$ stage Runge–Kutta method with $s^{\left\{\mathfrak{s}\right\}}$ predictor stages and $s^{\left\{\mathfrak{s}\right\}}$ corrector stages. In the absence of the fast component, it is natural to expect the predictor and corrector stages to coincide and for the method to degenerate into the slow base scheme. Matching coefficients of Eq. 20 to those of the predictor stage of Eq. 19 gives the self-consistency conditions

[TABLE]

where $\mathbase{T}^{\left\{\mathfrak{s}\right\}}$ is the strictly lower triangular part of $\mathbase{A}^{\left\{\mathfrak{s},\mathfrak{s}\right\}}$ and

[TABLE]

Remark 3.3** (Repeated abscissae).**

When $\mathbase{c}^{\left\{\mathfrak{s}\right\}}_{i}=\mathbase{c}^{\left\{\mathfrak{s}\right\}}_{i-1}$, the fast function disappears from the ODE in Eq. 19b as it is scaled by zero. The corrector stage simplifies to

[TABLE]

*Clearly, an ODE solver is no longer needed to compute $Y_{i}$. This can be viewed as modifying (the slow part of) the initial conditions for the next step’s ODE. *

For component partitioned systems Eq. 8, an IPC-MRI-GARK step reads:

[TABLE]

where $\widetilde{\delta}_{i,j}\left(\tfrac{\theta}{H}\right)=\widetilde{\gamma}_{i,j}\left(\tfrac{\theta}{H}\right)+\widetilde{\psi}_{i,j}\left(\tfrac{\theta}{H}\right)$ and f^{\left\{\mathfrak{f}\right\}\ast}_{j}\coloneqq f^{\left\{\mathfrak{f}\right\}}\big{(}T_{j},Y^{\left\{\mathfrak{f}\right\}\ast}_{j},Y^{\left\{\mathfrak{s}\right\}\ast}_{j}\big{)}.

3.2 Order conditions

Following Section 2.2, we look to utilize GARK order condition theory to derive order conditions for IPC-MRI-GARK methods. Again, we apply an arbitrarily accurate Runge–Kutta method $\left(A^{\left\{\mathfrak{f},\mathfrak{f}\right\}},b^{\left\{\mathfrak{f}\right\}},c^{\left\{\mathfrak{f}\right\}}\right)$ to discretize the ODEs and recover the GARK stages and GARK tableau. We use the labels ${\rm Phys.~{}Rev.~{}E}$ and to denote predictor and corrector stages, respectively. Also we define $Y^{\left\{\mathfrak{f},i\right\}}_{k}$ to be the $k$-th stage of the discretized ODE between abscissae $\mathbase{c}^{\left\{\mathfrak{s}\right\}}_{i-1}$ and $\mathbase{c}^{\left\{\mathfrak{s}\right\}}_{i}$. Now, the $i$-th step of Eq. 19 is composed of the GARK stages

[TABLE]

with f^{\left\{\mathfrak{f},i\right\}}_{j}\coloneqq f^{\left\{\mathfrak{f}\right\}}\big{(}T_{i-1}+{\Delta\mathbase{c}}^{\left\{\mathfrak{s}\right\}}_{i}\,c^{\left\{\mathfrak{f}\right\}}_{j},Y^{\left\{\mathfrak{f},i\right\}}_{j}\big{)} and f^{\left\{\sigma,\nu\right\}}_{j}\coloneqq f^{\left\{\sigma\right\}}\big{(}T_{j},Y^{\left\{\sigma,\nu\right\}}_{j}\big{)} for $\sigma\in\{\mathfrak{f},\mathfrak{s}\}$ and $\nu\in\{{\rm Phys.~{}Rev.~{}E},$. Now, we simplify $Y^{\left\{\mathfrak{f},\right\}}_{i}$ to obtain:

[TABLE]

The stages of the discretized ODE simplify to

[TABLE]

The coefficients appearing in the stages can be organized into the following GARK tableau:

[TABLE]

The unspecified entries are

[TABLE]

with

[TABLE]

and $L\in\mathbbm{R}^{s^{\left\{\mathfrak{s}\right\}}\times s^{\left\{\mathfrak{s}\right\}}}$ is a lower shift matrix with entries $L_{i,j}=\delta_{i,j+1}$.

3.2.1 Internal Consistency

Theorem 3.4** (Internal consistency conditions).**

An IPC-MRI-GARK method Eq. 19 fulfills the “internal consistency” conditions

[TABLE]

for any fast method iff the following conditions hold:

[TABLE]

Proof 3.5**.**

All internal consistency equations are automatically satisfied except

[TABLE]

*It is easy to confirm Eq. 24 is sufficient to satisfy this condition, and thus, internal consistency. Since the equality must hold for all $A^{\left\{\mathfrak{f},\mathfrak{f}\right\}}$, it must hold when all $A^{\left\{\mathfrak{f},\mathfrak{f}\right\}}c^{\left\{\mathfrak{f}\right\}\times k}$ are linearly independent. Matching powers of the left- and right-hand sides proves the necessity of Eq. 24. *

Like with SPC-MRI-GARK methods, internal consistency and a slow base method of order two guarantees an IPC-MRI-GARK method is order two [32].

3.2.2 Fourth order conditions

In this section, we derive order conditions of the IPC-MRI-GARK schemes for up to order four.

Lemma 3.6** (Intermediate matrix products).**

The coefficients of the GARK tableau Eq. 23 satisfy

[TABLE]

Theorem 3.7** (Fourth order coupling conditions).**

An internally consistent IPC-MRI-GARK method Eq. 19 satisfying Eq. 21 has order four iff the slow base scheme has order at least four, and the following coupling conditions hold:

[TABLE]

Proof 3.8**.**

An internally consistent GARK scheme is order four iff the base methods are order four and the 12 coupling conditions up to order four are satisfied [32]. We proceed with checking each coupling condition.

Condition 3a

By using Eq. 25c, the first third order condition gives Eq. 26a:

[TABLE]

Condition 3b

The other third order condition is expanded with Eq. 25b to get Eq. 26b:

[TABLE]

Condition 4a

By using Eq. 25c, the first fourth order condition gives Eq. 26c:

[TABLE]

Condition 4b

We derive Eq. 26d with Eq. 25b:

[TABLE]

Condition 4c

We derive Eq. 26e with Eq. 25c:

[TABLE]

Condition 4d

We derive Eq. 26f with Eq. 25b:

[TABLE]

Condition 4e

We derive Eq. 26g with Eq. 25c:

[TABLE]

Condition 4f

We derive Eq. 26h with Eq. 25b:

[TABLE]

Condition 4g

We derive Eq. 26i with Eq. 25a:

[TABLE]

Condition 4h

We derive Eq. 26j with Eq. 25b:

[TABLE]

Condition 4i

We derive Eq. 26k with Eq. 25c:

[TABLE]

Condition 4j

This condition is equivalent to condition 4d since the fast base method has an arbitrarily large stage order, and thus, $\mathbf{A}^{\left\{\mathfrak{f},\mathfrak{f}\right\}}\mathbf{c}^{\left\{\mathfrak{f}\right\}}=\frac{1}{2}\mathbf{c}^{\left\{\mathfrak{f}\right\}\times 2}$:

[TABLE]

3.3 Linear stability analysis

3.3.1 Scalar stability analysis

We revisit the scalar linear test problem Eq. 14 now for IPC-MRI-GARK methods. The predictor stages in vector form are given by

[TABLE]

The internal ODEs become

[TABLE]

Integrating and substituting in the predicted stages Eq. 27 gives

[TABLE]

Applying Eq. 19 to Eq. 14 gives the scalar stability function

[TABLE]

with

[TABLE]

It can be verified that Eq. 16a still holds. As the slow part of the test problem becomes infinitely stiff, however, we would like the stability function to be bounded. One natural way to enforce this is by ensuring $\mathfrak{M}$, and thus $\mathfrak{M}^{-1}$, remains bounded in the limit. The last two terms in $\mathfrak{M}$ are problematic since they are $\order{z^{\left\{\mathfrak{s}\right\}}missing}$. If $\mathbase{A}^{\left\{\mathfrak{s},\mathfrak{s}\right\}}$ is invertible, the following condition ensures these terms cancel in the limit:

[TABLE]

Note that $\mu$ and $\nu$ are sums over linearly independent $\varphi$ functions. By matching terms in this summation, we arrive at the stability simplifying assumption

[TABLE]

If $\Psi$ and $\Gamma$ are degree zero polynomials, then Eq. 21 automatically ensures Eq. 29 is satisfied.

3.3.2 Matrix stability analysis

Now we consider the component partitioned IPC-MRI-GARK method Eq. 22 applied to the matrix test problem Eq. 17. First, we define the following intermediate quantities:

[TABLE]

The predictor stages become

[TABLE]

The fast internal ODEs become

[TABLE]

The solution to the system of ODEs gives the corrector stages

[TABLE]

The combined fast and slow corrector stages are

[TABLE]

The stability matrix is given by

[TABLE]

with

[TABLE]

3.4 Construction of practical methods

We develop new implicit IPC-MRI-GARK methods up to order four which are presented in Section A.2. The second order base methods are reused from SPC-MRI-GARK, but the third and fourth order base methods are custom due to the nondecreasing abscissae constraint. Upon deriving a parameterized family of $\Gamma$ and $\Psi$ coefficients that satisfy the coupling order conditions, we use free coefficients to satisfy the stability simplifying assumption Eq. 29. Any remaining parameters are used to optimize the size of the stability region. Plots of the scalar and matrix stability regions are provided in Figs. 8 and 9, respectively. Compared to SPC-MRI-GARK methods, we found it significantly more challenging to achieve large stability regions at high orders.

4 Numerical results

In this section, we present the numerical tests performed on the SPC-MRI-GARK and IPC-MRI-GARK methods.

4.1 Additive partitioning: the Gray–Scott model

The first test problem considered is the Gray–Scott reaction-diffusion PDE [26]:

[TABLE]

It is solved over the 2D spatial domain $[0,1]\times[0,1]$, which is discretized with second order finite differences. The timespan is taken to be $[0,30]$, and the model parameters are $\varepsilon_{u}=0.0625$, $\varepsilon_{v}=0.0312$, $\mathfrak{k}=0.0520$, and $\mathfrak{f}=0.0180$. The linear diffusions terms of Eq. 31 make up the slow partition while the nonlinear reaction terms make up the fast partition.

MATLAB is used to carry out the convergence experiments. ODEs that appear within the integrators are solved using ode45 with the tolerances abstol = reltol = 1e-10. Convergence diagrams for the new methods presented in Appendix A are shown in Fig. 2. The numerical orders of accuracy are consistent with theoretical orders.

4.2 Component partitioning: the KPR problem

For a component partitioned test problem of the form Eq. 8, we use the KPR system [8] as a multi-scale extension to the scalar Prothero-Robinson [6, 17, 27] problem. We define the system as:

[TABLE]

The tests are performed from $t=0$ to $t=5\pi/2$ with the initial condition coming from evaluating Eq. 32b at $t=0$. From the exact solution we can also see that the differences in the fast and slow time scales are driven by $\omega$ and not $\lambda^{\left\{\mathfrak{f}\right\}}$ and $\lambda^{\left\{\mathfrak{s}\right\}}$.

The fast integration Eq. 22b is also carried out using ode45 solver with abstol = reltol = 1e-10. The convergence diagrams reported in Fig. 3 indicate that the methods perform at their theoretical orders for this problem.

4.3 Multirate performance: the inverter chain problem

We also consider the inverter chain model of [24] given by the equations

[TABLE]

with $U_{0}=0$, $U_{op}=5$, $U_{T}=1$, $\Gamma=100$, and

[TABLE]

The initial conditions of the system are

[TABLE]

and the input signal is taken to be

[TABLE]

For the numerical experiments, we use $m=500$ and a timespan of $[0,100]$. As the signal propagates through the circuit, only a small percentage of the inverters experience a change in voltage while the other inverters maintain a constant voltage. A componentwise partitioning of Eq. 33 is used where the fast components come from a sliding window that follows the signal, and the remaining components form the slow partition.

A C implementation of Eq. 33 is used to measure the performance gains provided by SPC-MRI-GARK and IPC-MRI-GARK over a single rate base method of the same order. For order two, we compare SPC SDIRK2(1)2 from Section A.1.1 and IPC SDIRK2(1)2 from Section A.2.1 to their shared base method SDIRK2(1)2. The results are plotted in Fig. 4(a). Figure 4(b) compares SPC SDIRK3(2)4 from Section A.1.3 to its base method SDIRK3(2)4 and to IPC SDIRK3(2)5 from Section A.2.3. Finally, Fig. 4(c) compares SPC ESDIRK4(3)6 from Section A.1.6 to its base method ESDIRK4(3)6 and to SPC SDIRK4(3)5 from Section A.1.5. We did not include results for IPC SDIRK4(3)6 as it was only stable for timesteps much smaller than those used for the SPC methods. This observation is consistent with the stability regions presented in Figs. 8 and 9.

Fixed timesteps were used for the coupled MRI-GARK methods, as well as for the method ESDIRK5(4)7[2]SA2 from [20] used to solve the internal ODEs. In the experiments, ten timesteps were taken to solve these internal ODEs, except for IPC SDIRK2(1)2 and SPC SDIRK3(2)4 where five and 15 steps were used, respectively.

In all of the performance results presented in Fig. 4, the multirate methods are able to achieve a desired accuracy in significantly less time than the single rate schemes. The best results occurred at order two where the speedup ranged from 8 to 60. This can be attributed to the excellent multirate characteristics of the inverter chain problem as well as the flexibility of mutirate infinitesimal methods to use any method to solve the modified fast ODEs.

5 Conclusions and future work

This work extends the class of multirate infinitesimal GARK schemes developed in [30] to include coupled methods. Such methods compute (some of) the stages by solving implicit systems that involve both the fast and the slow components, which gives their “coupled” character. The coupled approach allows us to construct multirate infinitesimal schemes with improved stability for stiff systems with multiple scales, at the additional cost of solving more complex, or larger, nonlinear systems.

Two approaches to formulating the coupling are studied herein. Both of them employ a predictor-corrector structure. The first approach, named step predictor-corrector MRI-GARK, starts with computing all predictor stages in a coupled fashion. The predicted stages are then used to formulate a modified fast ODE, and a single infinitesimal integration is carried out to correct the fast component of the system. The second approach, named internal stage predictor-corrector MRI-GARK, alternates prediction and correction stages. Specifically, each discrete predictor stage is followed by a corrector stage, which integrates a modified fast ODE system and corrects the fast components of that stage.

Elegant formulations of the order conditions for both families of methods are developed, and stability requirements for practical methods are analyzed. Methods of order up to four are constructed. Numerical tests verify the orders of convergence on additive and component partitioned cases. Finally, we demonstrate computational efficiency of these multirate methods when compared to their single rate counterparts. Our numerical experiments indicate a performance edge for both MRI-GARK strategies compared to single rate ones, with SPC-MRI-GARK methods having slightly better performance in general. Our analysis also shows larger stability regions for SPC-MRI-GARK methods compared to IPC-MRI-GARK.

The succession of discrete and infinitesimal integration stages in the IPC-MRI-GARK schemes requires non-decreasing abscissae for the slow base method. This requirement, in conjunction with stability, can become difficult to satisfy for high order methods. One solution to alleviate this restriction is to construct MRI-GARK methods that compute their own initial conditions for the infinitesimal stage integration. The authors plan to study these extensions in future works.

Appendix A New MRI-GARK methods

Here, we present the newly derived SPC-MRI-GARK and IPC-MRI-GARK methods. In some cases, the exact representation of the method coefficients is too long to fit on a page, so the first 16 digits are provided. Exact coefficients are available in the supplementary materials. The stability regions are computed according to [30]:

[TABLE]

A.1 SPC-MRI-GARK methods

We use the following tableau representation for SPC-MRI-GARK methods:

[TABLE]

A.1.1 SDIRK2(1)2

This method is based on the two stage, second order method in [1].

[TABLE]

A.1.2 ESDIRK2(1)3

This method is based on TR-BDF2 in [5].

[TABLE]

A.1.3 SDIRK3(2)4

This method is based on SDIRK3M in [19].

[TABLE]

A.1.4 ESDIRK3(2)4

This method is based on the optimal four stage, third order ESDIRK method described in [19].

[TABLE]

A.1.5 SDIRK4(3)5

This method is based on SDIRK4M in [19].

[TABLE]

A.1.6 ESDIRK4(3)6

This method is based on ESDIRK4(3)6L[2]SA in [19] and has the property that the first and second column of coefficients are identical.

[TABLE]

A.1.7 Stability plots

Figures 5 and 7 show the scalar and matrix stability regions, respectively.

A.2 IPC-MRI-GARK methods

We will using the following tableau representation for IPC-MRI-GARK methods:

[TABLE]

A.2.1 SDIRK2(1)2

This method is based on the two stage, second order method in [1].

[TABLE]

A.2.2 ESDIRK2(1)3

This method is based on TR-BDF2 in [5].

[TABLE]

A.2.3 SDIRK3(2)5

[TABLE]

A.2.4 SDIRK4(3)6

[TABLE]

A.2.5 Stability plots

Figures 8 and 9 show the scalar and matrix stability regions, respectively.