The Interaction of Iteration Error and Stability for Linear Partial Differential Equations Coupled through an Interface

We investigate properties of algorithms that are used to solve coupled evolutionary partial differential equations posed on neighboring, nonoverlapping domains, where the solutions are coupled by continuity of state and normal flux through a shared boundary. The algorithms considered are based on the widely used approach of iteratively exchanging boundary condition data on the shared boundary at each time step. There exists a significant and sophisticated numerical analysis of such methods. However, computations for practical applications are often carried out under conditions under which it is unclear if rigorous results apply while relatively few iterations are used per time step. To examine this situation, we derive exact matrix expressions for the propagation of the error due to incomplete iteration that can be readily evaluated for specific discretization parameters. Using the formulas, we show that the universal validity of several tenants of the practitioner’s conventional wisdom are not universally valid.


Introduction
The class of multiphysics problems in which one physical process operates in one domain while a second physical process operates in a neighboring domain, and the solutions in the component subdomains are coupled by continuity of state and normal flux through a common boundary, which is central to a number of applications [1].Examples include conjugate heat transfer between a fluid and solid [2][3][4][5], Stokes-Darcy flow in karst aquifers [6] and in catalytic converters [7], modeling of the core and edge plasma flows in a tokamak reactor [8,9], and flow in heterogeneous porous media [10][11][12][13].
In applications, the processes in the component subdomains are often themselves represented by complex, multiscale, and multiphysics models and the component models are solved with different discretization methods utilizing significantly different scales.Such coupled problems present enormous challenges for efficient high performance computing [1].There are very strong incentives to use existing "legacy" codes for the component models and to treat component physics solvers as "black boxes." For these reasons the most commonly encountered solution strategy is a simple iterative approach [1][2][3][4][14][15][16][17].In this approach, the models in each of the component subdomains are associated with one of the two boundary coupling conditions and subsequently solved independently except for data from the coupling through the boundary.The coupling data for the boundary conditions for one component model are provided by the solution of the other component model from the previous iteration.In a time dependent problem, this iteration is carried out on each time step.
This type of coupled physics problem has been extensively studied by the numerical analysis community, with the goals of deriving accurate, stable numerical methods, efficient, and accurate coupling algorithms and deriving rigorous a priori analyses, for example, in [10-15, 17, 18], as well as a posteriori analyses [2][3][4]19].However, in practice, as mentioned, the component physics is usually complex and there are large differences in the discretization strategies employed in 2 Advances in Mathematical Physics the component subdomains necessitating significant processing of transferred information.Thus, the application of current rigorous mathematical methods and analysis is problematic if not impossible.This is partially the reason that there is a very signicant gulf between the kinds of practices carried out in high performance computational science and the "best practices" determined by the mathematical community [3].
In addition, the computational strategies employed in practice are often rationalized by a body of conventional wisdom which asserts that (1) stability is equivalent to accuracy and (2) the use of unconditionally stable implicit solvers for the component physics generally stabilizes the entire coupled simulation as long as the component solutions are stable.This is often rationalized by some sort of (Neumann) linear stability analysis applied to the simplied case of two coupled linear parabolic problems.Based on these ideas, simulations are deemed acceptable provided they do not exhibit obvious instabilities such as unbounded behavior or rapid oscillation.Moreover, due to the prohibitive computational cost, such conclusions are often based on the single computation at hand, rather than through an exhaustive study.
This paper attempts to address the difficulties that face extending rigorous convergence analysis to the complex models and discretizations that occur in practice and illustrate issues arising from the conventional wisdom.Instead of focusing on deriving conditions that guarantee good behavior as is usual for a standard convergence analysis, we adopt a different approach and develop a rigorous computable error representation that can be evaluated for any given choice of discretization parameters, thus allowing the conventional wisdom to be verified or not.We consider the canonical problem of two linear coupled parabolic equations and formulate the iterative solution method as a fixed point iteration for a block system on the discrete level.We then derive an exact formula for the error in the iterated solution.This formula is related to the Neumann series for the inverse of the full system matrix.The argument has several virtues: it is elementary, it subsumes any sort of ad hoc linear stability analysis, and it is general in the sense that it holds for a variety of discretization methods and range of scales for the component physics.Importantly, it allows for easy evaluation of various discretization parameter choices, for example, step size, space mesh, and number of iterations per time step.
We then present a detailed study of the canonical problem of two linear coupled parabolic equations that amply demonstrates shortcomings of the conventional wisdom.Firstly, we are able to dispel the notion that stability implies accuracy.In particular, we demonstrate that a divergent iteration can be part of a stable algorithm if the number of iterations used is low, while the resulting approximation is inaccurate.In this case, the user might not be aware that the iteration is divergent if they only consider whether or not obvious instability occurs.This case is particularly interesting since a serious problem is being masked by seemingly "reasonable" results.Secondly, we demonstrate that there is no guarantee that an "unconditionally stable" time integration scheme like backward Euler remains unconditionally stable if the system is not solved exactly at each time step.We explore cases where using a limited number of iterations leads to a conditional stability, which depends on the size of the time step.The values of time steps which provide stability can occupy a range with a maximum and minimum value, rather than just a maximum.The remainder of this paper is organized as follows.Section 2 introduces the problem and the notation associated with discretization and iteration.Section 3 derives the primary results of the paper regarding the stability of the iterative solution.Section 4 provides numerical examples.Section 5 addresses the multirate case, in which the time step in one component is an integer multiple of the time step in the other component.Multirate numerical examples are included.A brief conclusion is given in Section 6.

Problem Formulation and Discretization
We consider a system posed on a domain consisting of two neighboring, nonoverlapping, convex, and polygonal subdomains Ω 1 and Ω 2 in 1, 2, or 3 spatial dimensions that share a common polygonal interface boundary Γ (see [1,2,8,9] for applications).We illustrate in Figure 1.
The general form of the PDE system is together with the following additional boundary and initial conditions: (i) boundary conditions for U 1 on Ω 1 \ Γ and initial conditions for where and n is a unit vector normal to Γ and pointing from Ω 1 to Ω 2 .We next introduce notation for the discretized problem.The analysis is based on the backward Euler method in time but can accommodate various spatial discretization schemes [3].For the numerical results included later in the paper, we use a finite volume method in space [13,18,19], and the details of that discretization are given in Appendix A. In order to complete the boundary data for each subdomain, we adopt the common strategy of using the solution on the first subdomain to determine Dirichlet boundary data on Γ for the problem on the second subdomain and, likewise, use the solution on the second subdomain to provide Neumann boundary conditions for the problem on the first subdomain; for example, see [2-4, 8-12, 14, 15, 17, 19].Thus, each problem is provided with a full set of boundary data, but the subdomains are not treated symmetrically.By convention, we assume that component 2 provides state data for component 1, and component 1 provides flux data for component 2. This leads to the system of discrete equations: where the superscript is the time step index.Note that the form of the operators  21 and  12 depends on the choice of discretization.Some discretizations, for example, finite volume, require forming these operators using some combination of averaging, extrapolation, and interpolation [1].Other discretizations provide natural definitions.For example, in the case of the mortar method [18,19], the coupling between the subdomains has a formal structure based on Lagrange multiplier variables on the interface.In this case  21 and  12 correspond to the off-diagonal blocks of the Schur complement [20] that are formed by eliminating the Lagrange multiplier variables.
In all of the above cases, we now write (2) as a single set of linear equations.Namely, where 2.1.Block Iteration.As mentioned, the implicitly coupled system (3)-( 4) is not formed or solved exactly in practice.
The common approach is to use an iterative block solution approach that involves a sequence of solutions of each component problem with alternating exchange of coupling boundary data; for example, see [2-4, 8-12, 14, 15, 17, 19].This is easily described using the concept of matrix splitting [20][21][22].We define a matrix splitting  =  −  so that Starting with some initial guess  0 , the following equation defines a fixed point iteration which may converge to û+1 .
Note  +1,−1 represents the exchange of information between the components, and the subsequent inversion of  is the solution of the individual components.In the case of (3), the splitting  =  −  is chosen to separate the solution of the problems on the subdomains [20][21][22].This is motivated by the fact that the coupling of the subdomains occurs on a manifold of lower dimension.This splitting can be accomplished in a "Gauss-Jacobi" sense: or in a "Gauss-Seidel" sense: If this iteration converges, it converges to the implicitly coupled solution û, which is the unique fixed point of the iteration.By differencing ( 5) and ( 6), we get where  +1, =  +1, −û +1 .This leads to the well-known result [20,21]: The iteration converges if the spectral radius of  −1  is less than one, and it diverges if the spectral radius is greater than one [20][21][22].If the iteration converges, then we can get as close as desired to the implicitly coupled solutions û by iterating.If for any reason we are unable to iterate until convergence at each time step, the situation is considerably more complicated.If û cannot be obtained at a given time step, then it cannot be used as the initial condition for the next time step, and the iterates   may wander away from the implicitly coupled solutions û .To investigate how the iterates wander from their implicitly coupled counterparts, we must account for both the error from incomplete iteration at every time step and for the error associated with the inherited initial value at every time step.

Stability in Time for a Single Iteration.
In this section, we derive an error formula for the case of a single iteration at each time step.The motivation is practical both as this is often encountered in practice and because the relatively simple result shows how the general result proceeds.The error is defined to be the difference between the computed solution and the implicitly coupled solution at each time step.Recall the implicitly coupled solution satisfies For the case of a single iteration per time step, we have  = 1 and we can suppress the iteration index.The iteration is simply If we define   =   − û and Δû +1 = û − û+1 , then (12) becomes Using (11), we obtain Substituting  = 0 and  = 1 into ( 14) and rearranging gives Substituting ( 15) into (16) gives In general, the error at time step  is While the form of ( 18) is not as concise as (10), it is clear that the stability in time depends on the spectral radius of

Stability in Time for 𝑛 Iterations.
In this section, we derive a more general version of (18) for the case of  iterations at each time step.Beginning with (11), the iterates satisfy The error is now defined as  , =  , − û for  = 1, . . ., .We are interested in obtaining a formula for  , .
The following relationships are useful in deriving the final result:  +1,0 =  , = û +  , = û+1 + Δû +1 +  , ,  +1, = û+1 +  +1, ,  = 1, . . ., . (20) Substituting these into (19) gives the following, for  = 1: and, for  = 2, . . ., , Substituting into (11) gives the following, for  = 1: and, for  = 2, . . ., , Starting with  0 =  0 − û0 , we can use ( 23) and ( 24) to find  1, and  2, and discover the general form for  , .The process is tedious, but the pattern is quickly apparent.The final form for the error at time step  where  iterations have been taken per time step is The form of   is as follows.For  = 1, and, for  > 1, the following recursion relationship holds: For example, Note that the results derived above assume that  is being inverted exactly.A modification of (25) for the case of inexact inversion is given in Appendix B. A second modification of (25), incorporating the use of a weighted Jacobi method in which the new iterate is taken to be a weighted average of the old iterate and the Jacobi iterate, is given in Appendix C. We conclude of course that the stability in time is determined by the spectral radius of the matrix   .However, given the conventional wisdom, it is important to point out that stability does not imply accuracy.If we examine (25), it is clear from the term [  ]   0 that if   has a spectral radius of greater than one then even tiny machine rounding errors will grow out of control [20], and this is one way to describe the notion of instability.If the spectral radius of   is less than one, this can not occur and the method can then be considered stable.However, the summation term Δû  shows that it is still possible for error to accumulate in time, and the size of this error depends on the size of the vectors Δû as well as the effect of  −1  and   on Δû.We consider these errors to be an issue of accuracy, not stability.This calls into question the conventional wisdom that if the result does not blow up, then the method is stable and therefore accurate.
We emphasize that the error formulas derived above measure the difference between the computed numerical solution and the idealized numerical solution obtained by solving the implicitly coupled discrete equations.The idealized numerical solution is itself a discrete approximation to the continuous solution of the PDE.The accuracy of the implicitly coupled solution relative to the continuous solution is a separate matter and would have some order in ℎ and , where ℎ and  represent the cell width and time step, such as (ℎ 2 ) + () assuming a typical second order accurate spatial discretization and first order accurate time discretization is used.Preserving this ideal accuracy requires that discrete effects such as solution error and errors arising from projecting between component discretizations are minimal [19].

Relationship between 𝑇 𝑛 and 𝐴
This form suggests a power series representation for an approximate inverse.Rewriting the splitting  =  −  as Using a power series for the inverse in brackets [20][21][22], assuming that ‖ −1 ‖ < 1, we get Hence, provided ‖ −1 ‖ < 1,   becomes an increasingly accurate approximation of  −1 as the number of iterations  increases.Specifically, the first term in (29) approaches (31) as  increases, while the second term in (29) goes to zero.

Numerical Examples
We can employ the error formulas derived above by forming the matrix   for any given set of discretization parameters, for example, step size, space mesh, and number of iterations per time step, and then easily examine the stability of the overall algorithm for particular choices.In this section, we present numerical experiments designed to examine aspects of the conventional wisdom.
We consider a one dimensional domain with  ∈ [0, 1] and the interface between the two subdomains located at  = 1/2, and we pose the heat equation with constant diffusivity equal to one in each subdomain ( 1 =  2 = 1).The outer boundary conditions are homogeneous Neumann at  = 0 and homogeneous Dirichlet at  = 1.Error measurement is facilitated by constructing the problems so that the exact continuous solution is U(, ) = (1 −  2 ) − .The complete statement of the PDE is We use finite volume in space and backward Euler in time [13,18,19,21,22].The details of the discrete equations are given in Appendix A. The cell-centered finite-volume method provides state values at cell centers.Fluxes at cell boundaries are calculated via differencing, while linear extrapolation is used to compute coupling values at the boundary between subdomains.Component 1 receives a Dirichlet condition at the interface that is a linear extrapolation of the two state values that are the nearest to the interface in component 2 and, similarly, component 2 receives a Neumann condition at the interface which is a linear extrapolation of the two flux values that are the nearest to the interface in component 1 [8,9].The extrapolation scheme is depicted in Figure 2.
In the examples, we use a Jacobi iteration with  = 1.Once the exchange of information is complete, each subdomain is solved to within machine precision by direct inversion of the matrix .At each time step, the code produces both the implicitly coupled solution, û, by inverting the coupled system matrix  and the iterative solution, , by performing the specified number of iterations.The implicitly coupled and iterative solutions can then be compared at any time step.It is easily verified that the implicitly coupled solution is unconditionally stable and has accuracy (ℎ 2 ) + ().In the following examples, we concentrate on the error between the implicitly coupled and iterative solutions.We introduce the relative error, which is based on the standard discrete two-norm:

Advances in Mathematical Physics
Finally, we conduct similar experiments with a wide range of spatial mesh sizes.The qualitative results are the same in every case, although the specific threshold values may of course vary.

Case 1: Stability of the Algorithm with 𝑛 Iterations Using a Fixed Time
Step, .In the first experiment, we set the grid size to be 16 cells in each subdomain (ℎ = 1/32) and fix the time step at  = 1/40.We plot the spectral radius of   ,  −1 , and  −1 for various values of  in Figure 3.Note that the spectral radius of  −1  is less and one, so the iteration is convergent, and the spectral radius of  −1 is less than one, so the implicitly coupled scheme would be unconditionally stable if we had the luxury of inverting  −1 at each time step.However, the method is unstable for certain values of .If one observes the progress of the iterative solution at a given time step, it is clear that although the iteration is convergent, certain early iterates contain significantly more error than the previous iterate.In other words, the error in the computational iterate does not decrease monotonically.For the particular problem examined here, it is every fourth iterate, starting with the second iterate, that has the increased error.Figure 3 shows that for the first few  values in this sequence ( = 2, 6, 10, 14) the reduced quality of the iterate results in the entire algorithm being unstable, despite the fact that the iteration is convergent.As the number of iterations increases, this effect is reduced in magnitude and eventually disappears.We verify the instability by solving to  = 2 (80 time steps) and listing the relative error at  = 2 for several values of  in Table 1.
The values  = 200 and  = 400 are included in Table 1 to verify that, since the iteration is convergent, the relative error approaches zero as the number of iterations becomes very large.The reduction in relative error for large  is fairly slow, since the spectral radius of  −1  is approximately .975, as indicated in Figure 3.
The results of this experiment serve to illustrate one of our main points that a low number of iterations can interfere with the unconditional stability of the time discretization.In addition, it shows that a lower number of iterations can be stable, while a higher number is unstable.The plot of the spectral radius of   in Figure 3 indicates that using 2, 6, 10, or 14 iterations leads to an unstable method.The numerical values in Table 1 confirm the instability for  = 2 and  = 6 ( = 10 and  = 14 yield analogous results).Backward Euler is considered "unconditionally stable" [21,22]; however, this fact is based on the assumption that the system is solved exactly at each time step.For the iterative algorithm used here, this unconditional stability is only valid in the limit as the number of iterations goes to infinity.This example shows how sensitive the stability of the algorithm is to the number of iterations used.Note that this example also provides another demonstration of the fact that stability does not imply accuracy.Notice in Table 1 that the relative error for the stable case  = 4 is around ten times that of the nearby stable cases,  = 1, 3, 5. Despite the fact that the spectral radius of   is less than one for all the stable cases, it is still possible for the summation term in (25) to accumulate significant error in time.This term makes a much larger contribution for the case of  = 4 than it does for the other stable cases.This is not evident from the spectral radius of   alone; note that ( 1 ) > ( 4 ).

Case 2: Stability of the Algorithm with Time
Step , Using a Fixed Number of Iterations, .In the second experiment, we use identical conditions except that we fix the number of iterations  and vary the time step .Figure 4 shows the spectral radii of the relevant matrices for  = 1 and values of  ranging from 0 to 0. of time step  corresponds to a range for the ratio /ℎ 2 of approximately 1 to 100. Figure 4 implies that the method with  = 1 is stable for all  < 0.1.Carrying out the calculation confirms this, and the relative errors at the end of the simulation are given in Table 2.
Figure 5 shows the spectral radii of the relevant matrices for  = 2.The implication is that the method with  = 2 is stable only for very small values of .Table 3 gives relative errors for the  = 2 case.
For the case of  = 6, the plot and table are shown in Figure 6 and Table 4.This case illustrates that instability can occur for a limited range of time steps with a minimum and a maximum, not simply for time steps above a certain threshold, as one might expect.The  = 6 case makes it clear that the method is stable for sufficiently small and sufficiently large time steps and is unstable in between.(In fact this is also true for the  = 2 case, and if the horizontal axis in Figure 5 is carried out to much larger time steps the spectral radius of  2 does eventually drop below one.)This is a surprising result since one would not expect reducing the time step to promote instability, nor increasing the time step to restore stability.Finally, in Figure 7 we provide two further plots for the  = 10 and  = 18 cases.These show the relationship between   and  −1 discussed in Section 3.3.The spectral radius of  −1 is less than one for all values of , and it decreases monotonically with increasing .We know from the expansions in Section 3.3 that as  becomes large then   approaches  −1 .However, there is no reason to expect the spectral radius of   to be monotonic in  for small values of , and the plots in this section show that it is not.There is a characteristic bump in the plots of the spectral radius of   versus , and it is in this range of  values that the spectral radius of   may exceed 1, leading to an unstable method.Naturally, as  becomes large, the size of this bump is reduced, since   is driven closer to  −1 .In the  = 18 case, the spectral radius of   is less than 1 for all time steps.

The Multirate Case
Multirate integration using different time steps for the different components is a common practice [1,8,9].In this case, the implicitly coupled system and error formulas are changed slightly.We assume the longer time step is an integer multiple of the shorter time step, so component 2 will take  steps for every one step in the component 1.The outline of the derivation for  = 2 is given here, but results easily generalize.
The large time steps correspond to the index  and the small time steps correspond to fractional indices.First, component 1 equation is straightforward since it only involves one time step: Next, component 2 equation requires compounding two steps: After carrying out the algebra and generalizing to larger , we obtain the system: The fractional index is used to indicate the small time steps.Now that the implicitly coupled system has been defined, we can choose a splitting that defines  and  in the error formulas above.The error formula for the multirate case is very similar to (25), but we must define matrices  , where  is the number of iterations per time step and  is the integer number of time steps taken in component 2 for each time step taken in component 1.Let where  1 is an identity matrix of the same size as  1 .The recursion relationship for the  , is The alternative form for  , , analogous to (29) is This means that as  goes to infinity,  , approaches  −1   (note that  1 = ).Finally, the multirate error formula is Figure 9: Spectral radii for multirate example 2.  = 3,  = 2 (a), and  = 100,  = 2 (b).Black is  −1 , dotted black is  −1   , green is  −1 , blue is   , and red marks the value one.The horizontal axis shows the size of the large time step.
In the first multirate example,  is set to 1, and  is set first to 2 and then to 10. Figure 8 shows the spectral radii of the relevant matrices.
Values of  larger than 10 produce no visible change in the plot, so Figure 8(b) can be taken to represent the case of "many" time steps in component 2 within each time step for component 1.There is now a range of time steps for which the iteration diverges.As might be expected, the algorithm is unstable when the iteration is divergent.In other words, when the spectral radius of  −1  is greater than one, the spectral radius of   is greater than one.
As the number of iterations, , approaches infinity, a divergent iteration must result in an unstable algorithm.However, for the case of finite iteration, this is not guaranteed, and the next experiment illustrates this point.

Multirate Example 2:
Varying  with Constant .In the second multirate example,  is set to 2, and  is set first to 3 and then to 100. Figure 9 shows the spectral radii for the relevant matrices.The results show that the method may be stable, even though the iteration is divergent.Figure 9(a) include a region in which the spectral radius of  −1  is more than one, yet the spectral radius of   is less than one.In Figure 9(a), where only 3 iterations are used per large time step, the entire range of  values for which the iteration is divergent results in a stable algorithm.In Figure 9(b), where the number of iterations is raised to 100, the region of divergence and the region of instability are very nearly the same.time step is small, since iterating many times with a divergent iteration would certainly lead to an unstable algorithm.This is an interesting case, since a stable algorithm could contain a divergent iteration and the user might not know.

Multirate Example in
The purpose of this example is to examine the changes in the spectral radii of the relevant matrices resulting from the change to two spatial dimensions.We use the standard finite volume discretization with 16 × 16 cells in each subdomain.Since the functions   ,   ,   , and   have no impact on the relevant matrices, we do not state them explicitly.Plots of the spectral radii of the relevant matrices are given in Figures 10 and 11.The spectral radius of  −1 is much smaller in the two space dimension case.Nonetheless, the plots show that the two space dimension case exhibits the same qualitative behavior as the one space dimension case.In particular, there is a range of time steps for which the iteration is divergent and also a range of time steps for which the overall method is unstable.

Conclusion
By selecting a coupling strategy comprising space and time grids and an associated rule by which information is exchanged between components, we define an implicitly coupled system.At each time step, we seek a solution of this implicitly coupled system through a block iterative strategy.Since only a limited number of iterations can be performed at each time step, there will be some iteration error which separates the final iterate from the implicitly coupled solution.These errors are propagated to the next time step in the form of errors in the initial condition and are compounded as the incomplete iteration process repeats itself at each time step.The cumulative effect can manifest itself as conditional stability, meaning the solution is only stable for certain values of  despite the unconditional stability of the implicitly coupled solution.
We have derived formulas for this error which show that stability hinges on the spectral radius of the matrix   and further show that there is a mechanism for error to accumulate in time even if the algorithm is stable.By examining the spectral radius of   for a variety of discretization choices on a simple model problem we demonstrate several important points.Firstly, the unconditional stability of the time integration method is not necessarily retained if a low number of iterations per timestep are used.Instead, a conditional stability may occur, where the method is stable for a range of values of .In additon, we show by example that it is possible for the method to be stable even if the iteration is divergent, provided a low number of iterations is used.These results seem to contradict intuition, which suggests that a convergent iteration will produce a stable algorithm and a divergent iteration will produce an unstable algorithm.However, both of these ideas hold only in the limit as the number of iterations goes to infinity, and in the case of finite iteration we have shown that there are exceptions.
Finally, the question of accuracy must be considered in addition to the issue of stability.The error formulas derived above show that the iterative solution may wander away from the implicitly coupled solution over time even if the method is stable in the traditional sense.where  +1 is the known Dirichlet boundary condition at  =   at time  + 1.Note that (A.1) and (A.2) are valid for both of our subdomains provided that  +1 and  +1 are interpreted as given boundary conditions on the outer boundaries of the entire domain and are determined by the extrapolation strategy depicted in Figure 2 at the interface between the subdomains.

B. Inexact Component Solutions (Inexact Inversion of 𝑀)
The error formulas derived in Sections 3.1 and 3.2 can easily be modified to include the effects of inexact inversion of .We first adjust (18)  Notice that, with multiple iterations per time step, the  term needs a double index because an error due to inexact inversion is made at every iteration.Equation (B.4) contains both terms from (25) and has an additional term which represents the error caused by the inexact inversions.The matrix   is the key to the growth of all three terms and therefore remains the key to stability.

C. Weighted Jacobi
It is worth noting that the error formula (25) can easily be altered to allow for a "relaxed" Jacobi iteration [20][21][22], in which a weighted average of the new and old iterates of a standard Jacobi iteration is taken to be the next iterate (see, e.g., [23]).If we rewrite (25) as With this definition of   , (25) holds for any value of  ∈ [0, 1], with  = 1 corresponding to standard iteration.In some cases where the iterative convergence at each time step is slow, the iterates tend to oscillate around the implicitly coupled solution at that time step.In these cases, adjusting  can drastically accelerate iterative convergence.

Figure 2 :
Figure 2: Diagram of the coupling strategy used in the numerical experiments, in which linear extrapolation from the last two available state or flux values in one subdomain are used to compute a boundary value for the other subdomain.

Figure 3 :
Figure 3: Spectral radius for various values of  with fixed .Black is  −1 , green is  −1 , blue is   , and red marks the value one.

Figure 7 :
Figure 7: Spectral radius for various values of  with  = 10 (a), and  = 18 (b).Black is  −1 , green is  −1 , blue is   , and red marks the value one.

Figure 10 :
Figure 10: Spectral radii for multirate example 1 in 2 spatial dimensions. = 1,  = 2 (a), and  = 1,  = 10 (b).Black is  −1 , dotted black is  −1   , green is  −1 , blue is   , and red marks the value one.The horizontal axis shows the size of the large time step.

Figure 11 :
Figure 11: Spectral radii for multirate example 2 in 2 spatial dimensions. = 3,  = 2 (a), and  = 100,  = 2 (b).Black is  −1 , dotted black is  −1   , green is  −1 , blue is   , and red marks the value one.The horizontal axis shows the size of the large time step.
Let there be  cells in the 1D spatial grid in a given subdomain,  ∈ [  ,   ].Let    be the discrete solution on cell  at time ,  the constant time step, and ℎ the constant cell width.The symbol    represents the exact integral over cell  of the right hand side of the PDE evaluated at time .The symbol   +1/2 represents the diffusivity function evaluated at the boundary between cell  and cell  + 1, at time .The system of equations below, when placed into matrix form, describes the terms in (3) as they were implemented for the numerical tests. +1 is the known Neumann boundary condition at  =   at time  + 1.

Table 1 :
Relative error at  = 2 for various numbers of iterations.

Table 2 :
Figure 4: Spectral radius for various values of  with  = 1.Black is  −1 , green is  −1 , blue is   , and red marks the value one.Figure 5: Spectral radius for various values of  with  = 2. Black is  −1 , green is  −1 , blue is   , and red marks the value one.Relative error at  = 2 for  = 1 with various values of .

Table 3 :
Relative error at  = 2 for  = 2 with various values of .Note that for larger values of , fewer time steps are needed to reach  = 2, so the relative error does not grow as large.Figure6: Spectral radius for various values of  with  = 6.Black is  −1 , green is  −1 , blue is   , and red marks the value one.

Table 4 :
Relative error at  = 2 for  = 6 with various values of .
Such cases are not unique to multirate examples and have also been observed in single rate examples.Of course these cases can only occur when the number of iterations per