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

Authors: B. Sheehan; D. Estep; S. Tavener; J. Cary; S. Kruger; A. Hakim; A. Pletzer; J. Carlsson; S. Vadlamani
Journal: Advances in Mathematical Physics
 (2015)

---

## Abstract

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.


---

## Body

## 1. 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–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–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–4, 14–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–4, 19]. However, in practice, as mentioned, the component physics is usually complex and there are large differences in the discretization strategies employed in 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. Section2 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.


## 2. 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.Figure 1
Illustration of neighboring subdomains.The general form of the PDE system is(1)∂∂tU1+L1U1=S1on⁡Ω1×0,T,∂∂tU2+L2U2=S2on⁡Ω2×0,T,U1=U2on⁡Γ×0,T,n^·a1∇U1=n^·a2∇U2on⁡Γ×0,T,
together with the following additional boundary and initial conditions: (i)
boundary conditions forU1 on ∂Ω1∖Γ and initial conditions for U1 on Ω1,(ii)
boundary conditions forU2 on ∂Ω2∖Γ and initial conditions for U2 on Ω2,where L1v=∇·(a1∇v) and L2v=∇·(a2∇v) are linear, coercive elliptic operators, S1, S2 are given functions in L2, 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:
(2)L1u^1k+1=u^1k+S1k+1+P21u^2k+1,L2u^2k+1=u^2k+S2k+1+P12u^1k+1,
where the superscript is the time step index. The symbols L1 and L2 are the discrete versions of the operators (∂/∂t+L1) and (∂/∂t+L2), except that u^k has been moved to the right hand side. The symbols S1 and S2 are the discrete versions of the source terms S1 and S2 and also take into account the known boundary data on ∂Ω1∖Γ and ∂Ω2∖Γ. P21 and P12 are the projections by which the solution in one subdomain is used to provide boundary data on Γ for the other subdomain. The subscripts “12” and “21” indicate information passing from “Ω1 to Ω2” and “Ω2 to Ω1.” We refer to u^1 and u^2, which solve (2), as the implicitly coupled solution. We emphasize that we analyze this method because it is one of the most common approaches used in practice.Note that the form of the operatorsP21 and P12 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 P21 and P12 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,
(3)Au^k+1=u^k+Sk+1,
where
(4)A=L1-P21-P12L2,S=S1S2,u^=u^1u^2.


### 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–22]. We define a matrix splitting A=M-N so that
(5)Mu^k+1=u^k+Sk+1+Nu^k+1.
Starting with some initial guess u0, the following equation defines a fixed point iteration which may converge to u^k+1. For the iterates, we use a double superscript, k for the time step and i for the iteration index. We use n for the number of iterations to be performed at each time step:
(6)Muk+1,i=uk,n+Sk+1+Nuk+1,i-1,kkikkii=1,…,n,uk+1,0=uk,n.
Note Nuk+1,i-1 represents the exchange of information between the components, and the subsequent inversion of M is the solution of the individual components. In the case of (3), the splitting A=M-N is chosen to separate the solution of the problems on the subdomains [20–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:
(7)M=L100L2,N=0P21P120,
or in a “Gauss-Seidel” sense:
(8)M=L10-P12L2,N=0P2100.If this iteration converges, it converges to the implicitly coupled solutionu^, which is the unique fixed point of the iteration. By differencing (5) and (6), we get
(9)Mek+1,i+1=Nek+1,i,
where ek+1,i=uk+1,i-u^k+1. This leads to the well-known result [20, 21]:
(10)ek+1,n=M-1Nnek+1,0.
The iteration converges if the spectral radius of M-1N is less than one, and it diverges if the spectral radius is greater than one [20–22]. If the iteration converges, then we can get as close as desired to the implicitly coupled solutions u^k by iterating. If for any reason we are unable to iterate until convergence at each time step, the situation is considerably more complicated. If u^k 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 uk may wander away from the implicitly coupled solutions u^k. 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.


## 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–22]. We define a matrix splitting A=M-N so that
(5)Mu^k+1=u^k+Sk+1+Nu^k+1.
Starting with some initial guess u0, the following equation defines a fixed point iteration which may converge to u^k+1. For the iterates, we use a double superscript, k for the time step and i for the iteration index. We use n for the number of iterations to be performed at each time step:
(6)Muk+1,i=uk,n+Sk+1+Nuk+1,i-1,kkikkii=1,…,n,uk+1,0=uk,n.
Note Nuk+1,i-1 represents the exchange of information between the components, and the subsequent inversion of M is the solution of the individual components. In the case of (3), the splitting A=M-N is chosen to separate the solution of the problems on the subdomains [20–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:
(7)M=L100L2,N=0P21P120,
or in a “Gauss-Seidel” sense:
(8)M=L10-P12L2,N=0P2100.If this iteration converges, it converges to the implicitly coupled solutionu^, which is the unique fixed point of the iteration. By differencing (5) and (6), we get
(9)Mek+1,i+1=Nek+1,i,
where ek+1,i=uk+1,i-u^k+1. This leads to the well-known result [20, 21]:
(10)ek+1,n=M-1Nnek+1,0.
The iteration converges if the spectral radius of M-1N is less than one, and it diverges if the spectral radius is greater than one [20–22]. If the iteration converges, then we can get as close as desired to the implicitly coupled solutions u^k by iterating. If for any reason we are unable to iterate until convergence at each time step, the situation is considerably more complicated. If u^k 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 uk may wander away from the implicitly coupled solutions u^k. 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.


## 3. Analysis of Stability of the Iterative Solution

### 3.1. 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(11)Mu^k+1=u^k+Sk+1+Nu^k+1.
For the case of a single iteration per time step, we have n=1 and we can suppress the iteration index. The iteration is simply
(12)Muk+1=uk+Sk+1+Nuk.
If we define ek=uk-u^k and Δu^k+1=u^k-u^k+1, then (12) becomes
(13)Mu^k+1+ek+1=u^k+ek+Sk+1+Nu^k+1+Δu^k+1+ek.
Using (11), we obtain
(14)Mek+1=ek+NΔu^k+1+ek.Substitutingk=0 and k=1 into (14) and rearranging gives
(15)e1=M-1+M-1Ne0+M-1NΔu^1,(16)e2=M-1+M-1Ne1+M-1NΔu^2.
Substituting (15) into (16) gives
(17)e2=M-1+M-1N2e0+M-1+M-1N×M-1NΔu^1+M-1NΔu^2.In general, the error at time stepk is
(18)ek=M-1+M-1Nke0+∑j=1k‍M-1+M-1N(k-j)M-1NΔu^j.
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 [M-1+M-1N].


### 3.2. Stability in Time forn Iterations

In this section, we derive a more general version of (18) for the case of n iterations at each time step. Beginning with (11), the iterates satisfy
(19)Muk+1,i=uk,n+Sk+1+Nuk+1,i-1,kkkkkii=1,…,n,uk+1,0=uk,n.
The error is now defined as ek,i=uk,i-u^k for i=1,…,n. We are interested in obtaining a formula for ek,n.The following relationships are useful in deriving the final result:(20)uk+1,0=uk,n=u^k+ek,n=u^k+1+Δu^k+1+ek,n,uk+1,i=u^k+1+ek+1,i,i=1,…,n.
Substituting these into (19) gives the following, for i=1:
(21)Mu^k+1+ek+1,1=u^k+ek,n+Sk+1+Nu^k+1+Δu^k+1+ek,n,
and, for i=2,…,n,
(22)Mu^k+1+ek+1,i=u^k+ek,n+Sk+1+Nu^k+1+ek+1,i-1.
Substituting into (11) gives the following, for i=1:
(23)Mek+1,1=ek,n+NΔu^k+1+ek,n,
and, for i=2,…,n,
(24)Mek+1,i=ek,n+Nek+1,i-1.Starting withe0=u0-u^0, we can use (23) and (24) to find e1,n and e2,n and discover the general form for ek,n. The process is tedious, but the pattern is quickly apparent. The final form for the error at time step k where n iterations have been taken per time step is
(25)ek,n=Tnke0+∑j=1k‍Tnk-jM-1NnΔu^j.
The form of Tn is as follows. For n=1,
(26)T1=M-1+M-1N,
and, for n>1, the following recursion relationship holds:
(27)Tn=M-1+M-1NTn-1.
For example,
(28)T2=M-1+M-1NM-1+M-1N,T3=M-1+M-1NM-1+M-1NM-1+M-1N.
Note that the results derived above assume that M 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 matrixTn. 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 [Tn]ke0 that if Tn 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 Tn is less than one, this can not occur and the method can then be considered stable. However, the summation term ∑j=1k[Tn]k-j(M-1N)nΔu^j 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 Δu^ as well as the effect of M-1N and Tn on Δu^. 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 inh and dt, where h and dt represent the cell width and time step, such as O(h2)+O(dt) 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].


### 3.3. Relationship betweenTn and A-1

An alternative form for theTn is
(29)Tn=∑i=0n-1‍M-1NiM-1+M-1Nn.
This form suggests a power series representation for an approximate inverse. Rewriting the splitting A=M-N as (M-N)=M[I-(M-1N)], then
(30)M-N-1=I-M-1N-1M-1.
Using a power series for the inverse in brackets [20–22], assuming that M-1N<1, we get
(31)M-N-1=∑i=0∞‍M-1NiM-1.
Hence, provided M-1N<1, Tn becomes an increasingly accurate approximation of A-1 as the number of iterations n increases. Specifically, the first term in (29) approaches (31) as n increases, while the second term in (29) goes to zero.


## 3.1. 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(11)Mu^k+1=u^k+Sk+1+Nu^k+1.
For the case of a single iteration per time step, we have n=1 and we can suppress the iteration index. The iteration is simply
(12)Muk+1=uk+Sk+1+Nuk.
If we define ek=uk-u^k and Δu^k+1=u^k-u^k+1, then (12) becomes
(13)Mu^k+1+ek+1=u^k+ek+Sk+1+Nu^k+1+Δu^k+1+ek.
Using (11), we obtain
(14)Mek+1=ek+NΔu^k+1+ek.Substitutingk=0 and k=1 into (14) and rearranging gives
(15)e1=M-1+M-1Ne0+M-1NΔu^1,(16)e2=M-1+M-1Ne1+M-1NΔu^2.
Substituting (15) into (16) gives
(17)e2=M-1+M-1N2e0+M-1+M-1N×M-1NΔu^1+M-1NΔu^2.In general, the error at time stepk is
(18)ek=M-1+M-1Nke0+∑j=1k‍M-1+M-1N(k-j)M-1NΔu^j.
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 [M-1+M-1N].


## 3.2. Stability in Time forn Iterations

In this section, we derive a more general version of (18) for the case of n iterations at each time step. Beginning with (11), the iterates satisfy
(19)Muk+1,i=uk,n+Sk+1+Nuk+1,i-1,kkkkkii=1,…,n,uk+1,0=uk,n.
The error is now defined as ek,i=uk,i-u^k for i=1,…,n. We are interested in obtaining a formula for ek,n.The following relationships are useful in deriving the final result:(20)uk+1,0=uk,n=u^k+ek,n=u^k+1+Δu^k+1+ek,n,uk+1,i=u^k+1+ek+1,i,i=1,…,n.
Substituting these into (19) gives the following, for i=1:
(21)Mu^k+1+ek+1,1=u^k+ek,n+Sk+1+Nu^k+1+Δu^k+1+ek,n,
and, for i=2,…,n,
(22)Mu^k+1+ek+1,i=u^k+ek,n+Sk+1+Nu^k+1+ek+1,i-1.
Substituting into (11) gives the following, for i=1:
(23)Mek+1,1=ek,n+NΔu^k+1+ek,n,
and, for i=2,…,n,
(24)Mek+1,i=ek,n+Nek+1,i-1.Starting withe0=u0-u^0, we can use (23) and (24) to find e1,n and e2,n and discover the general form for ek,n. The process is tedious, but the pattern is quickly apparent. The final form for the error at time step k where n iterations have been taken per time step is
(25)ek,n=Tnke0+∑j=1k‍Tnk-jM-1NnΔu^j.
The form of Tn is as follows. For n=1,
(26)T1=M-1+M-1N,
and, for n>1, the following recursion relationship holds:
(27)Tn=M-1+M-1NTn-1.
For example,
(28)T2=M-1+M-1NM-1+M-1N,T3=M-1+M-1NM-1+M-1NM-1+M-1N.
Note that the results derived above assume that M 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 matrixTn. 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 [Tn]ke0 that if Tn 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 Tn is less than one, this can not occur and the method can then be considered stable. However, the summation term ∑j=1k[Tn]k-j(M-1N)nΔu^j 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 Δu^ as well as the effect of M-1N and Tn on Δu^. 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 inh and dt, where h and dt represent the cell width and time step, such as O(h2)+O(dt) 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].


## 3.3. Relationship betweenTn and A-1

An alternative form for theTn is
(29)Tn=∑i=0n-1‍M-1NiM-1+M-1Nn.
This form suggests a power series representation for an approximate inverse. Rewriting the splitting A=M-N as (M-N)=M[I-(M-1N)], then
(30)M-N-1=I-M-1N-1M-1.
Using a power series for the inverse in brackets [20–22], assuming that M-1N<1, we get
(31)M-N-1=∑i=0∞‍M-1NiM-1.
Hence, provided M-1N<1, Tn becomes an increasingly accurate approximation of A-1 as the number of iterations n increases. Specifically, the first term in (29) approaches (31) as n increases, while the second term in (29) goes to zero.


## 4. Numerical Examples

We can employ the error formulas derived above by forming the matrixTn 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 withx∈[0,1] and the interface between the two subdomains located at x=1/2, and we pose the heat equation with constant diffusivity equal to one in each subdomain (a1=a2=1). The outer boundary conditions are homogeneous Neumann at x=0 and homogeneous Dirichlet at x=1. Error measurement is facilitated by constructing the problems so that the exact continuous solution is U(x,t)=(1-x2)e-t. The complete statement of the PDE is
(32)∂∂tU1-∂2∂x2U1=1+x2e-t,x∈0,.5,t∈0,T,∂∂tU2-∂2∂x2U2=1+x2e-t,x∈.5,1,t∈0,T,U1=U2,x=.5,t∈0,T,∂∂xU1=∂∂xU2,x=.5,t∈0,T,∂∂xU10,t=0,U21,t=0,U1x,0=U2x,0=1-x2.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.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.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 M. At each time step, the code produces both the implicitly coupled solution, u^, by inverting the coupled system matrix A and the iterative solution, u, 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 O(h2)+O(dt). 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:
(33)erel⁡=u-u^2u^2.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.


### 4.1. Case 1: Stability of the Algorithm withn Iterations Using a Fixed Time Step, dt

In the first experiment, we set the grid size to be 16 cells in each subdomain (h=1/32) and fix the time step at dt=1/40. We plot the spectral radius of Tn, M-1N, and A-1 for various values of n in Figure 3. Note that the spectral radius of M-1N is less and one, so the iteration is convergent, and the spectral radius of A-1 is less than one, so the implicitly coupled scheme would be unconditionally stable if we had the luxury of inverting A-1 at each time step. However, the method is unstable for certain values of n. 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 n values in this sequence (n=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 t=2 (80 time steps) and listing the relative error at t=2 for several values of n in Table 1.Table 1
Relative error att=2 for various numbers of iterations.
n
e
rel
⁡
1
.
05
2
4.4
*
1
0
9
3
.
09
4
.
65
5
.
05
6
8.5
*
1
0
5
200
1.8
*
1
0
-
4
400
5.6
*
1
0
-
7Figure 3
Spectral radius for various values ofn with fixed dt. Black is A-1, green is M-1N, blue is Tn, and red marks the value one.The valuesn=200 and n=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 n is fairly slow, since the spectral radius of M-1N 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 ofTn 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 n=2 and n=6 (n=10 and n=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 n=4 is around ten times that of the nearby stable cases, n=1,3,5. Despite the fact that the spectral radius of Tn 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 n=4 than it does for the other stable cases. This is not evident from the spectral radius of Tn alone; note that ρ(T1)>ρ(T4).


### 4.2. Case 2: Stability of the Algorithm with Time Stepdt, Using a Fixed Number of Iterations, n

In the second experiment, we use identical conditions except that we fix the number of iterationsn and vary the time step dt. Figure 4 shows the spectral radii of the relevant matrices for n=1 and values of dt ranging from 0 to 0.1. Since h=1/32, this range of time step dt corresponds to a range for the ratio dt/h2 of approximately 1 to 100. Figure 4 implies that the method with n=1 is stable for all dt<0.1. Carrying out the calculation confirms this, and the relative errors at the end of the simulation are given in Table 2.Table 2
Relative error att=2 for n=1 with various values of dt.
d
t
e
rel
⁡
.
005
.
01
.
05
.
08
.
1
.
22Figure 4
Spectral radius for various values ofdt with n=1. Black is A-1, green is M-1N, blue is Tn, and red marks the value one.Figure5 shows the spectral radii of the relevant matrices for n=2. The implication is that the method with n=2 is stable only for very small values of dt. Table 3 gives relative errors for the n=2 case.Table 3
Relative error att=2 for n=2 with various values of dt. Note that for larger values of dt, fewer time steps are needed to reach t=2, so the relative error does not grow as large.
d
t
e
rel
⁡
.
0025
.
01
.
005
5.7
*
1
0
16
.
05
2.5
*
1
0
4
.
1
19Figure 5
Spectral radius for various values ofdt with n=2. Black is A-1, green is M-1N, blue is Tn, and red marks the value one.For the case ofn=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 n=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 n=2 case, and if the horizontal axis in Figure 5 is carried out to much larger time steps the spectral radius of T2 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.Table 4
Relative error att=2 for n=6 with various values of dt.
d
t
e
rel
⁡
.
005
.
0135
.
05
117
.
1
.
0114Figure 6
Spectral radius for various values ofdt with n=6. Black is A-1, green is M-1N, blue is Tn, and red marks the value one.Finally, in Figure7 we provide two further plots for the n=10 and n=18 cases. These show the relationship between Tn and A-1 discussed in Section 3.3. The spectral radius of A-1 is less than one for all values of dt, and it decreases monotonically with increasing dt. We know from the expansions in Section 3.3 that as n becomes large then Tn approaches A-1. However, there is no reason to expect the spectral radius of Tn to be monotonic in dt for small values of n, and the plots in this section show that it is not. There is a characteristic bump in the plots of the spectral radius of Tn versus dt, and it is in this range of dt values that the spectral radius of Tn may exceed 1, leading to an unstable method. Naturally, as n becomes large, the size of this bump is reduced, since Tn is driven closer to A-1. In the n=18 case, the spectral radius of Tn is less than 1 for all time steps.Spectral radius for various values ofdt with n=10 (a), and n=18 (b). Black is A-1, green is M-1N, blue is Tn, and red marks the value one.
(a)
(b)


## 4.1. Case 1: Stability of the Algorithm withn Iterations Using a Fixed Time Step, dt

In the first experiment, we set the grid size to be 16 cells in each subdomain (h=1/32) and fix the time step at dt=1/40. We plot the spectral radius of Tn, M-1N, and A-1 for various values of n in Figure 3. Note that the spectral radius of M-1N is less and one, so the iteration is convergent, and the spectral radius of A-1 is less than one, so the implicitly coupled scheme would be unconditionally stable if we had the luxury of inverting A-1 at each time step. However, the method is unstable for certain values of n. 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 n values in this sequence (n=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 t=2 (80 time steps) and listing the relative error at t=2 for several values of n in Table 1.Table 1
Relative error att=2 for various numbers of iterations.
n
e
rel
⁡
1
.
05
2
4.4
*
1
0
9
3
.
09
4
.
65
5
.
05
6
8.5
*
1
0
5
200
1.8
*
1
0
-
4
400
5.6
*
1
0
-
7Figure 3
Spectral radius for various values ofn with fixed dt. Black is A-1, green is M-1N, blue is Tn, and red marks the value one.The valuesn=200 and n=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 n is fairly slow, since the spectral radius of M-1N 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 ofTn 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 n=2 and n=6 (n=10 and n=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 n=4 is around ten times that of the nearby stable cases, n=1,3,5. Despite the fact that the spectral radius of Tn 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 n=4 than it does for the other stable cases. This is not evident from the spectral radius of Tn alone; note that ρ(T1)>ρ(T4).


## 4.2. Case 2: Stability of the Algorithm with Time Stepdt, Using a Fixed Number of Iterations, n

In the second experiment, we use identical conditions except that we fix the number of iterationsn and vary the time step dt. Figure 4 shows the spectral radii of the relevant matrices for n=1 and values of dt ranging from 0 to 0.1. Since h=1/32, this range of time step dt corresponds to a range for the ratio dt/h2 of approximately 1 to 100. Figure 4 implies that the method with n=1 is stable for all dt<0.1. Carrying out the calculation confirms this, and the relative errors at the end of the simulation are given in Table 2.Table 2
Relative error att=2 for n=1 with various values of dt.
d
t
e
rel
⁡
.
005
.
01
.
05
.
08
.
1
.
22Figure 4
Spectral radius for various values ofdt with n=1. Black is A-1, green is M-1N, blue is Tn, and red marks the value one.Figure5 shows the spectral radii of the relevant matrices for n=2. The implication is that the method with n=2 is stable only for very small values of dt. Table 3 gives relative errors for the n=2 case.Table 3
Relative error att=2 for n=2 with various values of dt. Note that for larger values of dt, fewer time steps are needed to reach t=2, so the relative error does not grow as large.
d
t
e
rel
⁡
.
0025
.
01
.
005
5.7
*
1
0
16
.
05
2.5
*
1
0
4
.
1
19Figure 5
Spectral radius for various values ofdt with n=2. Black is A-1, green is M-1N, blue is Tn, and red marks the value one.For the case ofn=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 n=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 n=2 case, and if the horizontal axis in Figure 5 is carried out to much larger time steps the spectral radius of T2 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.Table 4
Relative error att=2 for n=6 with various values of dt.
d
t
e
rel
⁡
.
005
.
0135
.
05
117
.
1
.
0114Figure 6
Spectral radius for various values ofdt with n=6. Black is A-1, green is M-1N, blue is Tn, and red marks the value one.Finally, in Figure7 we provide two further plots for the n=10 and n=18 cases. These show the relationship between Tn and A-1 discussed in Section 3.3. The spectral radius of A-1 is less than one for all values of dt, and it decreases monotonically with increasing dt. We know from the expansions in Section 3.3 that as n becomes large then Tn approaches A-1. However, there is no reason to expect the spectral radius of Tn to be monotonic in dt for small values of n, and the plots in this section show that it is not. There is a characteristic bump in the plots of the spectral radius of Tn versus dt, and it is in this range of dt values that the spectral radius of Tn may exceed 1, leading to an unstable method. Naturally, as n becomes large, the size of this bump is reduced, since Tn is driven closer to A-1. In the n=18 case, the spectral radius of Tn is less than 1 for all time steps.Spectral radius for various values ofdt with n=10 (a), and n=18 (b). Black is A-1, green is M-1N, blue is Tn, and red marks the value one.
(a)
(b)


## 5. 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 p steps for every one step in the component 1. The outline of the derivation for p=2 is given here, but results easily generalize. The large time steps correspond to the index k and the small time steps correspond to fractional indices. First, component 1 equation is straightforward since it only involves one time step:
(34)L1u^1k+1=u^1k+S1k+1+P21u^2k+1.
Next, component 2 equation requires compounding two steps:
(35)L2u^2k+1/2=u^2k+S2k+1/2+P12u^1k+1,L2u^2k+1=u^2k+1/2+S2k+1+P12u^1k+1.After carrying out the algebra and generalizing to largerp, we obtain the system:
(36)L1-P21-∑j=0p-1‍L2-1jP12L2u^1k+1u^2k+1=u^1k+S1k+1L2-1p-1u^2k+∑j=0p-1‍L2-1jS2k+p-j/p.
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 definesM and N in the error formulas above. The error formula for the multirate case is very similar to (25), but we must define matrices Tn,p where n is the number of iterations per time step and p is the integer number of time steps taken in component 2 for each time step taken in component 1. Let
(37)Zp=I100L2-1p-1,
where I1 is an identity matrix of the same size as L1. The recursion relationship for the Tn,p is
(38)T1,p=M-1Zp+M-1N,Tn,p=M-1Zp+M-1NTn-1,p.
The alternative form for Tn,p, analogous to (29) is
(39)Tn,p=∑i=0n-1‍M-1NiM-1Zp+M-1Nn.
This means that as n goes to infinity, Tn,p approaches A-1Zp (note that Z1=I). Finally, the multirate error formula is
(40)epk,n=Tn,pke0+∑j=1k‍Tn,pk-jM-1NnΔu^j.


### 5.1. Multirate Example 1: Varyingp with Constant n

The multirate examples are the same as those described in Section4.2 with only the values of n and p altered. We use the same Jacobi style iteration, so based on (36) we have
(41)M=L100L2,N=0P21∑j=0p-1‍L2-1jP120.
In the first multirate example, n is set to 1, and p is set first to 2 and then to 10. Figure 8 shows the spectral radii of the relevant matrices.Spectral radii for multirate example 1.n=1, p=2 (a), and n=1, p=10 (b). Black is A-1, dotted black is A-1Zp, green is M-1N, blue is Tn, and red marks the value one. The horizontal axis shows the size of the large time step.
(a)
(b)Values ofp 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 M-1N is greater than one, the spectral radius of Tn is greater than one.As the number of iterations,n, 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.


### 5.2. Multirate Example 2: Varyingn with Constant p

In the second multirate example,p is set to 2, and n 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 M-1N is more than one, yet the spectral radius of Tn is less than one. In Figure 9(a), where only 3 iterations are used per large time step, the entire range of dt 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. 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 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.Spectral radii for multirate example 2.n=3, p=2 (a), and n=100, p=2 (b). Black is A-1, dotted black is A-1Zp, green is M-1N, blue is Tn, and red marks the value one. The horizontal axis shows the size of the large time step.
(a)
(b)


### 5.3. Multirate Example in Two Spatial Dimensions

In this section, the multirate tests in Sections5.1 and 5.2 are extended to a problem in two spatial dimensions posed on x1∈[0,1], x2∈[0,.5] with the interface located at x1=.5:
(42)kkkkkkk∂∂tU1-ΔU1=fLx1,x2,t,x1,x2∈0,.5×0,.5,t∈0,T,kkkkkkk∂∂tU2-ΔU2=fRx1,x2,t,x1,x2∈.5,1×[0,.5],t∈0,T,U1=U2,x1=.5,x2∈0,.5,t∈0,T,n^·∇U1=n^·∇U2,x1=.5,x2∈0,.5,t∈0,T,∂∂x1U1(x1,x2,t)=0,x1=0,x2∈[0,.5],t∈[0,T],U1x1,x2,t=0,x1=1,x2∈[0,.5],t∈[0,T],U2x1,x2,t=0,x1∈0,1,x2=0or.5,t∈0,T,U1x1,x2,0=gLx1,x2,t,x1∈0,1,x2∈0,.5,U2x1,x2,0=gRx1,x2,t,x1∈0,1,x2∈0,.5.
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 fL, fR, gL, and gR 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 Figures10 and 11. The spectral radius of A-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.Spectral radii for multirate example 1 in 2 spatial dimensions.n=1, p=2 (a), and n=1, p=10 (b). Black is A-1, dotted black is A-1Zp, green is M-1N, blue is Tn, and red marks the value one. The horizontal axis shows the size of the large time step.
(a)
(b)Spectral radii for multirate example 2 in 2 spatial dimensions.n=3, p=2 (a), and n=100, p=2 (b). Black is A-1, dotted black is A-1Zp, green is M-1N, blue is Tn, and red marks the value one. The horizontal axis shows the size of the large time step.
(a)
(b)


## 5.1. Multirate Example 1: Varyingp with Constant n

The multirate examples are the same as those described in Section4.2 with only the values of n and p altered. We use the same Jacobi style iteration, so based on (36) we have
(41)M=L100L2,N=0P21∑j=0p-1‍L2-1jP120.
In the first multirate example, n is set to 1, and p is set first to 2 and then to 10. Figure 8 shows the spectral radii of the relevant matrices.Spectral radii for multirate example 1.n=1, p=2 (a), and n=1, p=10 (b). Black is A-1, dotted black is A-1Zp, green is M-1N, blue is Tn, and red marks the value one. The horizontal axis shows the size of the large time step.
(a)
(b)Values ofp 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 M-1N is greater than one, the spectral radius of Tn is greater than one.As the number of iterations,n, 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.


## 5.2. Multirate Example 2: Varyingn with Constant p

In the second multirate example,p is set to 2, and n 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 M-1N is more than one, yet the spectral radius of Tn is less than one. In Figure 9(a), where only 3 iterations are used per large time step, the entire range of dt 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. 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 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.Spectral radii for multirate example 2.n=3, p=2 (a), and n=100, p=2 (b). Black is A-1, dotted black is A-1Zp, green is M-1N, blue is Tn, and red marks the value one. The horizontal axis shows the size of the large time step.
(a)
(b)


## 5.3. Multirate Example in Two Spatial Dimensions

In this section, the multirate tests in Sections5.1 and 5.2 are extended to a problem in two spatial dimensions posed on x1∈[0,1], x2∈[0,.5] with the interface located at x1=.5:
(42)kkkkkkk∂∂tU1-ΔU1=fLx1,x2,t,x1,x2∈0,.5×0,.5,t∈0,T,kkkkkkk∂∂tU2-ΔU2=fRx1,x2,t,x1,x2∈.5,1×[0,.5],t∈0,T,U1=U2,x1=.5,x2∈0,.5,t∈0,T,n^·∇U1=n^·∇U2,x1=.5,x2∈0,.5,t∈0,T,∂∂x1U1(x1,x2,t)=0,x1=0,x2∈[0,.5],t∈[0,T],U1x1,x2,t=0,x1=1,x2∈[0,.5],t∈[0,T],U2x1,x2,t=0,x1∈0,1,x2=0or.5,t∈0,T,U1x1,x2,0=gLx1,x2,t,x1∈0,1,x2∈0,.5,U2x1,x2,0=gRx1,x2,t,x1∈0,1,x2∈0,.5.
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 fL, fR, gL, and gR 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 Figures10 and 11. The spectral radius of A-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.Spectral radii for multirate example 1 in 2 spatial dimensions.n=1, p=2 (a), and n=1, p=10 (b). Black is A-1, dotted black is A-1Zp, green is M-1N, blue is Tn, and red marks the value one. The horizontal axis shows the size of the large time step.
(a)
(b)Spectral radii for multirate example 2 in 2 spatial dimensions.n=3, p=2 (a), and n=100, p=2 (b). Black is A-1, dotted black is A-1Zp, green is M-1N, blue is Tn, and red marks the value one. The horizontal axis shows the size of the large time step.
(a)
(b)


## 6. 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 ofdt 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 matrixTn 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 Tn 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 dt. 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.



---

*Source: 787198-2015-01-14.xml*