Collocation Methods for High-Order Well-Balanced Methods for Systems of Balance Laws

In some previous works, two of the authors introduced a technique to design high-order numerical methods for one-dimensional balance laws that preserve all their stationary solutions. The basis of these methods is a well-balanced reconstruction operator. Moreover, they introduced a procedure to modify any standard reconstruction operator, like MUSCL, ENO, CWENO, etc., in order to be well-balanced. This strategy involves a non-linear problem at every cell at every time step that consists in finding the stationary solution whose average is the given cell value. In a recent paper, a fully well-balanced method is presented where the non-linear problems to be solved in the reconstruction procedure are interpreted as control problems. The goal of this paper is to introduce a new technique to solve these local non-linear problems based on the application of the collocation RK methods. Special care is put to analyze the effects of computing the averages and the source terms using quadrature formulas. A general technique which allows us to deal with resonant problems is also introduced. To check the efficiency of the methods and their well-balance property, they have been applied to a number of tests, ranging from easy academic systems of balance laws consisting of Burgers equation with some non-linear source terms to the shallow water equations—without and with Manning friction—or Euler equations of gas dynamics with gravity effects.


Introduction
Let us consider a PDE system of the form: where U(x, t) takes values on an open convex set Ω ⊂ R N , f : Ω −→ R N is the flux function, S : Ω −→ R N , and H is a continuous known function from R to R (possibly the identity function H(x) = x). It is supposed that system (1) is strictly hyperbolic, that is, D f (U) = ∂ f ∂U (U) has N real distinct eigenvalues r 1 (U), · · · , r N (U) and associated eigenvectors v 1 , · · · , v N .
Systems of the form (1) have non-trivial stationary solutions that satisfy the ODE system: A numerical method is said to be well-balanced if it has discrete stationary solutions that approximate all the stationary solutions of the system or, at least, a relevant family of them. The design of numerical methods that have the previous property is of great importance in problems where small perturbations of a stationary solution are going to be considered, like in many geophysical problems. Methods that are well-balanced have been proposed by many authors: for example, References  and their references.
The aim of this work is to introduce a general methodology to build well-balanced high-order finite-volume numerical methods for (1) of the form: where: are the computational cells, whose length ∆x is supposed to be constant for simplicity; • U i (t) is the approximation of the average of the exact solution at the ith cell at time t, that is, U(x, t) dx; , where F is a consistent numerical flux, and U t,± i+ 1 2 are the reconstructed states at the intercells, i.e., where P t i (x) is a reconstruction of order p of the solution at the given ith cell, computed from the sequence { U i (t)}: where S i is the set of cell indices belonging to the stencil associated to P t i (x); • finally, Given a function U, the following notation will be used to represent its cell averages and the approximations to its cell averages or its point values at the intercells: Following References [36,37], the well-balanced property of the methods can be transferred to the reconstruction operator: Definition 1. A reconstruction operator P i (x) is said to be well-balanced for a function U * if where P i is the approximation of U * computed from the vector of cell-averages {Ū * i } of the function U * .
In Reference [37], it has been proved that the numerical method (3) is exactly wellbalanced according to Definition 2 if the reconstruction operator is well-balanced for every stationary solution U * and Definition 2. The numerical method (3) is said to be exactly well-balanced if the sequence of cell-averages {Ū * i } (or the sequence of their approximations { U * i } if a quadrature formula is used to compute them) of any stationary solution U * of (1) is an equilibrium of the ODE system (3).
Furthermore, in References [36,37], the authors propose a method to build a wellbalanced reconstruction operator from any standard one that will be recalled in the next section. As we shall see, the main difficulty in the application of this technique comes from the fact that, at every cell and at every time step, one has to solve a non-linear problem of the form: Find U such that where U i is an approximation of the average at the i-th cell of the sought solution of (1). Once this problem with given average has been solved at the cell, the solution of (7), that will be denoted by U * i , has to be extended to the whole stencil by solving two Cauchy problems of the form More precisely, (8) with initial condition U(x i+1/2 ) = U * (x i+1/2 ) (resp. with final condition U(x i−1/2 ) = U * (x i−1/2 ) ) has to be solved forward in space (resp. backward in space) in the cells at the right (resp. at the left) of the ith cell in its stencil.
The previous procedure has been applied to balance laws whose stationary solutions are known in implicit or explicit form, so that problems (7) and (8) can be easily solved: see, for example, References [32,38], and the references therein. However, when solving these problems by analytic procedures is not possible, they have to be solved numerically. Please note that the two main difficulties come from the condition on the average in (7) and from the fact that the ODE is not in normal form: special care has to be taken in resonant problems i.e., problems in which the Jacobian matrix of the flux, D f (U) becomes singular. When the Jacobian is regular, the ODE system (2) can be written in normal form and the solution of the problems is easier. At a sonic point, i.e., a point in which one of the eigenvalues vanishes, the ODE system may have no solution or may have more than one. In recent papers [39,40], a strategy based on the interpretation of (7) as a control problem has been proposed to deal with the first difficulty: the value of U * i at the leftmost intercell of the stencil is the control that has to be selected so that the condition on the average is satisfied. Newton's method is used to find the control, where the gradient is computed by solving the adjoint problem. A standard ODE solver, namely RK4, is used to solve the state and the adjoint equations, and to extend the solution to the whole stencil.
The first goal of this paper is to present a general framework to design and analyze well-balanced high-order numerical methods in which problems (7) and (8) are numerically solved: the definition of the well-balanced property and a general result that allows one to check it will be stated. The second goal is then to introduce a methodology based on the resolution of (7) and (8) using RK collocation methods that leads to numerical methods that are well-balanced according to the given definition. As we shall see, the condition on the average in (7) will naturally lead to non-linear problems that will be solved by an iterative procedure. Finally, the third goal is to propose a strategy to deal with resonant problems that is general, albeit problem-dependent: this strategy will be illustrated in the particular case of the shallow water system.
The organization of the article is as follows: Section 2 is devoted to state the general framework. Definitions of well-balanced methods and reconstruction operators adapted to the numerical solution of the local problems will be introduced and a general result will be stated showing that well-balanced reconstruction operators lead to well-balanced methods. In Section 3, we introduce a strategy based on the use of RK collocation methods to solve problems (7) and (8), and we show that it leads to numerical schemes that are well-balanced. Section 3 ends with the introduction of the technique that allows us to deal with resonant problems.
In Section 4, a number of numerical tests are presented to analyze the performance of the methods and to check their accuracy and well-balancedness. Some numerical tests for scalar and systems of balance laws are considered: Burgers' equation with a non-linear source term, the shallow water system-both with and without Manning friction-or the Euler equations of gas dynamics including gravity effects are considered. In addition, we also present some numerical results to check that such numerical methods introduced in this work behave correctly in the presence of critical states. Moreover, we show that the numerical methods are able to preserve subsonic and supersonic moving stationary solutions for the compressible Euler equations with gravitational force. Finally, some conclusions are drawn in Section 5, and further developments are also discussed.

Well-Balanced Methods
The following strategy to obtain a reconstruction operator that is well-balanced according to Definition 1, from a standard one Find, if possible, a stationary solution U * i (x) defined in the stencil of cell I i , ∪ j∈S i I j , such that: 1 ∆x Otherwise, take U * i ≡ 0.

2.
Compute the fluctuations {V j } j∈S i given by and compute the reconstruction operator: According to Reference [36], the reconstruction operator P i is well-balanced if Q i is exact for the null function, conservative if Q i is conservative 1 ∆x and P i is of the same order of accuracy p of Q i provided that the stationary solutions are smooth.
Observe that, if it is impossible to find a stationary solution defined in the stencil that satisfies (10), then the standard reconstruction is used. Please note that this choice does not spoil the well-balanced character of the numerical method: in this case, the cell values in the stencil cannot be the averages of a stationary solution (otherwise, there would be at least one solution U * i ); thus, there is no local equilibrium to preserve. On the other hand, if there is more than one stationary solution defined on the stencil that satisfies (10), a criterion is needed to select one of them. This criterion depends on the particular problem and it will be discussed later.
Notice that the use of quadrature formulas for the computation of the integral (3) appearing in the right-hand side of (4) may destroy the well-balanced property of the method. In Reference [37], the authors proposed to write the source term as follows: Again, it can be proved that this reconstruction operator is well-balanced and leads to numerical methods that are well-balanced according to Definition 2.
The main difficulty to implement the well-balanced reconstruction procedure given by Algorithm 2 comes from the first step in which the values at the quadrature points and the intercells of the stationary solution that satisfies (15) have to be computed in the stencil. When it is not possible to solve these problems by analytical procedure, a numerical one is required. In this case, the generic local discrete problem to be solved is the following: Problem 1 (Local problem (LP)). Given an index i and a state W ∈ Ω, find approximations Although the approximations to be found depend in general on ∆x and W, this dependency will not be explicitly written to avoid an excess of notation.
Once a numerical solver has been selected for the local discrete problems, the approximations of the stationary solution satisfying (16) at the quadrature points U * ,m i,j and at the intercells U * ,i±1/2 i are used to define the reconstruction operator as follows (Algorithm 3): Algorithm 3 Well-balanced methods. Reconstruction procedure.
Given a sequence of cell values { U i }, at every cell I i :

1.
Apply, if possible, the local solver at the i-th cell with W = U i to obtain

2.
Compute the fluctuations {V j } j∈S i given by and compute:

Define
Then, the numerical method is defined by (3) with where the dependency on t has not been explicitly written to avoid an excess of indices. The question is now what is the well-balanced property satisfied by the previous method. The following definition is considered: Definition 3. The numerical method (3) is said to be well-balanced with order q ≥ p if for every stationary solution U * of (1), and, for every ∆x, there exists a discrete stationary solution { U * ∆x,i }, i.e., an equilibrium of (3), such that In order to illustrate the definition, let us consider some examples. We consider first the equation: whose stationary solutions are The cell-averages will be computed by the mid-point value so that 1 ∆x We consider the upwind numerical flux and different discretizations of the source term: This is the first-order well-balanced numerical method based on the well-balanced piecewise constant reconstruction operator described in Algorithm 2. The discrete stationary solutions satisfy It is clear then that the set of discrete stationary solutions is given by the point values of the continuous ones: and the numerical method is the exactly well-balanced according to Definition (2). • Method 2: This is the first-order well-balanced numerical method based on the RK2 collocation method that will be introduced in Section 3. The discrete stationary solutions satisfy which leads to The set of discrete stationary solutions is given by: One has for every C: so that the numerical method is well-balanced with order 2. • Method 3: This is the numerical method that corresponds to the upwind treatment of the source term (see Reference [41]). The discrete stationary solutions satisfy which leads to The set of discrete stationary solutions is given by: One has for every C: so that the numerical method is well-balanced with order 1. • Consider now the equation whose stationary solutions are and u ≡ 0. Let us consider the numerical method: Let us see that there are no discrete stationary solutions that change their sign from negative to positive (as it happens with all the stationary solutions) Let us suppose that there exists a discrete stationary solution such thatū i−1/2 < 0, u i+1/2 > 0. Then, at the i-th cell, we would have: which is a contradiction.

-
Let us suppose now thatū i−3/2 < 0,ū i−1/2 = 0,ū i+1/2 > 0. Then, at the i-th cell, we have Sinceū i−1/2 = 0, we have that u i−1 = 0, as well. Therefore, at the (i − 1)-th cell, we have which is again a contradiction. Therefore, it is not possible to find discrete stationary solutions that change from negative to positive. The numerical method is not well-balanced in spite of the fact that it is consistent and stable.
Remark 1. Some authors impose q > p in the definition of well-balanced, i.e., stationary solutions have to be approximated with enhanced accuracy. Please note that, although any standard numerical method of order p will provide approximations of any smooth stationary solution with order of accuracy p, these approximations are not necessarily discrete stationary solutions. This is the case of the fourth example above: although the numerical method will provide first order numerical approximations of the stationary solutions of the form u * (x) = x + C with C < 0, these approximations cannot be discrete stationary solutions.

Remark 2.
From the practical point of view, it is important to know how to compute discrete stationary solutions in order to 'prepare well' the initial conditions. Think, for instance, of a numerical experiment for Equation (22) with initial condition where u * (x) is a stationary solution, say u * (x) = e x , and δ(x) a small perturbation with compact support. If the initial condition is computed by and Method 1 is used, the stationary solution will remain unperturbed in regions where the wave generated by the initial perturbation is not arrived. Nevertheless, if Method 2 is used, the numerical solution will move everywhere to fit a discrete stationary solution of the problem. If the amplitude of the initial perturbation is of order O(∆x), the wave can be lost among the numerical errors at the beginning of the experiment. If, instead, the initial condition is computed using the discrete stationary solution that approximates u * , the numerical solution will remain again unperturbed until the arrival of the wave.
According to Definition 3, in order to be well-balanced, a numerical method has to have discrete stationary solutions. Let us discuss what are the equilibria of (3)-(20) when the well-balanced reconstruction operator given by Algorithm 3 is used. To do this, let us first introduce the following discrete version of Definition 1: Definition 4. The reconstruction operator defined by Algorithm 3 is said to be well-balanced for a sequence of cell values { U i } if, for every i, where U * ,m i,j , U * ,i±1/2 i is the output given by the solver of local problems for W = U i at the ith cell.
In other words, the reconstruction procedure defined by Algorithm 3 is well-balanced for a sequence { U i } if (a) the cell values are recovered when the quadrature formula is applied to the approximations at the points given by the local solver and (b) the approximations at the intercells do not depend on the stencil.
Like in the continuous case, there is a direct relation between discrete stationary solutions and well-balanced reconstruction operators: Theorem 1. If the reconstruction operator described in Algorithm 3 is well-balanced for a sequence of cell values { U * i }, then the sequence is a discrete stationary solution, i.e., an equilibrium of the ODE system given by (3)- (20).
Proof. Let us show that the sequence { U * i } is an equilibrium of (3)- (20). Due to the well-balancedness of the reconstruction operator, at the ith cell, we have: so that Q i (x) = 0. Therefore, As a consequence, (20) reduces to On the other hand: In the previous expression, the consistency of the numerical flux and (28) have been used. Therefore, the right-hand side of (3) vanishes.
The previous result gives the condition under which a sequence of cell values defines a discrete stationary solution. In order to prove that the method is well-balanced, one has to prove that, given any stationary solution U * and ∆x, it is possible to construct a discrete one that approximates U * with order of accuracy greater or equal than p. If we observe the three first examples above, we see that the discrete stationary solutions satisfy a consistent and well-defined ODE solver of the equation satisfied by the continuous ones u x = u, the exact solver for Method 1, the trapezoidal scheme for Method 2, and forward Euler scheme for Method 3. This is not the case for the fourth example. The following result states that, if there exists an ODE solver of (2) with order of accuracy greater or equal than p that is consistent in some sense with the local problems solver, then the numerical method is well-balanced: Theorem 2. Let us suppose that there exists an ODE solver of (2) with order of accuracy q ≥ p that, given any stationary solution U * and any value of ∆x, provides approximations of the values of U * at the quadrature points and the intercells in such a way that, at every cell I i , the solution of the local discrete problem (16) with W equal to is given by Then, the numerical method (3)-(20) is well-balanced.
Proof. The sequence { U * ∆x,i } given by (29) approaches the cell values of U * with order at least p since both the orders of the ODE solver and the quadrature formula are greater or equal than p. On the other hand, it is obviously well-balanced for the reconstruction operator since (27) and (28) are satisfied. As a consequence of Theorem 1, { U * ∆x,i } is a discrete stationary solution, which completes the proof.
In a recent paper [39], the ODE solver chosen for (2) was RK4. The local problems were interpreted as control problems in which both the state and the adjoint local problems were solved using RK4, as well, which gave the consistency of both solvers and, thus, the well-balanced property of the method.

Collocation Methods
In this paper, we propose to build well-balanced numerical methods for (1) on the basis of collocation RK methods. Remember that these methods admit a double interpretation: on the one hand, they are standard RK methods based on a Butcher tableau .
Given an index i 0 , when these methods are applied to a Cauchy problem in a uniform mesh of nodes . . , the numerical solutions are updated as follows: where Here, K 1 i , . . . , K s i solve the non-linear system where On the other hand, they can be interpreted as follows: where P i is the only polynomial of degree s that satisfies: We will consider here Gauss-Legendre methods, in which x 1 i , . . . , x s i and b 1 , . . . , b s are, respectively, the quadrature points and the weights of the Gauss quadrature formula in the interval [x i−1/2 , x i+1/2 ]. This quadrature formula will be used to compute the averages at the cells. The order of accuracy of these methods is 2s.
Because of this double interpretation, it can be easily shown that Gauss methods are symmetric or reversible in the following sense (see Reference [42]):

Cauchy Problem Solver
Let us describe how the collocation RK methods will be adapted to solve problems of the form (8). The following algorithm will be used to approximate the solution of a Cauchy problem with initial condition forward and backward in space (Algorithm 4): Algorithm 4 Numerical solver for the Cauchy problem (8) using collocation RK methods.
With the notation of Theorem 2, this ODE solver provides the following approximations of a stationary solution U * at the quadrature points and intercells: . . , s, ∀i, and the discrete stationary solutions will be then given by The non-linear systems (38) and (39) are solved by using a fixed-point algorithm. If the Jacobian of the flux function is never singular, this algorithm is equivalent to apply the collocation RK method to the ODE system in normal form (9). When the Jacobian may become singular, special care has to be taken when solving the linear systems involved by the fixed-point algorithm: this difficulty will be discussed in Section 3.6.

Local Problem Solver
Let us suppose that the stencils are Remember that local problems (LP) consist in, given an index i and a state W ∈ Ω, finding approximations To solve the problem, we look first for an approximation of U * i in the cell I i using the collocation method, i.e., we look for a vector valued polynomial P i of degree s such that: Observe that this system has (s + 1) vector unknowns (the coefficients of the polynomials) and (s + 1) vector equations.
A fixed-point algorithm will be used to solve this system. Once it has been solved, the approximations of U * i at the intercells are computed by evaluating the polynomial P i : From these two values, the Cauchy solver described in the previous section will be used to compute the approximations of the stationary solution in the other cells of the stencil. Using the standard RK notation, the local problem solver is as follows (Algorithm 5): Algorithm 5 Numerical solver for the local problems (LP) using collocation RK methods.
• Apply Algorithm 4 backward in space from U i−1/2 to obtain for j = i − l, . . . , i − 1 The output of the solver is then: A fixed-point algorithm is used again to solve the non-linear system (42). A sensible initial guess is given by: where K is the solution of the linear system: where x i represents the mid-point of the cell I i . Let us mention again that special care has to be taken if the Jacobian of the flux function may become singular: this difficulty will be discussed in Section 3.6.

Well-Balanced Property
To prove that the high-order numerical methods are well-balanced, we have to check that the ODE solver and the local problem solver described in the two last subsections are consistent in the sense of the statement of Theorem 2. Let be the approximations of a stationary solution U * at the quadrature points and intercells provided by the ODE solver described in Section 3.1 in the uniform mesh of nodes x i+1/2 = x i−1/2 + ∆x, i = i 0 , i 0 + 1, . . . . Let P * i be the only polynomial of degree s such that: By the double interpretation of the method, we have: To solve now the local problem at the i-th cell with we have to look for a polynomial P i satisfying (41). Obviously, P * i solves this non-linear system; thus, P i = P * i , which implies, in particular, The equalities are then trivial, since these values are computed by applying Algorithm 4 forward in space from the same value U * ∆x,i±1/2 , and the equalities follow from the symmetry of the collocation method. Therefore, the hypothesis of Theorem 2 is satisfied, and the numerical method is then well-balanced provided that every stationary solution U * can be approximated by the RK collocation method.

First and Second Order Methods
For second order methods, the MUSCL reconstruction that uses 3-cell centered stencils (l = r = 1) and the second-order 1-stage collocation RK method whose Butcher tableau is will be used. The solution of local problems is particularly simple, since (42) reduces to find Therefore, it is enough to solve and then define Then, the standard method is used forward and bakward in space to compute U 1 −1 and U 1 1 . Finally, for first order numerical method, the trivial pointwise constant reconstruction is considered so that the stencils for the reconstructions consist only of one cell (l = r = 0). The second order collocation method (47) is also used.

Remark 3.
Notice that, since the 1-stage Gauss-Legendre method is second order, the schemes considered here are: • First order numerical methods which are well-balanced with order two. • Second order numerical methods which are well-balanced with order two.

Third-Order Methods
The third-order CWENO reconstruction operator will be used here to design thirdorder methods. This operator uses 3-cell centered stencils, so that l = r = 1. The fourthorder 2-stage Gauss-Legendre method corresponding to the choices, will be used to solve both the local and global problems.

Remark 4.
Notice that, since the 2-stage Gauss-Legendre method is fourth order, we consider third order numerical methods which are well-balanced with order four.

Resonant Problems
When a fixed-point algorithm is used to solve (38), (39), or (42), one or more linear systems of the form have to be solved at every step, where W is a known state, andx a quadrature point. If W is a sonic state, i.e., if one of the eigenvalues of D f (W) vanishes, the problem is said to be resonant and the system may not have solution or it may have infinitely many ones: does not belong to the image of the linear application defined by the matrix D f (W), the system has no solution.

2.
Otherwise, the system has infinitely many solutions where K * is a particular solution, and R is an eigenvector associated to the null eigenvalue.
Since, when applying RK collocation methods, we are looking for solutions of the differential Equation (2), a solution K of (50) is said to be admissible if there exists a C 1 stationary solution U * such that In case 1, there are no admissible solutions. If this situation arises when solving a local problem, it will be assumed that the problem has no solution and the standard reconstruction will be used in the corresponding cell.
In case 2, if K is an admissible solution and W is an isolated sonic point of the corresponding smooth stationary solution U * , the following equality is satisfied: for some ε > 0. Therefore, the limit has to exist. If the available information about the PDE system and its stationary solution allows one to compute this limit by any analytic procedure, then the possible values of (51) constitute the set of admissible solutions, and the solution of (50) will then be selected in this set. This technique will be illustrated in the particular case of the shallow water equations in Section 4.3.

Remark 5.
Observe that, if W is sonic and S(W) does not belong to the image of D f (W), case 2 can only happen if H x (x) = 0. In this case, if there exists a C 1 stationary solution U * such that U * (x) = W, then U * (x) would be an admissible solution of (50) which implies that necessarily H x (x * ) = 0, i.e., a smooth stationary solution of (1) can only reach a sonic state W such that S(W) does not belong to the image of D f (W) at a critical point of H. In this case, (51) will imply a 0/0 indeterminate form that can be solved by L'Hôpital rule. This is the case for many systems, including the shallow water equations, for which a smooth stationary solution can only reach a sonic (critical) state in the minimum points of the depth function H.

Numerical Experiments
Let us describe first the choices that have been made to implement the well-balanced scheme previously described: • The second-order 1-stage collocation RK method is chosen for the first and second order well-balanced numerical schemes and the 2-stage Gauss-Legendre collocation method for third order ones. • Concerning the quadrature rules, the midpoint rule for first and second order schemes are used, and the 2-point Gauss quadrature rule for third order schemes.
In all the numerical tests, Rusanov numerical flux is chosen and MUSCL (see Reference [43]) and CWENO reconstructions operators (see References [44,45]) are used, respectively, for the second and third order schemes. Finally, TVD Runge-Kutta schemes of first, second, and third [46] are used to integrate the ODE system (3)- (20).
In those problems where the initial condition is a steady state U * , the approximation of its cell averages are computed using a quadrature formula if the exact solution is known, or by approximating the solution of the Cauchy problem (2) with initial condition U(a) = U * (a), using the RK collocation method: see Remark 2.
Note that the procedure that we have been describing in this work is quite general and works for a wide class of f , S, and H.
In this section, we will denote the different numerical methods that we compare as follows: • NWBMi, i = 1, 2, 3: a standard non-well balanced numerical scheme of order i. The numerical tests have been run on a laptop with Intel(R) Xeon(R) CPU E3-1220 v3 @ 3.10 GHz with 8 Mb cache using one single core.

Problem 1: Burgers Equation with a Non-linear Source Term I
We consider the Burgers equation with a non-linear source term: This problem corresponds to the choice The stationary solutions satisfy the ODE whose general solution is u(x) = Ce x , C ∈ R.
As the expression of the steady states is known, the stationary solution u * i satisfying (15) could be easily computed for any quadrature formula Although in this test WBMi can be easily computed, we also consider CWBMi and CLWBMi to compare the efficiency of the methods. In what follows, the tolerance considered to stop Algorithm 5 is ε = 10 −15 .
Note that, at the discrete level, the cell-averages of u 0 (x) = e x have to be computed. We consider 4 different ways of computing them: (a) Exactly computed: (b) Using a quadrature formula: (c) The cell averages are approximated by where u * ,i h,m are the approximation the Cauchy problem computed with the standard RK4 method on the mesh selected to solve the control problems described in Reference [39]. The numerical methods are integrated until the final time t = 5, which is more than twice the travel time of a small perturbation of the steady state along the domain. L 1 errors between the initial and final cell-averages will be computed. According to Theorem 1, errors should be of the order of the machine precision (MP) if option (d) is considered to compute the initial condition for CLWBMi, i = 1, 2, 3. Figure 1 (right) shows the errors corresponding to CLWBMi, i = 1, 2, 3 for a 200-cell mesh. As expected, they are of the order of MP.
Notice that, in the other cases, if (a), (b), or (c) are used for CLWBMi, i = 1, 2, 3, the errors decrease with ∆x (see Table 1). Thus, if option (b) is chosen and a third order reconstruction operator is used, then the convergence rate is expected to be 4 as errors are generated by the use of the 2-stage collocation RK method for solving the local problems. If option (c) is used, the same convergence rate is expected, and the errors are due to the use of RK4 for computing the initial condition and the use of the 2-stage collocation RK method. Finally, in the case of option (a), the errors are due to the use of both numerical integration and the 2-stage collocation method in the solution of local problems: again, the errors are expected to be of order 4.
Similar results are obtained for first and second-order methods: the order of convergence is two if option (a), (b), and (c) are used since the mid-point method is second-order accurate, and we consider the 1-stage collocation RK method. Now, we compare the results obtained with NWBMi, WBMi, CWBMi, CLWBMi, i = 1, 2, 3. The cell-averages are computed using option (b) for NWBMi and WBMi, option (c) for CWBMi, and option (d) for CLWBMi. L 1 errors with respect to the initial condition are computed.  Figure 1 shows the differences between the stationary and the numerical solutions obtained with NWBMi and CLWBMi, i = 1, 2, 3. Similar results are obtanied for WBMi and CWBMi, i = 1, 2, 3. Tables 2 and 3 show the errors corresponding to NWBMi, WBMi, CWBMi, CLWBMi, i = 1, 2, 3.  Notice that the errors for NWBM decrease with ∆x at the expected rate. WBM preserve the exact solution up to MP, as well as CWBM and CLWBM, provided that the tolerance used in the Newton's method for the former and Algorithm 5 for the latter is small enough. The computational costs are shown in Table 4. Nonetheless, since the goal of carrying out these experiments is to prove that the introduced well-balanced technique works, the implementation is not optimized. It can be seen that the well-balanced reconstruction given by Algorithm 2 increases the computational cost by a factor ranging from 1.5 to 7.5 (see Table 4). On the other hand, the use of Algorithm 3 increases the computational cost by a factor ranging from 1 to 1.8 if collocation RK method are used for solving the local problems and a factor ranging from 1.1 to 3.8 if control techniques with N p = 1 (see Reference [39]) are used. Therefore, the CLWBM schemes introduced in this work are more efficient than the CWBM schemes presented in Reference [39].  100   1  20  30  70  30  2  30  60  140  110  3  40  190  240  230   200   1  20  60  230  100  2  40  190  330  280  3  110  480  530  520   400   1  50  180  520  220  2  100  530  1150  720  3  350  1680  1980  1870   800   1  140  570  2020  870  2  270  2040  3580  2820  3 1080 5540 6600 5960

Problem 2: Burgers Equation with a Non-linear Source Term II
We now consider Burgers equation with another non-linear source term: The stationary solutions satisfy the ODE Note that no simple expression of the stationary solutions is known so (10) has to be computed numerically.

Test 2.1
We consider again the same computational domain, and the system is integrated up to T = 5. The CFL parameter is imposed to be 0.9. As initial condition, we consider the stationary solution that corresponds to the solution of the Cauchy problem (58) with u(−1) = 2, that is approximated by the Gauss-Legendre collocation method. The left boundary condition is set to u(−1, t) = 2, and open boundary conditions are considered upstream. Figure 2 shows the differences between the stationary and the numerical solutions computed with NWBMi, i = 1, 2, 3 (left) and CLWBMi, i = 1, 2, 3 (right), (the figures corresponding to CWBMi, i = 1, 2, 3 are similar).
The errors corresponding to NWBMi, CWBMi, and CLWBMi, i = 1, 2, 3, respectively, are shown in Tables 5 and 6. Computational costs are shown in Table 7. Note that the collocation strategy increases the computational cost from 7.3 to 14, whereas the method proposed in Reference [39] increases the computational cost by a factor ranging from 27 to 40. After the study of these two experiments, one can conclude that the results obtained with the collocation strategy introduced in this paper are better than the ones got by CWBMi, i = 1, 2, 3 as described in Reference [39], considering that the errors are similar for both procedures, but the first strategy is much more efficient. Consequently, we will only consider the collocation strategy in the following experiments.   Now, we consider the evolution of a small perturbation of the previous stationary solution. Thus, the initial solution is set to where u * (x) is the previous stationary solution. Figure 3 shows the difference between the stationary and the numerical solutions computed with NWBMi, i = 1, 2, 3 and CLWBMi, i = 1, 2, 3 at times t = 0.3, 5 (the results are similar for CWBMi, i = 1, 2, 3). In black, we show a reference solution computed on a fine mesh with 12,800 cells using a first order well-balanced method.    Observe that, as expected, the well-balanced methods are able to recover the stationary solution after the perturbation left the domain, while non-well balanced methods introduce an error that remains in the domain. This could be also observed in Table 8, where the errors with respect to the stationary solution at time t = 5 s are shown for the 100-cell mesh.

Problem 3: Shallow Water Equations
Let us consider now the 1D shallow water system which corresponds to the choices N = 2, with where x corresponds to to the axis of the channel, and t is the time. The unknown h(x, t) and q(x, t) are the thickness and discharge, respectively; the function H(x) is the depth function, and g is the gravity; u = q/h is the depth-averaged velocity, and c = gh.
The eigenvalues of the system are λ ± = u ± c, and the Froude number defined by characterizes the flow regime: subcritical (Fr < 1), critical (Fr = 1), or supercritical (Fr > 1). The stationary solutions U * satisfy the ODE system : which can be written as follows: The stationary solutions are given in implicit form by: In Reference [38], a family of high-order well-balanced numerical methods was presented, based in the Algorithm 2, in which local problems were solved on the basis of this implicit form. Here, Algorithm 3 with collocation RK methods to solve the local problems is used.
The numerical treatment of resonant situations discussed in Section 3.6 will be illustrated here in the particular case of the shallow water system. It can be easily checked that system (50), which, in our case, becomes when U * is a critical state, i.e., when u 2 = gh, reduces to: Therefore, the system has solutions only if H x (x) = 0: in this case, the solutions are In fact, it can be shown (see Reference [38], for instance) that a smooth stationary solution can only reach a critical state at a minimum point x c of the depth function H. Differentiating the second relation of (62) with respect to x, and using the relation H x (x c ) = 0, one obtains, at With this information in mind, let us compute the limit to determine the value of h x at x c . This is a 0/0 indeterminate limit, and L'Hôpital's rule can be applied. To do this, first, the limit is rewritten as follows: Some easy computation leads to therefore: The above expression shows that, if q = 0 and H xx (x c ) = 0, it is h x (x c ) = 0, thus justifying the assumption used before. Then, the chosen solution of (50) will be More precisely, when applying Algorithm 5 or Algorithm 4, a threshold > 0 is selected to detect critical states. If a system (50) has to be solved in which |Fr(U * ) − 1| < , at a pointx, then • Ifx is not close to a minimum point of H, it is assumed that is not possible to find a smooth stationary solution that solves the problem and the algorithm is stopped. • Otherwise, the selected solution of (50) is , if h is decreasing close tox, and the algorithm goes on.

Test 3.1: Transcritical Smooth Stationary Solution
The objective with this test is to check the correct behavior of the numerical methods of this work in presence of critical states. Following Reference [38], we consider here a smooth transcrital stationay solution for the 1D shallow water system on the interval [0, 3]. The integration time is T = 1, and the bathymetry is given by otherwise.
The initial condition corresponds to the solution of the Cauchy problem (61) with initial condition h(0) = 1.67750727 and q(0) = 2.5. It can be checked that it has a critical state at x c = 1.5. The shallow water system is integrated up to time T = 1, setting q downstream, and the water height upstream. The difference between the stationary and the numerical solutions computed with NWBMi, i = 1, 2, 3 and CLWBMi, i = 1, 2, 3 are shown in Figure 4. L 1 Errors are shown in Tables 9 and 10. Notice that Reference [38] obtained similar results. The computational costs are shown in Table 11.   100   1  20  30  2  30  130  3  70  320   200   1  30  60  2  80  370  3

Test 3.2: Perturbation of a Smooth Transcritical Stationary Solution
In this test, we consider the evolution of a small perturbation of the previous transcritical stationary solution. More precisely, a small perturbation of the water thickness of amplitude ∆h = 0.02 is set at the interval [1.1, 1.2], that is in the subcritical region. Therefore, this perturbation splits in two wawes, one moving downstream and the other upstream, as shown in Figures 5 and 6 at time t = 0.1 s. Figures 5 and 6 show the difference between the stationary and the numerical solutions computed with NWBMi, i = 1, 2, 3 and CLWBMi, i = 1, 2, 3 at times t = 0.1, 5 s. A reference solution is also computed using a first order well-balanced method with fine mesh (2000 cells).
Observe that, as in the previous examples, only the well-balanced methods are able to preserve the stationary solution. The non-well-balanced method generates strong perturbations that are bigger than the initial perturbation. Differences between the stationary and the numerical solutions in L 1 -norm at time t = 5 s for a 200-cell mesh are shown in Table 12. Finally, let us point out that the third order well-balanced method presents a small spurious oscillation in the location of the critical point, that disappears with time. This is not the case for the first and second order well-balanced methods.

Test 3.3: Perturbation of a Smooth Subcritical Stationary Solution
Let us consider evolution of a small perturbation of a smooth subcritical stationary solution. The initial condition U 0 (x) = (h 0 (x), q 0 (x)) T is given by: The depth function is given again by (69). Figure 7 shows the differences between the stationary solution and the numerical solutions obtained with NWBMi, i = 1, 2, 3 and CLWBMi, i = 1, 2, 3 at times t = 0.1, 5 for h (the graphs are similar for q). A reference solution has been computed with a first order well-balanced scheme on a fine mesh (3200 cells).
Notice that, since the perturbation amplitude is small enough, the well-balanced methods perform much better than the non-well-balanced methods, that generates spurious waves of amplitudes bigger than the original perturbation. In Table 13, the differences with respect to the stationary solution at time t = 5 s are shown for the 100-cell mesh. In order to check the efficiency of the methods, Figure 8 shows the errors in L 1 -norm with respect to the reference solution versus the CPU times in milliseconds for the NWBM2 and CLWBM2 at time t = 1 s, when the perturbation is still in the domain. Similar results have been obtained for the first and third order methods. Notice that the CLWBM are computationally more efficient. See, for instance, that, in order to obtain an error of approximately 2 × 10 −4 , the computational cost for the NWBM2 is increased by a factor of about 2.2 with respect to the CLWBM2.

Problem 4: Shallow Water Equations with Manning Friction
Let us consider the shallow water equations with Manning friction: In the previous, k is the Manning friction coefficient, and η is set to 7 3 . The stationary solutions satisfy the ODE system In what follows, we consider some experiments taken from Reference [47].

Test 4.1
Let us consider first a moving equilibria test with constant bathymetry. Therefore, the smooth stationary solutions satisfy As in Reference [47], we consider the interval the space interval [0, 1], discretized with 200 cells. Now, we consider a subcritical steady state with q(0) = −1 and h(0) a positive root of the non-linear function ξ(h), obtained by integrating (73), given by In this experiment, k = 1, which enforces the stiff character of the friction term. As boundary conditions, q is set downstream and h upstream, and the system is integrated up to t = 1. Due to the stiff character of the friction term, the explicit non-well-balanced methods turn out to be unstable for small values of h and/or big values of k, unless a very small time step is used. If an implicit discretization of the stiff source term is considered, NWBMi, i = 1, 2, 3 methods do not explode, yet there is not convergence to any stationary solution. The stability restriction is much less severe for the well-balanced schemes. For these reasons, only the numerical solutions obtained with well-balanced schemes are shown. Figure 9 shows the discrepancies between the stationary and the numerical solutions at t = 1 s with CLWBMi, i = 1, 2, 3; Table 14 shows the L 1 errors for a mesh with 200 cells.
As in Reference [47], we consider a mesh with 100 cells, and the model is run up to t = 9 s. As in the previous test case, only the well-balanced schemes will be considered. Figure 10 shows the discrepancies between the stationary and the numerical solutions for the water depth computed with CLWBMi, i = 1, 2, 3 at the time steps t = 0.06 s and t = 9 s (the figures are similar for q). In black, we plot a reference solution computed with a 1st-order well-balanced scheme with a 1600-cell mesh.  As expected, the stationary solution is preserved by the well-balanced schemes after the perturbation left the domain. Table 15 shows the L 1 errors with respect to the steady state at t = 9 s for 100-cell mesh. We end this subsection with two numerical experiments whose objective is show the efficiency of the schemes to preserve steady states that involve a varying bathymetry and friction. In this case, k takes the value of 0.01 and the space domain is again the interval [0, 1]. The bathymetry is given by The first experiment concerns the preservation of the supercritical stationary solution corresponding to (72) with initial conditions q(0) = 1 and h(0) = 0.3.
Following Reference [47], we take 100 discretization cells, and the numerical simulation of this experiment is run until t = 1 s. Figure 11 shows the numerical solutions at t = 1 s for NWBMi, i = 1, 2, 3 and CLWBMi, i = 1, 2, 3. Notice that NWBM1 is completely wrong: due to the numerical diffusion, the supercritical regime is lost at the right side of the domain, which causes a shock wave traveling to the left (of course, this behavior disappears when the mesh is fine enough). Therefore, in Figure 12, where the differences with the stationary solution are plotted, only the results for NWBMi, i = 2, 3 and for CLWBMi, i = 1, 2, 3 are shown. Table 16 shows the errors corresponding to the different methods for the 100-cell mesh. We also consider a reference solution computed with a well-balanced first order scheme on a 1600-cell mesh.    Following Reference [47], the very last experiment in this subsection focuses on a perturbation of the aforementioned steady state. The initial condition U 0 (x) = [h 0 (x), q 0 (x)] T is given by: where U * (x) = [h * (x), q * (x)] T is the stationary solution considered in the previous experiment. As in Reference [47], we use a 100-cell mesh for the numerical simulation on the domain [0, 1], and the computations are carried out until the final time t = 2 s. Observe that there are not big differences among the solutions during for small times. Nevertheless, for the NNWBM1, the supercritical regime is lost on the right of the domain, causing a shock traveling to the left, as in Test 4.3. Figure 13 shows then the differences between the numerical solutions and the stationary solution for h with NWBMi, i = 1, 2, 3 at time t = 0.015 and i = 2, 3 at time t = 2, and with CLWBMi, i = 1, 2, 3 at times t = 0.015 and t = 2 (the plots are similar for q). Again, in black, we plot a reference solution computed on a 1600-cell mesh with a first order well-balanced scheme.  Once more, well-balanced schemes are able to recover the stationary states after the perturbation left the domain. Table 17 shows the differences in L 1 -norm with respect to the stationary solution at time t = 2 s for a mesh with 100 cells. Finally, let us consider the compressible Euler equations for gas dynamics with gravity: Here, u is the velocity, p ≥ 0 the pressure, ρ ≥ 0 the density, q = ρu the momentum, E the total energy per unit volume, and H(x) is the gravitational potential. The internal energy e is defined by ρe = E − 1 2 ρu 2 , and pressure could be determined from the internal energy using the equation of state. In this work, we suppose an ideal gas; therefore, γ > 1 being the adiabatic constant. Here, it is set to γ = 1.5. Note that (76) is a particular case of (1) with N = 3, Supposing that the system is strictly hyperbolic and considering (77), the stationary solutions satisfy the ODE system: where c = γp/ρ is the sound speed, which is valid for stationary regular solution. Notice that, for regular solutions of the Euler equations, the entropy is constant along material lines. So that, for stationary solutions, us x = 0, where s denotes the entropy density. Therefore, if u does not vanish, it follows that these stationary solutions are isentropic. Test 5.1 In this test, we consider the space domain [−1, 1], and the system is integrated up to time t = 5 s. The gravity potential is H(x) = x.
As initial condition, a supersonic stationary solution corresponding to the solution of the Cauchy problem is computed with the Gauss-Legendre collocation method. Density, momentum, and energy is set downstream, and open free boundary conditions upstream. CFL parameter is set to 0.9. Figures 14-16 show the differences between the stationary and the numerical solutions for NWBMi and CLWBMi, i = 1, 2, 3 for ρ, q, and E, respectively; Table 18 shows the errors for for the 100-cell mesh with NWBMi and CLWBMi, i = 1, 2, 3 at time t = 5 s. Table 19 shows the computational cost for NWBMi, CWBMi, and CLWBMi, i = 1, 2, 3. Notice that CLWBM are significantly more efficient than CWBM. We have carried out similar experiments with a perturbation of this stationary solution as initial condition, and we obtained the expected results, i.e., only the well-balanced methods are able to exactly recover the steady state solution. For the sake of conciseness, we do not present the tests here.

Conclusions
Following Reference [36], we propose a family of well-balanced high-order numerical schemes that can be applied to general 1D balance laws. Note that the main difficulty of the method proposed in Reference [36] comes from the fact that local stationary solutions must be computed at every cell that may require the knowledge of the stationary states. When solving the ODE (2) for the stationary solution is not possible or too costly, the first step has to be computed numerically: what we propose here is to build both the stationary solutions and the local solvers for the first step on the basis of the collocation RK methods. The well-balanced property of these numerical methods is precisely stated, and we have also introduced a general technique that allows us to deal with resonant problems for any 1D systems of balance laws.
We have considered a good number of examples: Burgers equation, shallow water system with friction, and Euler equations with gravity to check the well-balanced property of the methods and to test their efficiency. In some of these test cases, the stationary solutions are known either in implicit or explicit form, while, in others, the only information comes from the ODE (2) that the stationary solutions satisfy: in the former case, we can compare the efficiency of the new implementation, while, in the latter, we can show the generality of the methods.
Although well-balanced schemes are more costly, they are more effective than standard non-well-balanced methods when they are applied to the propagation of small perturbation around an equilibria or when long-time integration is required. Furthermore, the tests show that the strategy based on the collocation RK methods introduced in this work is more efficient than the one presented in Reference [39].
The methods considered here can be applied in principle to any system of balance laws. For instance, in the context of Euler equations, in Reference [48], a first and second order numerical scheme that is well-balanced for particular 1D radial solutions for the Euler equations in cylindrical coordinates has been developed following the strategy proposed by two of the authors that was later reviewed in Reference [37]. It is possible to use the technique described here to extend the methods in Reference [48] to more general source terms and higher-order. This will be considered in a future work.
Futher developments will include applications to: • Balance laws (1) with non-regular source terms, that is systems for which H has jump discontinuities. • Multidimensional problems.

Conflicts of Interest:
The authors declare no conflict of interest.