{'fraction_non_alphanumeric': 0.0406754008418862, 'fraction_numerical': 0.01558435415976919, 'mean_word_length': 4.514736567553469, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Relaxed Multirate Infinitesimal Step Methods for Initial-Value Problems

This work focuses on the construction of a new class of fourth-order accurate methods for multirate time evolution of systems of ordinary differential equations. We base our work on the Recursive Flux Splitting Multirate (RFSMR) version of the Multirate Infinitesimal Step (MIS) methods and use recent theoretical developments for Generalized Additive Runge-Kutta methods to propose our higher-order Relaxed Multirate Infinitesimal Step extensions. The resulting framework supports a range of attractive properties for multirate methods, including telescopic extensions, subcycling, embeddings for temporal error estimation, and support for changes to the fast/slow time-scale separation between steps, without requiring any sacrifices in linear stability. In addition to providing rigorous theoretical developments for these new methods, we provide numerical tests demonstrating convergence and efficiency on a suite of multirate test problems.


Introduction
Increasingly, computational science requires large-scale simulations that consistently and accurately couple distinct physical processes. While the mathematical models for specific processes often have a well-known type (hyperbolic, parabolic, etc.) and are suitable for classical numerical integrators, the same cannot be said for the coupled models. These multiphysics models are often of mixed type, may have limited differentiability, and involve processes that evolve at dissimilar rates. As such, many multiphysics simulations require more flexible time integrators that may be tuned for these complex problems.
In this work, we consider systems of autonomous ordinary differential equations that may be characterized by two distinct dynamical time scales, one "fast" and the other "slow", y(t 0 ) = y 0 ∈ R n .
We note that this additive form of the problem does not prohibit multirate applications where the solution variables are also partitioned between the two scales, In this work, we consider the development of multirate methods where the slow operator is integrated explicitly and the fast operator is integrated either explicitly or implicitly. To simplify our presentation, we denote the characteristic time-scales of the slow and fast components of the multirate problem as dt s and dt f , respectively. We define the time-scale separation between these as m = dt s /dt f . We select the time-step sizes that are used to integrate these components as dt s = h and dt f = h/m, where m is the smallest integer satisfying m < m. We note that other authors weaken this restriction, by allowing for variable step sizes at the fast time scale [11]. While the methods proposed here may be easily extended to such situations, the focus of this paper is on construction of higher-order methods; extensions to allow for adaptive fast substepping are left for future work.
As with single-rate problems, the order of accuracy of a multirate method describes how the error in the combined system behaves as h → 0. Within the multirate context, the overall error arises from multiple sources: the error in evolution of the fast components, the error in evolution of the slow components, and the coupling error between the two scales. For a variety of multirate methods, the fast component method is a composition of m steps of a fast base method, while the slow component method is one step of a slow base method -in this case the multirate method order depends on the order of both base methods and on the coupling order.
Multirate methods have been in use for some time. The first such methods were introduced by Gear and Wells in the early 1980's [12,13,14]; these methods considered integration of only a single time scale at a time, using linear interpolation from one scale to provide data for integrating the other. This work was extended by Constantinescu and Sandu to construct second order multirate methods based on explicit linear multistep base methods [15,16]. More recently, Fok and Rosales constructed linearly fourth order accurate multirate methods using higher-order interpolation between the slow and fast stages of a particular Runge-Kutta base method [17].
A more general approach to construction of third-order multirate methods was introduced in 1999 by Kvaernø [18], and was further elaborated by many authors [19,20,21,22]. In these approaches, one assumes a partitioned formulation with variables split into active (fast) and latent (slow). From here, higher-order interpolants between the two components are obtained by inter-leaving the interpolation processes between the fast and slow scales. A particular strength of this work is their derivation of a set of simplifying assumptions, so that once these are satisfied and the base method is third order, the resulting multirate method will also be third order.
An alternate approach for higher-order multirate schemes has been pursued through extrapolation of lower-order methods. The first such work is from Engstler and Lubich's 1997 paper [23], that generated higher-order solutions by extrapolating first-order methods. The efficiency of these methods was subsequently improved through use of dense output formulas [23,24,25]. More recent methods of this type, including explicit fast integration and explicit slow integration (explicit/explicit), as well as implicit fast integration and explicit slow integration (implicit/explicit), have been explored by Constantinescu and Sandu [26,27], that investigated extrapolation from first-order accurate methods. These multirate extrapolation approaches increase the possible order of the resulting method at the cost of repeated evaluations of the same time step. Hence, although these methods are reasonable for extrapolating up by one or two orders of accuracy, they may become cost-prohibitive when generating highorder methods.
Our work builds most closely off of recent work in Multirate Infinitesimal Step (MIS) methods. Originally developed by Wensch, Knoth and Galant, MIS methods were constructed as a generalization of split-explicit methods in the context of numerical weather prediction [28,29,8]. A key contribution from [8] was the development of a systematic approach to the order conditions for splitexplicit methods based on PRK theory [8], allowing the development of second and third order MIS methods for a variety of applications [6,7,8,1,2,3,4,5]. In particular, we highlight the the development of the Recursive Flux-Splitting Multirate (RFSMR) methods in [1] based on this MIS theory.
Coincident with the above developments on MIS and RFSMR methods, Sandu and Günther introduced a theory for Generalized-Structure Additively-Partitioned Runge-Kutta (GARK) methods [30,10]. As the name suggests, this theory may be applied to understand a wide range of Runge-Kutta-like time integration methods (including Kvaernø's methods, MIS and RFSMR), and notably provides both order conditions and a linear stability theory for methods in this general format. For simplicity, they formulate GARK methods for autonomous systems of ODEs written in additively-split form, Here, we consider the case of N = 2, corresponding to our target problems of the form (1), i.e.
GARK methods for this problem are then uniquely determined by coefficients in an expanded Butcher tableau, where internal stages and solution update are given by In particular, we note that this constitutes a "hybrid" of partitioned and additive Runge-Kutta formulations: there are distinct stages for each right-hand side component (PRK-like), and all right-hand side functions are used to update each stage (ARK-like). As is typical for additive Runge-Kutta methods, Sandu and Günther make the simplifying assumption of internal consistency, i.e.
Order conditions for methods up to fourth order (and for GARK tableau of arbitrary size) are provided in matrix-vector form in [30,10], and additionally in elementwise form in [30]. As our work extends the aforementioned MIS methods to fourth order, we reproduce the more concise matrix-vector form here, where we assume internal consistency (5): In these formulas, σ, ν and µ each range over the values {f, s}, and we denote a column vector of ones in R s {σ} as 1 {σ} . Throughout this manuscript we denote standard matrix and vector multiplication with adjacent objects (e.g., b c is an inner product), and we use the × operator or exponents to denote componentwise multiplication (e.g., . . , n). We note that of the above equations, (6) corresponds to the order 1 conditions, (7) corresponds to the order 2 conditions, (8)-(9) correspond to the order 3 conditions, and (10)- (13) correspond to the order 4 conditions. Hence for a two-component splitting, a fourth order method requires 28 order conditions to be met. As in [31], we will categorize these order conditions into two groups: all conditions that include b {f } (i.e., the 14 conditions (6)- (13) with σ = f ) are considered as "fast order conditions", and all conditions including b {s} (i.e., the 14 conditions (6)-(13) with σ = s) are considered as "slow order conditions."

GARK representation of MIS methods
Prior to introducing our proposed RMIS methods, we first summarize the analysis of MIS methods using the GARK formalism, first shown in [9,Theorem 4]. In that analysis the problem is considered to be in autonomous form in order to simplify the exposition; however, we note that in [5] the MIS approach is extended to non-autonomous problems.
MIS methods are typically constructed using a pair of "base" Runge-Kutta methods: a variation of the "outer" method A O , c O , b O is applied to the slow time scale, and a variation of the "inner" method A I , c I , b I is applied to the fast time scale. The outer method is an explicit Runge-Kutta method with s O stages, with the requirement that The MIS scheme makes no assumption about the structure of the inner method, aside from the requirement that it is a one-step method. In general, MIS methods may be run so that the inner method takes several "fast" time steps to evolve between the slow stages, e.g. over the interval t n + c O i−1 h, t n + c O i h , i = 2, . . . , s O . As in [9], both here and in our new theoretical work in Section 2, we consider the case when only a single step of the inner method is taken to evolve between each slow stage. Since any sequence of m steps of a q-order, one-step method may be equivalently written as a single step of a corresponding q-order method, this in no way limits the applicability of the current theory. We do note, however, that if the abcissae of the outer method are not evenly spaced, i.e.
then the Butcher tables for the inner one-step methods corresponding to the outer stages i and j will not be identical.
As shown in [9,Theorem 4], the GARK tables corresponding to the MIS scheme are given by A where the block matrices and vectors comprising these are defined using standard matrix-vector notation, We note that for MIS methods the number of slow stages matches the outer table, s s = s O , and the number of fast stages equals the product of the outer and inner tables, s f = s O s I . We also note the slight difference in presentation of A {f,s} from that shown in [9,Theorem 4]; in equation (16) above the generic entry denoted by i corresponds to the i-th block row of the matrix. Based on these Butcher tables for the GARK coefficients, Günther and Sandu prove a number of particularly beneficial properties of MIS methods [9,32]: (i) The coefficients (14)-(18) satisfy the simplifying internal consistency conditions (5).
(ii) If both the fast and slow methods have order at least two, then the overall multirate method is second order.
(iii) If both the fast and slow methods have order at least three, and if the outer method satisfies the additional condition then the overall MIS is third order (this guarantees that the MIS method satisfies the third-order "fast" coupling condition, b {f } A {f,s} c {s} = 1 6 ).

Relaxed Multirate Infinitesimal Step Methods
In this work, we extend the MIS schemes to allow for methods with a more 'relaxed' definition of b {f } . All other components remain identical to those in MIS methods, namely the construction of A {f,f } , A {f,s} , A {s,f } , A {s,s} and b {s} . As with MIS methods, we consider a pair of base methods, T I and T O , corresponding to the inner and outer base Runge-Kutta tables, having s I and s O stages, respectively, where T O is explicit and T I may be either explicit or implicit; however, we require an additional assumption that T O has an explicit first stage. Due to this structural similarity between RMIS and MIS methods, we may immediately leverage all aspects of the theory presented in [9,32]  (ii) If both T I and T O have order at least two, then the RMIS method satisfies the "slow" second-order conditions (i.e. (6)-(7) with σ = s).
(iii) If both T I and T O have order at least three, then the RMIS method satisfies the "slow" third-order conditions (i.e. (6)-(9) with σ = s and ν = {f, s}).
The overall order of accuracy for a RMIS multirate method will therefore depend both on the base methods and on the choice of coefficients b {f } . To this end, we first extend the analysis from Günther and Sandu [9,32] to consider the "slow" fourth-order conditions for MIS methods.
Theorem 1. Assume that the inner base method T I is at least third order, the outer base method T O is at least fourth order, and that both satisfy the rowsum consistency conditions (5). If T O is explicit and satisfies the additional where
We examine these in order: where we have used the fact that T I is third order, and T O is both fourth order and explicit.
Our final slow fourth-order condition becomes where we have relied on the assumption (20) and the fact that T I is at least second-order accurate. 2 The above result, in combination with [32, Theorem 3.1], guarantees that when using any third-order inner method T I , and any explicit fourth order outer method T O that satisfies (20), the MIS and RMIS methods will automatically satisfy all of the "slow" fourth-order order conditions. We now turn our attention to the fast solution coefficients, b {f } ∈ R s f . Assuming that we select T O and T I according to the above criteria, then we must only select b {f } to satisfy the 14 "fast" order condition equations, where ν, µ = {f, s}, to guarantee that the overall RMIS method is fourthorder accurate. We note that each of these equations depends linearly on b {f } , indicating that if s f = s O s I is "large enough", we may be able to solve this under-determined linear system for valid sets of coefficients b {f } . We further note that for the RMIS method to be truly "multirate", s f ∝ m s s , indicating that as the problem becomes more multirate (and m increases), the fourth-order conditions become even simpler to satisfy. Of the potentially infinite options for b {f } that satisfy the above criteria, we select the choice We note that this particular structure allows for us to simplify the order conditions (21)-(28) dramatically, as seen in the following result.
Lemma 1. Suppose that the coefficients b {f } are chosen as in equation (29), and that the inner Butcher table T I has explicit first stage (i.e. the first entry of c I and the first row of A I are identically zero). Then the following identities hold: where These follow from direct application of the vector-matrix products: Proof of the results (33) and (34) are essentially identical to those above.  16) and (29) is fourth-order accurate.
We may immediately use Theorem 1 to ensure that all of the "slow" fourthorder conditions are met. We also note that due to Lemma 1, all of the "fast" fourth-order conditions are equivalent to "slow" order conditions, which are al-ready satisfied by Theorem 1: We note that the above theorem does not include the requirement (19) on T O for third-order accuracy of MIS methods, due to the alternate structure of b {f } within RMIS methods. We further note that based on the proof above, if all of the assumptions from Theorem 2 are satisfied except for (20), then the RMIS method will be third-order accurate.

RMIS with MIS embedding
Since RMIS methods and MIS methods share the same algorithmic structure, the same coefficients and only differ in their selection of b {f } , one may naturally question whether an MIS method could be used as an embedding within a RMIS method to enable temporal error estimation. For this combination of RMIS and MIS methods to work, however, the base methods T O and T I must be compatible with both approaches. Specifically, if these base methods satisfy the requirements: (a) T O is explicit and at least order four, (b) T I has explicit first stage and is at least order three,  (29) and MIS coefficients will result in a fourth-order method with third-order embedding. While the above criteria may be met by a variety of base methods, we have identified some specific candidates for use in the numerical results in Section 4. To this end, we relied on Butcher's derivation of families of explicit 4th order base Runge-Kutta methods [33] when determining T O . His general solution for a 4-stage fourth order explicit Runge-Kutta method depends only on two free variables, c 2 and c 3 , and is given by , For T O with this structure, the 3rd-order MIS condition (19) is equivalent to (37) and the 4th-order RMIS condition (20) is equivalent to We plot the solutions (c 2 , c 3 ) ∈ [0, 1] 2 to these equations in Figure 1, and point out that the intersection of these curves denote choices of T O that satisfy all  (37) and (38). Values of c 2 and c 3 not on these curves result in MIS and RMIS methods of 2nd and 3rd order, respectively.
of our desired criteria (a)-(d). Of these, only two satisfy the MIS criteria that c 2 < c 3 : Both due to its historical elegance, as well as its equally-spaced c i , we use this base method for T O in our numerical results that follow.

Method comparisons
In the sections that follow, we examine the linear stability, and assess the attainable order of accuracy, for both the proposed RMIS and MIS algorithms. To this end, we examine subcycled, "telescopic" versions of both algorithms, wherein the inner table T I is comprised of ∼ m/s O substeps of the outer table, T O , thereby facilitating further recursion to support problems with three or more rates. We note that due to the patterned structure of both the MIS and RMIS methods, these subcycled algorithms may be implemented more efficiently than a generic GARK scheme having ms fast and s slow stages. Discussion on implementation-specific details for these algorithms is provided in Appendix A.
For these methods, we utilize T O as either the four-stage 3/8-Rule (41) that satisfies the 3rd-order MIS condition (19), is widely used as a basis for MIS methods in the literature [5,7,1,8], and in previous tests of all multirate methods listed in Section 1 we found it to be the most efficient (smallest error with least computational cost).
To be more precise regarding how we form T I from these two choices of T O to enable subcycling, we elaborate on the case of time-scale separation factor m = 100. Schlegel suggested using to measure the appropriate number of subcycles for each subcycling period [3]. This measure was used as a guide for selecting how many subcycles should be used for our methods, although we determined our comparisons would be fairest if the ratio of fast function evaluations to slow function evaluations was approximately constant.
The 3/8-Rule (41) has four evenly-spaced abcissae with final c O s O = 1. Hence the multirate implementation requires only three subcycling periods, and so we form T I using 34 subcycles of the T O table. The KW3 method has only three abcissae, but c O s O = 3 4 < 1, so although it has one fewer stage than the 3/8-Rule, it requires the same number of subcycling periods. Moreover, since these are unevenly spaced, our implementation could either use different A I per outer stage (to attempt nearly-identical substep sizes) or use a single A I for each outer stage, mapped to the largest subcycling interval. Choosing the latter for simplicity, we form T I using 35 subcycles of T O . This is a compromise between using the largest subcycling interval for substep sizes and making the work of fast functions vs slow functions comparable.
In addition to the RMIS and MIS algorithms, for some of the convergence tests we also utilize an "optimized" method, constructed using the approach outlined after equations (21)-(28) above. More precisely, starting with the 3/8-Rule for T O , and choosing a time-scale separation of m = 100, we solved the under-determined linear system (21)-(28) for a valid set of coefficients b {f } ∈ R 408 . With the resulting free variables, we minimized the error in the fifth-order conditions using the two-norm of the residual of all the fifth-order conditions, using MATLAB's fminsearch algorithm, which is based on a simplex search method [35].
In the ensuing sections we therefore compare five multirate methods: • "Opt-3/8" -this uses T O given by (41) and constructs b {f } via the optimization approach described above (should be O(h 4 ) accurate); • "RMIS-3/8" -this is our proposed RMIS method using T O given by (41) that satisfies the 4th-order condition (20) (should be O(h 4 ) accurate); • "RMIS-KW3" -this is our proposed RMIS method using T O given by (42) that does not satisfy (20) (should be O(h 3 ) accurate); • "MIS-3/8" -this is the MIS method using T O given by (41)  We note that although we have proposed an embedding of MIS within RMIS, optimal approaches for time-step adaptivity (varying both h and m) are not within the scope of this manuscript and are left for future work. Furthermore, of the five methods listed above, only one provides such an embedding for temporal error estimation. We therefore explore each of the above methods using fixed time-step sizes h.

RMIS and MIS linear stability analysis
In their paper introducing GARK methods [10], Sandu and Günther analyze linear stability using a modification of the scalar-valued Dahlquist test problem, where the real parts of both λ {f } and λ {s} must be negative. While that approach yields elegant definitions of the resulting GARK stability regions, we prefer an approach for multirate linear stability analysis that was first proposed by Kvaernø [22] in 2000, and has subsequently been used by a variety of authors [27,36,37,38,39]. In this approach, one instead considers the partitioned test problem y f y s = g 11 g 12 g 21 g 22 Here, y f and y s correspond to the fast and slow variables, respectively, with a resulting splitting of the right-hand side into fast and slow components as shown in (2). The benefit of this approach is that unlike the purely additive scalar version, this directly allows direct analysis of stability as a function of the strength of fast/slow coupling in the problem. Defining the matrix Z = hG, stability is defined as usual by first determining the amplification matrix S(Z) for one step of the numerical method, i.e.
The multirate method is then linearly stable stable if the eigenvalues of S have magnitude less than 1 for a given step size Z = hG. While more complex than standard linear stability regions (seemingly depends on five parameters, h, g 11 , g 12 , g 21 and g 22 instead of only two, h and λ), we may in fact determine the eigenvalues of S using only 3 real parameters, ξ, η, and κ [36,37,39]. The eigenvalues of G have negative real part if g 11 , g 22 < 0 and β = g12g21 g11g22 < 1, where here β is a measure of the off-diagonal coupling strength of the problem. We may then define the parameters where κ encodes the time-scale separation of the problem, ξ encodes the stiffness of the fast time scale (as hg 11 goes from 0 → −∞, ξ moves from 0 → −1), and η encodes the strength of coupling (no coupling at η = 0, increased coupling as |η| → 1, coupling-dominant as η → −1). With these three parameters, one may plot snapshots of the (ξ, η) stability regions for fixed problem time-scale separation values κ.
In the stability region plots that follow, we match the problem and method time-scale separation values, i.e. κ = m. We then create linearly-spaced arrays of 100 ξ values in (−1, 0) and 200 η values in (−1, 1). For each of the 20000 resulting (ξ i , η j ) combinations, we compute the eigenvalues of S(Z) to determine stability, and mark the stable regions in yellow yellow and unstable regions in blue.
In Figure 2 we plot the stability regions for the RMIS-3/8 and MIS-3/8 methods for the time-scale separation values κ = m = {10, 100}. We see that the stability regions for both methods are similar: for weakly coupled fast and slow time scales, these indicate stability for ξ values close to zero and instability for ξ close to -1, as would be expected for any explicit method. Additionally, both exhibit a region of increased stability for problems with stronger fast/slow coupling, followed by complete instability as η → −1 (where the coupling terms become infinitely large). Comparing the plots of κ = 10 versus 100, we see that the region of increased stability shifts toward stronger fast/slow coupling as the problem time-scale separation increases. Lastly, as expected for explicit methods, the maximum step size shrinks as the fast time scale decreases.
To see whether these stability properties translate to other RMIS and MIS methods (as opposed to only the ones based on the 3/8-Rule, above), we additionally examine the linear stability of the RMIS and MIS methods that use the other 4-stage base method T 0 satisfying the 3rd-order MIS and 4th-order RMIS conditions, given by the intersection point (40) from Figure 1. These plots are shown in Figure 3. While the precise shapes of these regions have changed slightly from those in Figure 2, the properties observed above are indeed retained for this alternate base method. As a final stability comparison between MIS and RMIS methods, we investigated the following question: given that there exist one-parameter families of 3rd-order MIS methods and 4th-order RMIS methods (as seen in Figure 1), can we find the MIS and RMIS methods with 'largest' linear stability region, and how would these compare against one another? To this end, we tested 100 base methods along this one-parameter family for each method and selected the one with largest stability region. The results of this investigation for κ = m = 10 are shown in Figure 4. We again note the similar stability regions as before, showing slight variations in the stability boundaries but the same overall shapes.
Based on the results shown in this section, we conclude that the proposed RMIS methods suffer no deterioration in stability as compared with their MIS 'cousins'. In the next section we investigate their convergence and efficiency on a variety of standard multi-rate test problems.

Convergence and Efficiency Tests
In this section we provide numerical results to verify the analytical properties of the proposed RMIS methods and compare their efficiency against the MIS methods of [1]. To that end, we consider three standard multirate test problems, corresponding to the invertor-chain test problem [18,21,19,27], a linear multirate test problem with strong fast/slow coupling from Kuhn and Lang [39], and the well-known brusselator test problem. For each problem we either compare against the analytical solution or a high-accuracy reference solution. These reference solutions were generated by using the 12th order Gauss implicit Runge-Kutta method, using fixed time-steps which are 4 times smaller than the smallest h value tested for our multirate methods.
We measure solution error in each test by computing the root-mean-square error over all solution components and over all time steps t n = t 0 + hn, where y n,h and y n,ref are our computed and reference solutions in R n at t n , respectively, · 2 is the standard vector 2-norm, and M = (t F − t 0 )/h is the total number of overall time steps of size h taken in each run over t ∈ [t 0 , t F ].
For each test, we present both plots of RMSerror versus h, as well as tables of the observed numerical order based on a least-squares fit these data where outliers are ignored, i.e., the numerical orders of convergence only include data where 10 −9 ≤ RMSerror ≤ 1.
To assess the efficiency of each method, we utilize the standard error-versuscost plots. Since these simulations were performed in MATLAB, where timings can provide rather poor predictions of runtimes for true HPC applications, in these tests we measure cost by counting the total number of ODE right-hand side function calls. For the subcycled and telescopic RMIS and MIS methods examined here, this may be easily calculated as where n i is the number of fast step subcycles required per slow stage i. In the context of our five methods (Opt-3   versus TotalFunctionCalls; "efficient" methods correspond to curves that are nearest the bottom-left corner.

Inverter-chain
The inverter-chain problem is a partitioned multirate ODE system that models a chain of MOSFET inverters, that has been used for testing multirate ODE solvers throughout the literature [21,19,36,40,41,38,42,43,44,1,27]. The form of the model that we examine is primarily based on the version used by Kvaernø and Rentrop [18], although we utilize an additional scaling term as used in [19,27]. The mathematical model is given by the system of ODEs for y(t) ∈ R n I , 0 ≤ t ≤ 7: where g k (t, y) = g(y in (t), y 1 (t), y 0 ), where y op = 5 V, y T = 1 V, y 0 = 0 V, γ is the scaling term for tuning the time-scale separation of the problem (we use γ = 100), and y in (t) is the forcing function, We note that y in (t) causes the problem (49) to be non-autonomous. This is easily handled since the MIS and RMIS methods are internally consistent, so we may identify 'stage times' t The time scales for these inverters decrease with index, so we partition this so that the first b equations are "fast" and the remainder are "slow", In our tests, we use n I = 100 total inverters, with the first b = 3 grouped into f {f } (t, y). We note that in other work using a similar setup, the first b = 20 were chosen as "fast"; however in our tests we obtained improved error and stability when using the smaller value of 3. In Figure 5 we show the solutions for this problem, as well as a zoom-in of the initial departure of the fast inverters from the larger group. In our tests we found that accuracy in these initial departures proved crucial for overall solution accuracy.
In the left portion of Figure 6 we plot the RMSerror versus step size h for each of the five methods tested. We highlight a few observations in these results. First, all methods demonstrate convergence as h → 0 to a point, beyond which convergence stagnates. This stagnation point is below 10 −9 for all but the Opt-3/8 method, that stagnates slightly earlier at around 10 −8 , indicating a reference solution accurate to approximately 10 −9 for this test. We hypothesize that the Opt-3/8 method is more susceptible to accumulation of floating-point round-off than the other methods, since unlike the others that are defined by small-valued coefficients, Opt-3/8 has more widely-varying coefficients b . We note that all show their expected rates of convergence except for Opt-3/8, which is likely due to its larger error floor, since prior to that point it is converging at least as rapidly as the 3rd-order methods. Accuracy alone provides only an incomplete picture of performance, since the methods using the 3/8-Rule require 25% more function calls per step than those using KW3 To this end, in the right portion of Figure 6 we also plot the RMSerror versus TotalFunctionCalls for each of the five methods. Although the blue and magenta curves for the KW3-based methods indeed shift further to the left in relation to the other methods, RMIS-3/8 is still the most efficient of all the methods at nearly all error values, and is only outperformed by Opt-3/8 for relatively loose error values (above ∼ 10 −4 ). This is most critical for errors below ∼ 10 −8 , where RMIS-3/8 requires ∼ 10 times less work than the other methods.

Strongly-coupled Linear Test
As a second test, we use a linear ODE system with strong fast/slow coupling that was used by Kuhn and Lang in their studies of multirate stability [39]; variants of this problem have been used by a variety of authors in testing multirate algorithms [22,39,27]. Our motivation for this test is to more rigorously explore accuracy and efficiency of the RMIS and MIS algorithms in the face of strongly-coupled problems, thereby exercising the coupling matrices A {f,s} and A {s,f } , and exploring whether the extreme sparsity of our b {f } in the RMIS algorithm causes trouble for strongly-coupled problems. The ODE system is identical to the partitioned test problem (43) examined in Section 3, where here the system-defining matrix is given by and the problem is evolved over the time interval t ∈ [0, 1]. The eigenvalues of G are complex conjugates, giving rise to the analytical solution which is shown in Figure 7. The time-scale separation in this problem arises from the strong fast/slow coupling, rendering the 'sin' term in y 1 approximately 100 times stronger than in y 2 , at least until the solutions decay to zero. We cast this problem into additive multirate form by again partitioning the right-hand side into fast and slow components, In Figure 8 we again plot both RMSerror versus h (left) and RMSerror versus TotalFunctionCalls (right). This problem exhibits very straight order of accuracy lines for all methods except Opt-3/8, which begins to level off at approximately 10 −11 , again likely due to increased sensitivity to accumulation of floating-point round-off errors. The best-fit convergence orders for these results are . We do not fully understand the higher-than-expected orders of accuracy for Opt-3/8 and RMIS-3/8, except that since the problem is linear then many of the higher order conditions are trivially met. As with the previous problem, the reduced cost per step of the KW3-based methods is insufficient to surpass the efficiency of the 4th-order methods, with Opt-3/8 and RMIS-3/8 demonstrating essentially-identical efficiency for errors larger than 10 −10 , and with RMIS-3/8 proving more efficient past that point.

Brusselator
We consider a system of stiff nonlinear ODEs that captures some of the physical challenges of the brusselator chemical reaction network problem, first described as a 1D PDE by Prigogine in 1967 [45]. More recent usages of this test problem have been seen in a multiphysics paper by Estep in 2008, and a general computational multiphysics review paper in 2013 [46,47]. Our version of the problem is a tunable two-rate initial-value problem represented as a system of three nonlinearly-coupled ODEs, where the parameters are chosen to be a = 1.2, b = 2.5, and ε = 10 −2 . As shown in the above partitioning, the fast function f {f } contains only the term which is scaled by ε; the problem time-scale separation is approximately 1/ε = 10 2 = 100. With this particular choice of parameters and initial conditions, the problem exhibits a rapid change in the solution at the start of the simulation for t < 0.2, with slower variation for the remainder of the time interval, as shown in Figure 9. In Figure 10 we plot the RMSerror versus h (left) and versus TotalFunctionCalls (right). These plots show perhaps the most interesting results of the three tests. Focusing first on the convergence plots we note an intersection point in the curves, showing that the optimal choice of method can depend more intimately on the desired solution error. Beginning at the large-error end, we see that the MIS-KW3 method has smallest error at the largest tested step sizes, although those errors are not significantly different than those for MIS-3/8 or RMIS-3/8. At smaller h and smaller error values, the increased accuracy of RMIS-3/8 and Opt-3/8 appear, rapidly achieving errors ∼ 100 times smaller than the 3rd-order methods at step sizes ∼ 2 × 10 −3 . Convergence of all methods stagnates at approximately 10 −9 , clearly indicating the accuracy of our reference solution. The corresponding best-fit estimates of the convergence orders for these results are   When examining the efficiency plot at the right of Figure 10, we observe that the decreased cost per step of the KW3-based methods shift those results to the left; making MIS-KW3 the clear winner for larger error values. However, at error values below ∼ 10 −5 , the higher order methods begin to outperform the lower order methods, with Opt-3/8 the most efficient, followed by RMIS-3/8.

Conclusions
In this work, we propose a variation of the existing MIS/RFSMR multirate methods [6,7,8,1,2,3,4,5]. The MIS methods that we extend demonstrate a number of very attractive properties for multirate integration. These methods are telescopic, allowing for recursion to any number of problem time scales, allow subcycling of the fast method within the slow, and are highly flexible in that they allow for varying the time-scale separation of the method (m) between steps. Finally, when constructed using inner and outer Runge-Kutta methods, T I and T O , of at least third order, MIS methods are at worst second order accurate and achieve third order when T O satisfies the auxiliary condition (19).
Recent theory by Sandu and Günther complements these MIS methods nicely. Specifically, they have proposed a general formulation for analyzing a wide range of Runge-Kutta-like methods, named "Generalized-structure Additive Runge-Kutta" (GARK) methods [10], that lays a strong theoretical foundation for understanding and extending MIS methods. More beneficial, however, is their subsequent work that directly analyzed MIS methods using their GARK framework [31], opening the door for subsequent extensions to this analysis.
In this context, we propose a new multirate algorithm which we name "relaxed multirate infinitesimal step" (RMIS), due to a relaxation from the MIS approach on how the fast stage solutions are combined to construct the overall multirate step solution (i.e., the coefficients b {f } ). This simple change, while leaving the remainder of the MIS approach intact (i.e., the algorithmic approach as well as the coefficients A {f,f } , A {f,s} , A {s,f } , A {s,s} and b {s} ), allows for a number of remarkable extensions to their work. First, we are able to construct up to fourth order, and at least third order, multirate methods with comparable stability and improved efficiency to MIS methods, where the determining factors for order four are: T I must be at least third order and have explicit first stage; T O must be explicit, at least fourth order, and must satisfy the condition (20). A key component of this analysis is our result from Lemma 1, where we show that the fast order conditions are equivalent to the slow conditions, so one only needs to prove half of the conditions to guarantee overall order of accuracy. These new methods show nearly identical stability properties, and retain all of the above "attractive properties" described above for MIS methods (telescopic, subcycling-ready, flexible to adjusting m). In addition, since MIS and RMIS only differ in their selection of the b {f } coefficients that form the time step solution from the stage computations, then it is straightforward to define MIS as an embedding within an RMIS method.
In addition to providing detailed proof of the fourth order conditions for the RMIS algorithm, and comparison of the linear stability between MIS and RMIS methods, we provide numerical comparisons of the performance of multiple RMIS and MIS methods on three standard multirate test problems in Section 4. These problems include the standard "inverter-chain" system of nonlinearly coupled ODEs, a linear multirate problem with strong fast/slow coupling, and the standard "brusselator" stiff nonlinear ODE system. Through these experiments, we note that our theoretical expectations regarding the order of accuracy for each method are borne out in the results for each problem, namely that the RMIS methods satisfying (20) exhibit fourth order convergence, and those that do not are only third order accurate, and that the MIS methods are at best third order accurate (even when using higher-order base methods). As a result, due to its consistently strong performance (stability, error and efficiency), the proposed RMIS-3/8 method is a clear contribution to the ever-growing 'stable' of multirate methods, particularly when solutions with significant accuracy are desired.
We note that there are many avenues for further extensions of this work.
On a theoretical level, we have yet to explore the 5th-order conditions for RMIS methods, an effort that will undoubtedly leverage Lemma 1 repeatedly, and likely result in a set of additional conditions on the outer table T O . Similarly, we note that although all of the theory presented in Sections 1.1 and 2 assume an autonomous ODE, the inverter-chain problem is non-autonomous and still shows the predicted orders of accuracy for all methods tested, hinting that this theory can be extended to the non-autonomous context as well. On a purely computational level, we plan to investigate mixed implicit/explicit multirate methods based on the RMIS structure (implicit T I with EDIRK or ESDIRK structure), including selection of T O and T I pairs for optimal efficiency. There are also numerous algorithmic extensions of this work, including (a) exploration of the RMIS/MIS embedding to perform adaptivity (in both h and m) for efficient, tolerance-based, calculations; (b) utilization of the telescopic property for three-rate problems with large timescale separations, and for n-rate problems with smaller timescale separations (as arise in explicit methods for hyperbolic PDEs posed on spatially adaptive meshes); and (c) exploration of RMIS-based "self-adjusting" methods [36,37,38,43,17], that use the temporal error estimate to automatically determine a fast/slow partitioning for the problem.