{'fraction_non_alphanumeric': 0.04201011614837018, 'fraction_numerical': 0.012012926189584114, 'mean_word_length': 4.227690047741461, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Conservative iterative methods for implicit discretizations of conservation laws

Conservation properties of iterative methods applied to implicit finite volume discretizations of nonlinear conservation laws are analyzed. It is shown that any consistent multistep or Runge-Kutta method is globally conservative. Further, it is shown that Newton's method, Krylov subspace methods and pseudo-time iterations are globally conservative while the Jacobi and Gauss-Seidel methods are not in general. If pseudo-time iterations using an explicit Runge-Kutta method are applied to a locally conservative discretization, then the resulting scheme is also locally conservative. However, the corresponding numerical flux can be inconsistent with the conservation law. We prove an extension of the Lax-Wendroff theorem, which reveals that numerical solutions based on these methods converge to weak solutions of a modified conservation law where the flux function is multiplied by a particular constant. This constant depends on the choice of Runge-Kutta method but is independent of both the conservation law and the discretization. Consistency is maintained by ensuring that this constant equals unity and a strategy for achieving this is presented. Experiments show that this strategy improves the convergence rate of the pseudo-time iterations.


Introduction
Conservation laws arise ubiquitously in the modelling of physical phenomena and their discretizations have been the subject of intense study; see e.g. [17, Chapter 1] and [18,Chapter 1]. They derive their name from the fact that they describe the conservation of some quantities of interest over time.
In computational fluid dynamics (CFD), these quantities are usually mass, momentum and energy. For convenience we always refer to the conserved quantity as the "mass" in this paper.
Many numerical schemes have been designed to mimic mass conservation. Throughout, we refer to such schemes as globally conservative. As special cases, locally conservative schemes dictate that local variations in the mass propagate from one computational cell to neighboring ones without "skipping" any cells (typically referred to as "conservative schemes" in the literature). In particular, finite volume methods are designed in part upon this principle, although many other schemes possess the same property [13]. A major result on the topic is the Lax-Wendroff theorem, see e.g. [16], which provides sufficient conditions for numerical solutions of convergent and locally conservative discretizations to converge to weak solutions of the corresponding conservation law.
It is natural to ask whether such approximate solutions satisfy the conservation properties upon which the discretizations are based, and in particular, whether a Lax-Wendroff type result is available. It is typically claimed that this is not the case, e.g. in [15], where a version of a Newton-Jacobi iteration was noted to violate global conservation. In this paper we set out to answer this question in greater detail for some well-known (families of) iterative methods. We consider conservative implicit space-time discretizations whose solutions in each time step are approximated by a fixed number of iterations with different iterative methods.
After introducing relevant concepts and definitions in section 2, we investigate the global conservation of iterative methods in section 3. It is shown that Newton's method, Krylov subspace methods and pseudo-time iterations using Runge-Kutta (RK) methods are globally conservative if the initial guess has correct mass. However, the Jacobi and Gauss-Seidel methods in general are not. In section 4 we focus on pseudo-time iterations using explicit RK methods. We show that such methods are locally conservative. By fixing the number of iterations and considering the limit of infinitesimal space-time increments, an extension of the Lax-Wendroff theorem is given. It turns out that the numerical solution in general converges to a weak solution of a modified conservation law for which the flux function is multiplied by a particular constant. An expression for this modification constant is given that only depends on the choice of Runge-Kutta method and the selected pseudo-time steps. A technique for ensuring that the constant equals one is presented, thereby ensuring consistency of the resulting scheme. Experiments indicate that the new technique improvess the convergence rate of the pseudo-time iterations. Numerical tests corroborate these findings and further suggest that the results hold for systems of conservation laws in multiple dimensions. Conclusions are drawn in section 5.
Throughout the paper we separate quantities that belong to a spatial discretization from those that do not. Any vector living on a grid with m cells is denoted by a lower case bold letter, e.g. u ∈ R m . Similarly, any matrix operating on such a vector is represented by a bold upper case letter, e.g. A ∈ R m×m . In Section 4, we consider s-stage Runge-Kutta methods. Vector quantities spanning these stages are denoted by bold and underlined lower case letters, e.g. b ∈ R s . Matrices operating on the stages are represented by bold and underlined upper case letters e.g. A ∈ R s×s .

Conservation and conservative discretizations
In this paper we consider numerical methods for conservation laws, meaning partial differential equations of the form possibly with additional boundary conditions. Here, the kth component u k of the vector u ∈ R q represents the concentration of a quantity, i.e. the total amount of that quantity in a domain is given by Ω u k dx. The conservation law states that this amount is only changed by flow across the boundary of the domain. Assuming that the net flow across the boundary is zero, we thus have for all t ∈ (t 0 , t e ] that Ω u(x, t)dx = Ω u 0 (x)dx. (2) In this paper we focus on scalar conservation laws in 1D, i.e. the case when d = q = 1 in (1). Numerical experiments in section 4 suggest that the results of the paper can be generalized to systems in multiple dimensions, however this task is deferred for future work. Further, we restrict our analysis to the Cauchy problem for (1), or the problem with periodic boundary conditions. Thus, the conservation law of interest becomes where either Ω = (−∞, ∞) or Ω = (a, b] and u(a) = u(b). Which case is considered will be clear from context.
A very successful line of research in computational fluid dynamics is to construct numerical methods that respect (2) on a discrete level. Herein, we predominantly consider finite volume methods that utilize the implicit Euler method as temporal discretization. Let the computational grid be given by (x i , t n ) = (i∆x, n∆t) with n ≥ 0. This grid may be either infinite or finite and periodic in space, i.e. x i+m = x i for some positive integer m. We consider discretizations on the form Here, u n i approximates the solution u(x, t n ) of (3) in the computational cell Ω i . The notationf i+ 1 2 (u) is used as a placeholder for a generic numerical flux of the formf (u n+1 i−p , . . . , u n+1 i+q ), where p, q ∈ N determine the bandwidth of the stencil. Multiplying (4) by ∆x and summing over all i reveals that Comparing this with (2) shows that the finite volume scheme discretely mimics conservation of mass.

Local conservation
It is the telescoping nature of the spatial terms in (4) that leads to the preceeding result. In the finite volume community, this property is typically referred to simply as conservation. However, to distinguish it from other mass conservative schemes, we will henceforth refer to this property as local conservation: Definition 1. A discretization of (1) that can be written in the form (4) is said to be locally conservative.
Local conservation plays a central role in the Lax-Wendroff theorem, which considers the Cauchy problem for (3). First, the notion of consistency must be specified: A numerical flux function is inconsistent if it fails to satisfy Definition 2. We will see examples of this in Section 4.
The Lax-Wendroff theorem applies to the Cauchy problem for (3) and considers locally conservative discretizations with consistent numerical flux. If the numerical solution of such a scheme converges to a function u in the limit of vanishing ∆x and ∆t, the theorem provides sufficient conditions for u to be a weak solution of the conservation law (1) [17,Chapter 12]. More precisely, consider a sequence of grids (∆x ℓ , ∆t ℓ ) such that ∆x ℓ , ∆t ℓ → 0 as ℓ → ∞. Let U ℓ (x, t) denote the piecewise constant function that takes the solution value u n i in (x i , x i+1 ] × (t n−1 , t n ] on the ℓth grid. We make the following assumptions: The first of these assumptions asserts that the numerical solution is convergent in the L 1 -norm with limit u. The second one states that the total variation of the numerical solution remains bounded independently of the grid. The Lax-Wendroff theorem may now be stated as follows: Theorem 1 (Lax-Wendroff). Consider the locally conservative discretization (4), suppose that the numerical fluxf is consistent and that Assumption 1 is satisfied. Then, u(x, t) is a weak solution of (1).
The theorem and its assumptions require that the linear or nonlinear systems arising in (4) are solved exactly. However, in practice these systems will be solved approximately using iterative methods. The natural question that now arises is: Do iterative solvers maintain the convergence to weak solutions?

Global conservation
To answer this question we must establish if iterative methods preserve local conservation. However, local conservation is a special case of the more general notion of global conservation. Let us therefore consider a finite and periodic grid. Discretizing (3) in space on a mesh with cells Ω i of volume |Ω i |, we arrive at an initial value problem with u(t),f ∈ R m .
Global conservation implies that (2) is fulfilled on a semidiscrete level by (6): Suppose next that a time integration method is used to solve (6). Note that the conservation property (5) derived for the locally conservative finite volume discretization (4) is a special case of Definition 4, where Ω i = ∆x for every i. Global conservation is thus necessary for local conservation, which in turn is necessary for the Lax-Wendroff Theorem.
In the CFD community, one hears from time to time the claim that implicit methods are not conservative. However, this is not true: All consistent linear multistep methods are globally conservative. To see this, consider a generic s-step linear multistep method s j=0 a j u n+j = ∆t where a j and b j are given method dependent coefficients. Suppose that (7) by |Ω i | and summing over all cells gives Consistent multistep methods satisfy a s = 1 and a 0 + · · · + a s−1 = −1.
Hence, global conservation is ensured.
Similarly, all Runge-Kutta (RK) methods are globally conservative. Consider an s-stage RK method where a j,l and b j are given by the method. Sincef is globally conservative it follows from Definition 3 that the mass of k j is zero and therefore that Hence, RK methods are globally conservative. Exponential integrators have been proven to be conservative as well [9], essentially due to the following useful lemma: Lemma 2. For any y ∈ R m the Jacobianf ′ (u) = ∂f ∂u ∈ R m×m of the globally conservative discretizationf (u) satisfies Proof. See [27].
In the remainder we also need the following result: Lemma 3. Letf ′ (u) ∈ R m×m be the Jacobian of a conservative discretizationf (u) and let α be any scalar. The solution x of a linear system of the form Proof. Multiplying the ith element of (10) by |Ω i | and summing over all i gives By Lemma 2 the final term on the right-hand side is zero.

Globally conservative iterative solvers
We begin by asking: Do iterative solvers respect global conservation? Since global conservation is necessary for local conservation, which in turn is used in Theorem 1, a negative answer will immediately exclude the possibility to extend the Lax-Wendroff theorem to incorporate a particular iterative method. Discretizing (6) in time using the implicit Euler method results in a nonlinear equation system of the form We concern ourselves with a set of iterates u (k) with limit u n+1 as k → ∞.
Definition 5. An iterative method with iterates u (k) is said to be globally conservative if In the following subsections we investigate the conservation properties of some of the most common iterative methods.

Newton's Method
Newton's method finds an approximate solution to (11) by solving the sequence of linear systems for k = 0, 1, . . . Sincef is conservative it has zero mass by Definition 3. If u (k) and u n have the same mass, the entire right-hand side of the first equation in (12) has zero mass. Hence, by Lemma 3, the increment ∆u also has zero mass. Thus, the mass of u (k+1) is the same as that of u (k) , which by assumption is the same as that of u n . Hence, Newton's method is globally conservative. In large scale practical applications, each linear system in (12) is solved approximately using other iterative methods. In the following, we consider a generic linear system that may either represent a Newton iteration or a linear discretization of a conservation law. In line with Lemma 2, it is assumed that the mass of Ay is zero for any vector y ∈ R m .

Stationary linear methods
Any iterative method for solving (13) that can be written in the form for some matrices M and N , is termed a stationary linear method. Well known examples include the Richardson, Jacobi and Gauss-Seidel iterations.

Richardson iteration
For the Richardson iteration, M = I − θ(I − αA) and N −1 = θI, where θ ∈ R is a fixed parameter. With some minor rearrangements, the iteration can thus be expressed as Multiplying the ith element by |Ω i | and summing over i results in By assumption, the matrix-vector product Au (k) has zero mass, hence the final term vanishes. Further, if u (k) has the same mass as b, then the second term on the right-hand side vanishes. By induction it follows that the Richardson iteration is globally conservative.

Jacobi iteration
Consider again the linear system (13). With the matrix decomposition I − αA = I − α(D + L + U ), where D, L and U are the diagonal, lower and upper triangular parts of A respectively, the Jacobi iteration is obtained by setting M = α(I − αD) −1 (L + U ) and N −1 = (I − αD) −1 . Inserting this into (14), multiplying by I − αD, then adding and subtracting αDu (k) , gives after some rearrangement Multiplying the ith element by |Ω i | and summing over i results in The final term vanishes by the assumption. However, the second term on the right-hand side is in general non-zero and thus introduces a conservation error The Jacobi method is therefore not globally conservative in general.
However, consider the special case when the diagonal D of A is a scalar multiple of the identity matrix, say D = aI. Then a i,i = a for i = 1, . . . , m. If u (k) has the same mass as b we find, after rearrangement of (15), that Thus, u (k+1) has correct mass, at least as long as αa = 1. Hence, the Jacobi iteration is conservative in this special case.

Gauss-Seidel iteration
The Gauss-Seidel iteration is defined by M = α(I − α(D + L)) −1 U and N −1 = (I −α(D +L)) −1 . Inserting this into (14), multiplying by I −α(D + L), then adding and subtracting αU u (k+1) gives Following the same procedure as previously, we multiply the ith element by |Ω i | and sum over i to obtain By assumption, the final term vanishes. However, the second term on the right-hand side is in general non-zero. Thus, the Gauss-Seidel method is not globally conservative and the conservation error is given by However, note that for problems where A is lower triangular the Gauss-Seidel iteration is indeed conservative. In fact, for such problems the method is by definition exact.

Krylov subspace methods
Given an initial guess u (0) for the solution of the linear system (13), the kth iteration of a Krylov subspace method belongs to the space Examples of Krylov subspace methods are GMRES, CG and BICGSTAB [24,.
We may express the kth iteration u (k) as for some coefficients c l , l = 0, . . . , k − 1 that depend on, and are chosen by the method. It follows that The final term vanishes by assumption. Note that the mass of r (0) is zero if the mass of u (0) equals that of b. Thus, the second term on the right-hand side also vanishes. All Krylov subspace methods are therefore conservative if the initial guess is chosen to have the same mass as the right-hand side. Typically a preconditioner P −1 will pre-or post-multiply A in Krylov subspace methods. We will not delve into this subject here but merely observe that if y and P −1 y have the same mass for any choice of y, then the above analysis applies without modification.

Multigrid methods
Multigrid methods combine two methods, namely a smoother and a coarse grid correction (CGC). The idea is to separate the residual into high and low frequency components. As a smoother, an iterative method is applied, designed to effectively damp the high frequency components.
The role of the CGC is to remove the low frequency components of the residual. The residual is mapped to a coarser grid using a restriction operator, R. Smoothing is then applied and a prolongation operator P is used to reconstruct the residual on the fine grid. Finally, the fine grid solution is updated by adding the correction to the previous iterate. The procedure can be applied to a hirearchy of grids, thereby resulting in a multigrid method. If we consider only two grid levels and assume that the system is solved exactly on the coarser one, then the coarse grid correction for the linear problem (13) has the form Here, the notation (·) ℓ and (·) ℓ−1 denote matrices operating on the fine and the coarse grid respectively.
The easiest way to make a multigrid method globally conservative is to choose its components to be globally conservative. The smoother can be any of the globally conservative iterations discussed so far. By Lemma 3, the inverted matrices in (16) are mass conserving. Thus, it remains to choose globally conservative restriction and prolongation operators R and P . This is achieved using agglomeration. The restriction is performed by agglomerating a number of neighboring cells by constructing a volume weighted average of fine grid values. Conversely, the prolongation is done by injecting coarse grid values at multiple fine grid points. In fact, this is the standard choice of multigrid method in CFD since it results in faster convergence than other alternatives. For details, see [3].

Pseudo-time iterations
Pseudo-time iterations are obtained by adding a pseudo-time derivative term to the algebraic equations (11) (or (13)), Here, u 0 is initial data that must be provided, e.g. u 0 = u n . The idea is that the system eventually should reach a steady state as τ → ∞ where the pseudo-time derivative vanishes [14], resuting in a solution to the original nonlinear system (11). Any time integration method can in principle be applied to the initial value problem (17). Explicit RK methods are globally conservative by construction. Implicit methods result in new systems of equations whose solutions are once again approximated using iterative methods [5]. Whether the resulting approximation is conservative depends on the choice of method. Any one of the conservative methods discussed so far can be applied in principle. We will return to pseudo-time iterations in Section 4, where more details are provided.

Numerical validation
To validate the results of the preceding section we consider two simple experiments. The first is the linear advection equation with periodic boundary conditions, We discretize using the finite volume scheme (4) with a central numerical This results in a linear system of the form Here, Tridiag(·) refers to a tridiagonal matrix with periodic wrap-around. Note that the diagonal of A is constant. We therefore expect the Jacobi method to be conservative for this particular discretization.
The second experiment is Burgers' equation with periodic boundary conditions, Again, we use the finite volume scheme (4), however with the upwind flux The result is a nonlinear system in the form (11) with . The Jacobian of this discretization is lower triangular except for a single element in the top right corner. Thus, we expect the Gauss-Seidel method to be nearly conservative if the grid is sufficiently fine.
Burgers' equation (19) is solved using Newton's method with initial guess u n . The resulting linear systems are either solved exactly or using the following iterative methods with 0 as initial guess: (R) The Richardson iteration using θ = 0.5.

(J)
The Jacobi method.
(CGC) A two-level coarse grid correction.
As prolongation and restriction operators for the CGC, agglomeration gives The same methods are used directly on the advection problem (18) with u n as initial guess. First we consider very coarse discretizations with ∆x = ∆t = 0.5 for both problems. A single iteration is used with each method and the mass error and residual is computed after one time step. The results are shown in Table  1. Here, any error smaller than 10 −15 is denoted as zero. The Richardson iteration, GMRES, CGC and the pseudo-time iterations are conservative for both problems. Similarly, the exact Newton method (Exact) is conservative. As expected, the Jacobi method is conservative for the advection problem but not for Burgers' equation. Gauss-Seidel is not conservative for any of the two problems.  The total mass is computed in each time step and compared with the mass of the initial data. Fig. 1 shows the results. The main difference to   (18) and Burgers' equation (19).
the previous test is that Gauss-Seidel has a non-detectible mass error for Burger's equation. This is in line with our expectations since the Jacobian only has a single nonzero element above the main diagonal.

Local conservation of pseudo-time iterations
Having covered global conservation in the previous section, we now switch focus to local conservation. The Jacobi and Gauss-Seidel methods cannot be locally conservative since they are not globally conservative in general.
For the other methods considered so far, the question remains open.
In what follows, we restrict our attention to the Cauchy problem for (3). We focus on pseudo-time iterations, show that these methods are locally conservative and prove an extension of the Lax-Wendroff theorem. Noting that the Richardson iteration is equivalent to pseudo-time iterations with explicit Euler, this method is also covered.

Pseudo-time iterations in locally conservative form
We once again consider the finite volume method (4) where the implicit Euler method is used as time discretization. In order to apply pseudo-time iterations to this scheme, we introduce a pseudo-time derivative, where the nonlinear function g i is given by Several different methods are available for iterating in pseudo-time [5,28]. Herein, we use an explicit s-stage Runge-Kutta method (ERK). Let (A, b, c) denote the coefficient matrix and vectors of the ERK method. We denote the kth pseudo-time iterate by u where the stage vectors U (k) As previously, p and q determine the bandwidth of the finite volume stencil.
In the remainder we always use a fixed number N of iterations. A step in physical time is taken by setting u n+1 Recall that the stability function φ(z) of an RK method (A, b, c) is given by (see e.g. [31,Chapter IV where I is the s × s identity matrix. The stability region of the RK method is defined as the subset of the complex plane for which |φ(z)| < 1. The proof of the following lemma is rather lengthy and is therefore deferred to Appendix A.
where the fluxĥ andf (k) Remark 1. The product in (25) is empty when k = N − 1. To handle this case we use the convention It follows immediately from Lemma 4 that the locally conservative nature of the discretization (4)   we see that it is additionally dependent on u n i so that we may writê h To establish the consistency we must therefore show thatĥ N (u, . . . , u; u) = cf (u). To do this, we first note from (20) that g i (u, . . . , u; u) = 0 due to the consistency off i+ 1 2 . Using the fact that A is lower triangular for any ERK method, it follows from (21) and (22) that U  (23) giveŝ where we have utilized the fact that the sum in the third equality is telescoping. Finally, note thatĥ Theorem 6. Consider a sequence of grids (∆x ℓ , ∆t ℓ ) such that ∆x ℓ , ∆t ℓ → 0 as ℓ → ∞. Fix N independently of ℓ, set u (0) i = u n i and terminate the pseudo-time iterations after N steps. Let ∆τ k,ℓ /∆t ℓ = µ k,ℓ = µ k be constants independent of ℓ for each k = 0, . . . , N −1. Suppose that the numerical fluxf in (4) is consistent with f and that Assumption 1 is satisfied. Then, u(x, t) is a weak solution of the conservation law Proof. We follow the proof of the Lax-Wendroff theorem given in [17,Chapter 12] with changes and details added where necessary. Throughout the proof we let U ℓ (x, t) denote the piecewise constant function that takes the solution value u n i in (x i , x i+1 ] × (t n−1 , t n ] on the ℓth grid. Similarly, for j = 1, . . . , s we let U (k) j ℓ (x, t) be the piecewise constant function that takes the value U (k) The discretization (26) can equivalently be expressed as Let ϕ ∈ C 1 0 be a compactly supported test function. Multiply (29) by ϕ(x i , t n ) and sum over all n and i. Due to the compact support of ϕ, these sums can be extended arbitrarily beyond the bounds of x and t, hence we obtain At this point we will use summation by parts on the sums in (30). Recall that for two sequences (a i ) and (b i ), the summation by parts formula can be expressed as Applied to the n-sum in the first term and to the i-sum in the second term in (30), this results in Here we have used the fact that ϕ has compact support in order to eliminate all boundary terms except the one at n = 0. We now let ℓ → ∞ and investigate the convergence of the terms arising in (31). The first and third terms are identical to those in the proof of the Lax-Wendroff theorem [17,Chapter 12]. Thus, it can immediately be concluded that the first term converges to It remains to investigate the second term in (31). Expandingĥ To this expression, add and subtract The bracketed part of this expression evaluates to c(µ 0 , . . . , µ N −1 ), as in the proof of Theorem 5. Thus, the first half of this expression equals Apart from the factor c(µ 0 , . . . , µ N −1 ), which is independent of ℓ, this term appears identically in the proof of the Lax-Wendroff theorem [17,Chapter 12] and converges to It remains to show that (32) vanishes in the limit ℓ → ∞. Since U ℓ and U (k) Here,f 2 (x, t) is placeholder notation for the vector whose jth element iŝ To establish that this term indeed vanishes, it suffices to show that tends to zero for almost every x. By the Cauchy-Schwarz and triangle inequalities, (33) is bounded by where the Euclidean norm in R s is used. The only term in (34) that depends and it therefore suffices to show that this term vanishes for almost every x.
Sincef i+ 1 2 is Lipschitz continuous in each argument there is a constant L such that the final expression is bounded by The second of these terms appear in the proof of the Lax-Wendroff theorem [17,Chapter 12] and vanishes in the limit for almost every x due to the bounded total variation of U ℓ . Further, from (22) and (21) and the fact that u A few remarks about Theorem 6 are in place: First, we demand from a useful iterative method that it converges to the correct solution as N → ∞. Thus, the pseudo-time steps should be chosen in a way that ensures that c → 1, or equivalently, The simplest way to do this is to choose each pseudo-time step so that |φ(−µ l )| < 1, i.e. to stay within the stability region of the RK method. Secondly, observe that if φ(−µ l ) = 0 for any l, then c = 1 irrespective of how many further iterations that are carried out. For some RK methods such a root exists while for others it does not. For instance, the explicit Euler method has stability function φ(−µ) = 1 − µ, hence µ = 1 is a root of φ. On the other hand, Heun's method has stability function φ(−µ) = 1−µ+µ 2 /2 > 0 for all µ ∈ R and therefore does not have any real roots. A strategy is thus to choose a RK method with a root, begin the pseudo-time iterations with a step that corresponds to this root, then resort to a conventional method for choosing the remaining pseudo-time steps. The initial iteration will not change the limit as N → ∞ if the remaining pseudo-time steps are selected within the stability region.

Numerical results
Next we validate Theorem 6 by numerically solving a series of linear and nonlinear conservation laws.

Linear advection
The first setting is the linear advection equation (18). The computational domain is x ∈ (−1, 1], t ∈ (0, 0.25] and periodic boundary conditions are used. The upwind fluxf i+ 1 2 = u n+1 i is used for the spatial discretization. The resulting finite volume method becomes Throughout the experiments the temporal and spatial increments are chosen to be equal: ∆t = ∆x. We validate Theorem 6 by studying the convergence of the numerical scheme to the solution of the original advection problem (18) as well as to the solution of the modified version u t + cu x = 0. As pseudo-time iteration we use the explicit Euler method, Heun's method and the third order strong stability preserving RK method SSPRK3 [26]. Theorem 6 predicts that these methods respectively will modify the propagation speed by the factor Here, we fix N = 4 and consider two different sequences of pseudo-time steps, one with constant and one with variable step sizes. The first is given by µ l = 1/20 and the second by µ l = 2 −l for l = 0, . . . , 3. The corresponding modification constants c(µ 0 , . . . , µ 3 ) are shown in Table 2. Note that with the second sequence, c = 1 for the explicit Euler method. This is due to the fact that µ 0 is a root of the stability polynomial in this case.  The advection problem is solved on a sequence of grids with grid spacing ∆x = 2/(40 × 2 j ) for j = 1, . . . , 12. The L 2 error is calculated with respect to the exact solution of the original and modified equations. The results are shown in Fig. 2. When constant pseudo-time steps are used (Fig. 2a), the lines overlap. None of the methods converge to the solution of the original problem since the numerical solution is moving with the incorrect speed. Instead, all three schemes converge to the solution of their respective modified conservation law. Similar results are seen when variable pseudotime steps are used (Fig. 2b).
Solutions of the advection equation (18) propagate from left to right at unit speed. In this setting it follows from Theorem 6 that the numerical solution will propagate with the speed c(µ 0 , . . . , µ N −1 ) in the limit ∆x, ∆t → 0. Thus, if c = 1 in each physical time step, the numerical solution will drift out of phase. However, it should be noted that Theorem 6 is an asymptotic result and does not necessarily imply that c(µ 0 , . . . , µ N −1 ) alone captures the propagation speed error on coarse grids. In fact, we may generally expect a  contribution to the speed from dispersion errors built into the discretization; see e.g. [19-22, 29, 30] for details and remedies.
In the next experiment we verify that c depends on N as predicted by Theorem 6 by measuring the propagation speedc of the numerical solution while varying N . To measurec we track the x-coordinate of the maximum of the propagating pulse through time. To get accurate measurements we extend the computational domain to x ∈ (−1/5, 1/5], t ∈ (0, 6].
For this experiment we fix ∆x = 0.003 and set µ l = µ to be fixed for each l = 0, . . . , N − 1. The measured propagation speed error 1 −c for different µ are shown in Fig. 3a for Heun's method and in Fig. 3b for SSPRK3. The solid lines show the theoretically predicted error, The measured and theoretical speed errors agree well for µ = 0.05 and µ = 0.2. For µ = 0.5 a discrepancy between theory and measurement is seen for small errors. This suggests that the dispersion error intrinsic to the finite volume scheme is starting to dominate the propagation speed error.

Burgers' equation
Next we consider a triangular shock wave propagating under the 1D Burgers' equation with periodic boundary conditions:  The exact solution to this problem is given by Modifying the conservation law to u t + (cu 2 /2) x = 0 while retaining the same initial condition modifies the exact solution to For this problem we therefore expect that pseudo-time iterations to affect both the speed and the amplitude of the shock front. As for the advection problem, we investigate the L 2 convergence of the numerical scheme to the original and modified conservation laws. We once again use an upwind numerical flux and implicit Euler in time; We run the simulation to time t = 0.1 with ∆t = ∆x using the sequence of grids ∆x = 1/(25 × 2 j ) for j = 1, . . . , 15. The explicit Euler method is used as pseudo-time iterator with N = 12 and µ l = 1/4 for l = 0, . . . , 11. The corresponding modification constant is c ≈ 0.9683. Fig. 4a shows that the numerical solution converges to the solution of the modified conservation law as expected.  Fixing ∆x = 0.004 and setting µ l = µ = 1/4, we run the simulation to time t = 1 with different choices of N . Fig. 5 shows the numerical solutions together with the initial data (dashed line) and the exact solution (dotted line). The locations of the tips of the shock waves as predicted by Theorem 6 are indicated as crosses in the figure. There is good agreement between theory and experiment, although the shocks appear slightly smeared due to the built-in dissipation in the numerical scheme.
Next, we repeat the experiments but change the initial data to a step function, Since this data is not periodic, we impose the boundary condition u(0, t) = 1. Theorem 6 does not treat boundary conditions and it is therefore interesting to see if it still provides useful predictions of the behaviour of the numerical solution in this setting. The exact solution is the initial step function travelling to the right. The shock speed is given by the Rankine-Hugoniot condition as where u l,r denote left and right states of the discontinuity respectively. The shock speed of the modified conservation law is instead c/2. We use the exact same grids and pseudo-time iterations as for the triangular shock and measure the L 2 error with respect to the exact and modified conservation laws. The results are shown in Fig. 4b. Convergence is once again seen towards the modified equation, thus verifying that Theorem 6 predicts the propagation speed error correctly despite the added boundary condition.
Next we extend the time domain to t ∈ (0, 1], fix ∆x = ∆t = 4∆τ = 1/100 and vary the number of iterations, N . The computed solutions using N = 1, N = 3 and N = 12 are shown in Fig. 6 together with the predicted shock locations (dashed lines). Once again, there is good agreement between prediction and experiment, although numerical dissipation smears the shock fronts somewhat.
As mentioned previously, the shock speed error can be eliminated entirely for the explicit Euler method by noting that µ = 1 is a root of the stability polynomial 1 − µ. To highlight the effect of this, we introduce two strategies for choosing the pseudo-time steps:  The first strategy is the same that we have used in the experiments so far.
The second one ensures that c = 1 by taking a large initial pseudo-time step. Note that both strategies correspond to integration in pseudo-time to the same point; τ = 3∆t. The relative residuals of the pseudo-time iterates in the first physical time step are shown in Fig. 7 for the two strategies. The large initial pseudo-time step in Strategy 2 results in a considerably greater residual reduction than the corresponding iteration using Strategy 1. Interestingly, subsequent iterations yield faster convergence of the residual using Strategy 2 as seen by the steeper gradient. This suggests that the incorrect shock speed makes the dominant and slowest converging contribution to the residual for this problem when using Strategy 1. We also conclude that the point to which we march in pseudo-time, here τ = 3∆t, have less of an impact on the convergence than the choice of pseudo-time steps used to reach this point.

The Euler equations
As a second nonlinear problem we consider the 2D compressible Euler equations,  where γ = 1.4 and e = E − (u 2 + v 2 )/2 is the internal energy density. The domain is taken to be periodic in both spatial coordinates. The setting is the isentropic vortex problem [25] with initial conditions where r = 1 − x 2 − y 2 . Here, ǫ = 5 is the circulation and M ∞ = 0.5 is the Mach number. As the solution evolves in time, the initial vortex propagates in the horizontal direction with unit speed. As in previous experiments, we use implicit Euler in time, yielding a finite volume scheme of the form Along the x-coordinate we use a fourth order centered flux and similarly along the y-coordinate. With this choice, the resulting problem violates the assumptions of Theorem 6 in three ways: (i) It is 2D, (ii) it is a system of equations, and (iii) the scheme is not total variation bounded. The exact solution of the isentropic vortex problem is given by the initial data (centred at the origin) translated to the point (x, y) = (T, 0), where T is the end point of the time domain. However, the exact solution of the modified conservation law is instead centred at (x, y) = (cT, 0).
The explicit Euler method is again used as pseudo-time iteration. We study the convergence of the numerical solution to the exact solutions of the original and modified conservation laws by measuring the L 2 error of the density component. The discretizations are chosen such that ∆x = ∆y = 4∆t on all grids. Two different strategies for choosing the pseudo-time steps are considered: Both strategies integrate in pseudo-time to the point τ = 1.8∆t in each time step. However, Theorem 6 predicts that Strategy 1 will give a speed modification c(µ 0 , . . . , µ 8 ) ≈ 0.866. On the other hand, Strategy 2 will give c(µ 0 , . . . , µ 4 ) = 1, i.e. the correct propagation speed, due to the large initial pseudo-time step. Fig. 8 shows the numerical solutions at time T = 10 for the case where ∆x = ∆y = 0.2, ∆t = 0.05. The predicted and observed vortex locations agree very well. Fig. 9a shows the convergence of the numerical solutions, measured at time T = 1. Convergence to the correct solution is observed for Strategy 2 (S2) but not for Strategy 1 (S1). This is expected due to the incorrect location of the vortex in the latter case. However, convergence to the solution of the modified conservation law is seen (S1 mod). Thus, Theorem 6 accurately predicts the behavior of the numerical solution despite the violated assumptions.
In practical applications, the pseudo-time iterations are terminated when the residual has decreased beneath some tolerance. It is interestig to see in what way the convergence is affected by the choice of strategy. Returning to the setting in Fig. 8, the residual in each pseudo-time iteration for all 200 physical time steps are shown for the two strategies in Fig. 9b. The 200 lines overlap nearly perfectly, suggesting that the residual behaves similarly in each physical time step. The initial iteration in Strategy 2 evidently has a large impact on the reduction of the residual that supercedes those of the other iterations put together. The remaining iterations appear to reduce the residual by comparable amounts for the two strategies, as seen by the similar gradients. In contrast to the shock problem considered previously, this suggests that the convergence rate of the residual is dictated by other   factors than the propagation speed for this particular problem. Nonetheless, a correct propagation speed is visibly very beneficial, here with a drop in relative residual of more than an order of magnitude.

Summary and Conclusions
In this paper we have studied conservation properties of a selection of iterative methods applied to 1D scalar conservation laws. The fact that conservation is a design principle behind many numerical schemes motivates such a study. We have established that Newton's method, the Richardson iteration, Krylov subspace methods, coarse grid corrections using agglomeration, as well as ERK pseudo-time iterations preserve the global conservation of a given scheme if the initial guess has correct mass. However, the Jacobi and Gauss-Seidel methods do not, unless the linear system being solved possesses particular properties. The stronger requirement of local conservation has been investigated for ERK pseudo-time iterations. We have shown that local conservation is preserved for finite volume schemes that employ the implicit Euler method in time. However, the resulting modified numerical flux may be inconsistent with the governing conservation law. An extension of the Lax-Wendroff theorem shows that this inconsistency leads to convergence to a weak solution of a conservation law modified by a particular constant. We give an exact expression for the modification constant, which depends only on the stability function of the ERK method and the pseudo-time steps. Depend-ing on the problem solved, this modification can alter both the propagation speed and amplitude of the numerical solution if the constant differs from unity. We present a strategy for ensuring that the constant equals one and show numerically that the strategy results in faster convergence. Experiments suggest that the results hold even for systems of conservation laws in multiple dimensions.
Proof. Note first that (21)- (22) can be equivalently written in the form . . , −1, 0, 1, . . . , where 1 = (1, . . . , 1) ⊤ ∈ R s . Here, we use the notation g i (U k ) = g i (U holds. We will show that it also holds for N + 1. From (20), (24) and (40) it follows that Solving for U Note that (I + µ N A) −1 exists since A is lower triangular. Evaluating g i U (N ) using (20) gives Rearranging, using the induction hypothesis and the expression (25)  It remains to show that the lemma holds when N = 1. To this end, recall that u (0) i = u n i and note from (20) and (40) that Using (20) it follows that Here we have once again used (41) in the final equality. Thus, u (1) i can be evaluated using (40) as Division by ∆t and rearrangement shows that the lemma holds when N = 1. By induction it holds for all N ≥ 1.