 # A Class of Multirate Infinitesimal GARK Methods

Differential equations arising in many practical applications are characterized by multiple time scales. Multirate time integration seeks to solve them efficiently by discretizing each scale with a different, appropriate time step, while ensuring the overall accuracy and stability of the numerical solution. In a seminal paper Knoth and Wolke (APNUM, 1998) proposed a hybrid solution approach: discretize the slow component with an explicit Runge-Kutta method, and advance the fast component via a modified fast differential equation. The idea led to the development of multirate infinitesimal step (MIS) methods by Wensch et al. (BIT, 2009.)Günther and Sandu (BIT, 2016) explained MIS schemes as a particular case of multirate General-structure Additive Runge-Kutta (MR-GARK) methods. The hybrid approach offers extreme flexibility in the choice of the numerical solution process for the fast component. This work constructs a family of multirate infinitesimal GARK schemes (MRI-GARK) that extends the hybrid dynamics approachin multiple ways. Order conditions theory and stability analyses are developed, and practical explicit and implicit methods of up to order four are constructed. Numerical results confirm the theoretical findings. We expect the new MRI-GARK family to be most useful for systems of equations with widely disparate time scales, where the fast process is dispersive, and where the influence of the fast component on the slow dynamics is weak.

Comments

There are no comments yet.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Many dynamical systems of practical interest consist of multiple components that evolve at different characteristic time scales. As a representative model we consider a two-way additively partitioned ordinary differential equation (ODE) driven by a slow process

and a fast process , acting simultaneously:

 y′=f(t,y)=f{s}(t,y)+f{f}(t,y),y(t0)=y0. (1)

Multirate methods are efficient numerical solution strategies for (1) that use small step sizes to discretize the fast components, and large step sizes to discretize the slow components. The approach goes back to the pioneering work on multirate Runge-Kutta methods of Rice  and Andrus [2, 3].

The main approach for the construction of multirate methods is to employ traditional integrators with different time steps for different components, and to couple the fast and slow components such as to ensure the overall stability and accuracy of the discretization. Following this philosophy multirate schemes have been developed in the context of Runge-Kutta methods [9, 19, 20, 30, 31], linear multistep methods [15, 26, 39], Rosenbrock-W methods , extrapolation methods [10, 11, 14, 40], Galerkin discretizations , and combined multiscale methodologies .

Multirate time integration has been adopted by many applications including the simulation of electronic circuits [20, 8], subsurface flows , ocean modeling [12, 46], energy conversion , multibody dynamics , electromagnetic fields , and fluids in complex geometries , to name just a few.

In a seminal paper Knoth and Wolke  proposed a hybrid approach to solve (1). First, the slow component is discretized with an -stage explicit Runge-Kutta method with increasing abscissae, . Next, the solution is advanced between the consecutive stages of this method by solving a modified fast ODE system. For autonomous systems (1) the solution process reads:

 Y{s}1\coloneqqyn,⎧⎨⎩v′i=f{f}(vi)+∑i−1j=1ai,j−ai−1,jci−ci−1f{s}(Y{s}j),τ∈[0,(ci−ci−1)H],with vi(0)=Y{s}i−1,Y{s}i\coloneqqvi((ci−ci−1)H),i=2,…,s;v′s+1=f{f}(v)+∑sj=1bj−as,j1−csf{s}(Y{s}j),τ∈[0,(1−cs)H],with vs+1(0)=Y{s}s;yn+1\coloneqqvs+1((1−cs)H). (2)

The modified fast ODE (2) is composed of the original fast component from (1) plus a piece-wise constant “slow tendency” term. Wensch, Knoth, and Galant  generalized (2) by adding linear combinations of stages (, ) to both the initial condition and to the slow tendency term. Since the fast ODE is solved exactly (or, equivalently, is solved numerically with an infinitesimally small step size) they aptly named the resulting schemes “multirate infinitesimal step” methods. The methods were interpreted as exponential integrators, and order conditions up to order three were derived. Schlegel, Knoth, Wolke, and Arnold casted MIS schemes as traditional partitioned Runge-Kutta methods , and furthered the theoretical and practical understanding of this family [28, 42, 34, 43].

The hybrid dynamical approach has proved fruitful not only in the construction of new numerical schemes, but also in the analysis of existing ones. For example, a hybrid approach has been used to abstract away the details of subsystem integration and study only the coupling aspects of waveform relaxation methods [5, 6]. The hybrid dynamical approach is also related to the engineering concept of co-simulation [17, 45].

Günther and Sandu  explained MIS schemes in the framework of General-structure Additive Runge-Kutta (GARK) methods . A multirate GARK (MR-GARK) method  integrates the slow component with a Runge–Kutta method and a large step size , and the fast component with another Runge–Kutta method and a small step size . Here represents the (integer) number of fast steps that are executed for each of the slow steps. One step of the method computes slow stages, denoted by , and fast stages, denoted by : equationparentequation

 f{f,λ}j \coloneqqf{f}(Y{f,λ}j),f{s}j\coloneqqf{s}(Y{s}j); Y{s}i=yn+Hs{s}∑j=1a{s,s}i,jf{s}j+hM∑λ=1s{f}∑j=1a{s,f,λ}i,jf{f,λ}j,1≤i≤s{s}, (3a) Y{f,λ}i=˜yn+λ−1M+Hs{s}∑j=1a{f,s,λ}i,jf{s}j+hs{f}∑j=1a{f,f}i,jf{f,λ}j,1≤i≤s{f}, (3b) ˜yn+λM =˜yn+λ−1M+hs{f}∑i=1b{f}if{f,λ}i,λ=1,…,M, (3c) yn+1 =˜yn+MM+Hs{s}∑i=1b{s}if{s}i. (3d)

The Butcher tableau for the MR-GARK method (3) is :

 (4)

This work constructs a family of multirate infinitesimal GARK schemes (named MRI-GARK) that extend the hybrid dynamics approach of [29, 48] in several ways. Time dependent “slow tendency” terms are used to construct the modified fast system. The MRI-GARK general order conditions theory is developed by leveraging the GARK accuracy theory . This allows the construction of the first fourth order multirate infinitesimal methods, as well as the first multirate infinitesimal schemes that are implicit in the slow component. Matrix stability analyses are carried out using a new simplified two-dimensional linear test problem.

The paper is organized as follows. The new family of MRI-GARK schemes is defined in Section 2. The order condition theory is developed in Section 3, and the stability analysis in Section 4. Practical explicit methods of order up to four are presented in Section 6, and decoupled implicit methods in Section 7. Numerical results are reported in Section 8, and conclusions are drawn in Section 9.

## 2 Multirate Infinitesimal GARK Methods

We start with a “slow” Runge-Kutta base method with non-increasing abscissae: equationparentequation

 (5a) and define the increments between consecutive stages of the base method: (5b)
###### Definition 1 (MRI-GARK methods for additively partitioned systems)

A Multirate Infinitesimal GARK (MRI-GARK) scheme applied to the additively partitioned system (1) advances the solution from to as follows: equationparentequation

 Y{s}1=yn, (6a) (6b) yn+1=Y{s}s{s}+1. (6c)

Linear combinations of the slow function values are added as forcing to the modified fast ODE system (6b); in order to use only already computed slow stages one needs for .

###### Definition 2 (Slow tendency coefficients)

It is convenient to define the time-dependent combination coefficients as polynomials in time, and to consider their integrals:

 γi,j(τ)\coloneqq∑k≥0γki,jτk,˜γi,j(x)\coloneqq∫x0γi,j(τ)dτ=∑k≥0γki,jxk+1k+1,¯¯¯γi,j\coloneqq˜γi,j(1). (7)

###### Remark 1 (Equal abscissae)

When two consecutive abscissae are equal, , the computation (6b) reduces to an explicit slow Runge-Kutta stage:

 Y{s}i+1=Y{s}i+Hi+1∑j=1(∫1τ=0γi,j(τ))f{s}(Y{s}j)=Y{s}i+Hi+1∑j=1¯¯¯γi,jf{s}(Y{s}j),

whenever for . We note that it is possible to choose , in which case (6b) is equivalent to a singly diagonally implicit Runge-Kutta stage.

###### Remark 2 (Embedded method)

The base slow scheme (5a) computes the main solution using the weights , and the embedded solution using the weights . The MRI-GARK scheme computes the main solution (6c) via the last stage (6b), that advances to by solving a modified fast system with the slow weights . An embedded solution is computed via an additional stage that advances to using with the modified weights .

###### Definition 3 (MRI-GARK methods for component partitioned systems)

Consider the component partitioned system:

 (8)

An MRI-GARK scheme applied to (8) computes the solution as follows: equationparentequation

 (9a) (9b) (9c)

For we require that such that only previously computed slow stage values are used in the formulation of the modified fast ODE. The additive (6) and component-based (9) formulations are equivalent to each other.

###### Example 1 (Second order MRI-GARK methods)

Consider the explicit midpoint method and the implicit trapezoidal method, each paired with a first order embedded scheme:

The explicit midpoint method is the slow component (5a) of the following second order explicit MRI-GARK scheme (6):

 (10)

The implicit trapezoidal method is the slow component (5a) of the following second order implicit MRI-GARK scheme (6):

 (11)

## 3 Order conditions

For accuracy analysis we replace the continuous fast process by a discrete Runge Kutta method of arbitrary accuracy and having an arbitrary number of stages . This approach casts the MRI-GARK scheme (6) into the multirate GARK framework (3), where each sub-step is carried out between one slow stage and the next. We denote the fast sub-steps by , where each advances the fast system from to .

###### Remark 3 (Matrix notation)

It is convenient to gather the gamma coefficients in the matrices , and express (7) as follows:

 Γ(τ)=∑k≥0Γkτk,˜Γ(x)=∑k≥0xk+1k+1Γk,¯¯¯¯Γ=∑k≥01k+1Γk.

The structure of the coefficient matrices is lower triangular, in the sense that:

 (12)

###### Lemma 1 (The MRI-GARK method (6) as a particular instance of a GARK method)

Consider the MRI-GARK scheme (6) where the continuous fast process is replaced by a discrete Runge Kutta method of arbitrary accuracy and having an arbitrary number of stages . The resulting computational process is a multirate GARK method (3) with the following Butcher tableau (4) components.

##### Fast component

equationparentequation

 = (13a) b{f} = (13b) c{f,f} = (13c)

Here is the total number of stages of the method and is the Kronecker product.

##### Slow component

equationparentequation

 A{s,s} ≡ (14a) b{s}T ≡ \mathbaseb{s}T=11{s}T¯¯¯¯Γ, (14b) c{s,s} = A{s,s}11{s}≡\mathbasec{s}=E¯¯¯¯Γ11{s}. (14c)
##### Slow-fast coupling

equationparentequation

 = (15a) A{s,f,i} = {Δ\mathbasec{s}igi+1b{f}T,i=1,…,s{s}−1,0s{s}×s{f},i=s{s}, c{s,f} = A{s,f}11s×1=\mathbasec{s}, (15b)

where for each we define:

##### Fast-slow coupling

equationparentequation

 A{f,s} = (16a) c{f,s} = (16b)

###### Proof

Replacing the -th continuous stage (6b) by the equivalent -th fast discrete sub-step to advance from to leads to: equationparentequation

 (17a) where V{f,i}k is the k-th stage of the fast discrete sub-step i. The corresponding slow stages are advanced as follows: (17b) Iterating after the slow stages (17b) yields: (17c)

where the last equation writes (17c) as a slow MRI-GARK stage (3a). From this we identify the slow method coefficients as:

 a{s,s}i,j \coloneqq ⎧⎪⎨⎪⎩0i=1,∑i−1p=max(j−1,1)b{f}Tγp,j(c{f})1≤j≤i,0i+1≤j≤s{s}, b{s}j \coloneqq s{s}∑p=jb{f}Tγp,j(c{f}).

Using (7) and the fact that for any due to the arbitrary accuracy of the fast scheme, we have that:

 b{f}Tγp,j(c{f})=∑k≥0γkp,jb{f}Tc{f}×k=∑k≥0γkp,jk+1=¯¯¯γp,j.

Consequently, the slow method coefficients satisfy:

 a{s,s}i,j = ⎧⎪⎨⎪⎩0i=1,∑i−1p=max(j−1,1)¯¯¯γp,j1≤j≤i, i≥2,0i≤j≤s{s},,b{s}j=s{s}∑p=j¯¯¯γp,j,

which establishes the equations (14).

We identify the slow-fast coupling coefficients as follows:

 a{s,f,λ}i,j \coloneqq {Δ\mathbasec{s}λb{f}j1≤λ≤i−1,0i≤λ≤s{s}, A{s,f,λ} = =Δ\mathbasec{s}λgλ+1b{f}T∈Rs{s}×s{f},λ=1,…,s{s},

which establishes (15a). Moreover:

 c{s,f} = s{s}−1∑λ=1A{s,f,λ}11{f}=s{s}−1∑λ=1Δ\mathbasec{s}λgλ+1, c{s,f}i = s{s}−1∑λ=1(\mathbasec{s}λ+1−\mathbasec{s}λ)(gλ+1)i=i−1∑λ=1(\mathbasec{s}λ+1−\mathbasec{s}λ)=\mathbasec{s}i,

which proves (15b).

Replacing (17c) into the discrete fast stage equations (17a) leads to:

 V{f,i}k = = yn+Hi∑λ=1s{f}∑j=1a{f,f,i,λ}k,jf{f}(V{f,λ}j)+Hi∑j=1a{f,s,i}k,jf{s}(Y{s}j).

where the last equation is the generic MR-GARK fast stage (3b). We identify the fast scheme coefficients as:

 A{f,f,i,λ} \coloneqq

which shows (13a), and

 (18)

which establishes (13c). Moreover:

 b{f,i}T =

which proves (13b).

For the fast-slow coupling we have that:

which establishes (16).

### 3.1 Order conditions

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

###### Definition 4 (Some B-series coefficients)

Consider the following Butcher tree :

 tk\coloneqq[τ,…,τk times]∈T,

where is the tree of order one and is the operation of joining subtrees by a root, and the following B-series coefficients of the exact solution of the fast sub-system: equationparentequation

 (19a) For a fast Runge-Kutta method of arbitrary accuracy it holds that: (19b)

Using (19) we define the following matrix:

 A{s,s}\coloneqq\mathbaseA{s}+∑k≥0ζkΓk. (20)

#### 3.1.1 Internal consistency

###### Theorem 1 (Internal consistency conditions)

The MRI-GARK method fulfills the “internal consistency” conditions:

 c{s,f}=c{s,s}≡\mathbasec{s}andc{f,f}=c{f,s} (21)

for any fast method iff the following conditions hold:

 Γ0⋅11{s}=Δ\mathbasec{s}andΓk⋅11{s}=0∀k≥1. (22)

###### Proof

From (14c) and (15b) we see that the first equation of (21) is automatically satisfied. Next, we write the second internal consistency equation (21) in the following equivalent form by equating (18) and (16b):

This relation has to be satisfied for any independently of the choice of the fast discretization method. In order to achieve this we equate separately different powers to obtain:

 (Γ011{s})i=Δ\mathbasec{s}i;(Γk11{s})i=0,  ∀k≥1;i=1,…,s{s},

which establishes (22).

###### Remark 4

Assuming that both the fast and the slow methods have order at least two, the internal consistency conditions (22) imply that the overall scheme is second order.

#### 3.1.2 Fast order conditions

Since the fast component is solved exactly (or equivalently, using a base Runge-Kutta scheme of arbitrary accuracy) the fast order conditions for any order are implicitly satisfied.

#### 3.1.3 Slow order conditions

The slow component of the MRI-GARK method (14) is the Runge-Kutta scheme . The MRI-GARK slow order conditions are fulfilled by selecting a slow base method of order .

###### Remark 5 (Coefficients of the slow base method)

The integrated gamma coefficients (7) satisfy:

 (23)

###### Remark 6 (Coefficients of the slow embedded method)

Consider now the embedded slow method , and the coefficients of corresponding MRI-GARK method given by the equations (14a) and (14b). From the structure of the error controller discussed in Remark 2 the coefficients of the main and the embedded methods differ only in the last row. Using (23), equation (14b) for the embedded method reads:

 ˆγki,j=γki,j  for  i=1,…,s{s}−1,ˆ¯¯¯γs{s},j=ˆb{s}j−a{s}s{s},j. (24)

###### Lemma 2 (Some useful formulas)

We start with several formulas that will prove useful in the derivation of order conditions. From (15a) and (13c) we have:

 (25)

From (16a):

 (26)

From (13a) and (13c) we have that:

 (27)

From (15a) it follows that:

 b{s}TA{s,f,i}=Δ\mathbasec{s}i(\mathbaseb{s}Tgi+1)b{f}T=Δ\mathbasec{s}i(s{s}∑ℓ=i+1\mathbaseb{s}ℓ)b{f}T. (28)

Finally, from (16a) we have that:

 b{f}TA{f,s}=Δ\mathbasec{s}TA{s,s}. (29)

#### 3.1.4 Third order coupling conditions

###### Theorem 2 (Third order coupling condition)

An internally consistent MRI-GARK method (6) has order three iff the slow base scheme has order at least three, and the following coupling condition holds:

 Δ\mathbasec{s}TA{s,s}\mathbasec{s}=16, (30)

where the coefficients are defined in (19), and in (20).

###### Proof

For internally consistent MR-GARK schemes there are two third order coupling conditions [21, 41]. We proceed with checking them. If the slow base method is of at least third order, then (25) implies that the first coupling condition is automatically satisfied:

 16 =

Using (29) the second coupling condition reads:

 16 = b{f}TA{f,s}\mathbasec{s}=Δ\mathbasec{s}TA{s,s}\mathbasec{s},

which establishes (30).

#### 3.1.5 Fourth order conditions

###### Theorem 3 (Fourth order coupling conditions)

An internally consistent MRI-GARK method (6) has order four iff the slow base scheme has order at least four, and the following coupling conditions hold: equationparentequation

 (31a) (31b) d{s}TA{s,s}\mathbasec{s}=124, (31c) (31d) (31e) with the coefficients ζk,ωk,ξk defined in (19), A{s,s} defined in (20), and: z{s} \coloneqq [\mathbasec{s}i+12−\mathbasec{s}i2]1≤i≤s{s}, (31f) d{s} \coloneqq (31g) t{s} \coloneqq (31h)

###### Proof

For internally consistent MR-GARK schemes there are ten coupling conditions for order four . We proceed with checking each of them.

Condition a.

 18 = b{f}T(c{f}×A{f,s}c{s}) = s{s}∑i=1Δ\mathbasec{s}ib{f}T(11{f}\mathbasec{s}i\mathbaseA{s}i,:\mathbasec{s}+∑k≥0(A{f,f}c{f}×k)c{s}iΓki,:\mathbasec{s} = s{s}∑i=1Δ\mathbasec{s}i(\mathbasec{s}i+1+\mathbasec{s}i)2\mathbaseA{s}i,:\mathbasec{s}

Using notation (31f) for the first term one obtains equation (31a).

Condition b. Using (25) and the fourth order of the slow base method one checks that the second MR-GARK coupling condition is satisfied automatically:

 18 =

Condition c. Using (29) the third coupling condition reads:

 112 = b{f}TA{f,s}\mathbasec{s}×2=Δ\mathbasec{s}TA{s,s}\mathbasec{s}×2,

which proves equation (31b).

Condition d. Using (15), (13c), and (28) the fourth coupling condition reads:

 112 = b{s}TA{s,f}c{f}×2 = =

We have that:

 (32)

for order four slow base methods with . Consequently, the fourth MR-GARK coupling condition is automatically satisfied.

Condition e. Using (25) one verifies that the fifth coupling condition always holds:

 124 =

Condition f. Using (28) and (26):

 124 = b{s}TA{s,f}A{f,s}\mathbasec{s} =

Using the notation (31g) and (20) the above establishes (31c).

Condition g. Using (27) and (28):

 124 = b{s}TA{s,f}A{f,f}c{f} = s{s}∑i=1⎛⎝Δ\mathbasec{s}i(s{s}∑ℓ=i+1\mathbaseb{s}ℓ)b{f}T⎞⎠⋅ ⋅(12\mathbasec{s}i211{f}+c{s}iΔ\mathbasec{s}ic{f}+Δ\mathbasec{s}i2A{f}c{f}) =

Using (32) we see that this condition is automatically fulfilled.

Condition h. Using (26):

 124 = b{f}TA{f,f}A{f,s}c{s} =

With notation (31h) this establishes (31d).

Condition i. Using (16a):

 124 = b{f}TA{f,s}A{s,s}c{s} = s{s}∑i=1Δ\mathbasec{s}i(\mathbaseA{s}i,:+∑k≥0ζkΓki,:)\mathbaseA{s}\mathbasec{s},

which proves (31e).

Condition j. Using (16a) and (29) the tenth coupling condition reads:

 124 = b{f}TA{f,s}A{s,f}c{f}=12Δ\mathbasec{s}TA{s,s}\mathbasec{s}×2,

which is the same as condition (31b).

## 4 Linear stability analysis

### 4.1 Scalar stability analysis

For additively partitioned systems (1) we consider the following scalar model problem:

 y′=λ{f}y+λ{s}y,λ{f},λ{s}∈C−. (33)

Let and . The stage computation (6b) applied to (33) solves exactly a linear ODE with the following analytical solution:

where the coefficients are functions of the fast variable:

and are defined using the following family of analytical functions:

 φ0(z)\coloneqqez,φk(z)\coloneqq∫10ez(1−t)tk−1dt,φk+1(z)=kφk(z)−1z. (34)

The MRI-GARK (9) applied to (33) reveals a stability function that depends on both slow and fast variables:

###### Definition 5 (Scalar stability)

The scalar slow stability region is defined as:

 (35)

Thus ensures that the solution is stable for any system (33) with in an -wedge in the left complex semi-plane, and of absolute value bounded by .

###### Example 2 (Stability functions of second order methods)

The stability function of the explicit midpoint MRI-GARK scheme Eq. 10 is a quadratic polynomial in , with coefficients analytical functions of :

Similarly, the stability function of the implicit trapezoidal MRI-GARK scheme of Eq. 11 is a rational function in , with coefficients analytical functions of :

### 4.2 Matrix stability analysis

Following Kværnø , for component partitioned systems (8) we consider the following model problem:

 (36)

where

The eigenvalues of

are linear combinations of the slow and fast diagonal terms, and . For the fast sub-system has a weak influence on the slow one; the first eigenvalue is slow and the second one is fast. For the slow sub-system has a weak influence on the fast one, and the first eigenvalue is fast.

###### Remark 7 (Scale considerations)

In order to have the slow and fast contributions to of similar magnitude, and the contributions to of similar magnitude as well, one needs a coupling coefficient inversely proportional to the scale ratio of the system, .

The general system (36) is complex, and in order to ease the analysis process some simplifications are needed. First, the factor represents a scaling of the fast variable, as can be seen from rewriting the system (36) in terms of the variables . To simplify the analysis we take , i.e., we assume that the multirate method is applied to the system (36) after rescaling the fast variables. It is possible, however, that the stability of some multirate schemes is affected by the scaling, in the sense that the same scheme, applied to systems with different scalings, has different stability properties.

In order to further simplify analysis we restrict ourselves to considering the case where