Well-balanced algorithms for relativistic fluids on a Schwarzschild background

We design well-balanced numerical algorithms with first-, second-, and third-order of accuracy for the spherically symmetry evolution of a compressible fluid on a Schwarzschild black hole. We treat both the relativistic Burgers-Schwarzschild model and the relativistic Euler-Schwarzschild model. Our schemes follow the finite volume methodology and preserve the stationary solutions and allow us to investigate the global asymptotic behavior of such flows and reach definite conclusions about the behavior of the mass density and velocity field.


Introduction
We design a finite volume algorithm for computing relativistic fluid flows and apply it to study the global dynamics of a flow evolving around a (spherically symmetric) Schwarzschild black hole. We consider the Euler equations for inviscid and compressible fluids and, throughout, we assume that the flow enjoys spherical symmetry. Our purpose in this paper is to design shock-capturing schemes that are high-order accurate and well-balanced in the sense that they preserve the spatially homogeneous solutions. We build upon earlier investigations on this problem by LeFloch et al [6,10,11] and extend to the present problem the well-balanced methodology in Castro and Parés [4], in order to properly take the Schwarzschild curved geometry into account. For earlier work on well-balanced schemes we also refer to [2,14,15] and, concerning the design of geometry-preserving schemes, we refer for instance to [1,6,7,8,13,17] and the references therein.
We will treat first a simplified model, namely the relativistic Burgers-Schwarzschild model [11]: Here the unknown is a scalar function which represents the normalized velocity of the flow. This equation clearly takes the form v t + F (v, r) r = S(v, r),

3)
M > 0 is a coefficient representing the mass of the black hole. We will next consider the relativistic Euler-Schwarzschild model in which the unknowns are the fluid density ρ and the normalized velocity v(t, r) ∈ (−1, 1), defined for all r > 2M , while k ∈ (−1, 1) denotes the (constant) speed of sound. (The limiting values v = ±1 will be reached at the boundary r = 2M only.) We can write this system in the form V t + F (V, r) r = S(V, r), (1.5) where where This is a strictly hyperbolic system and the eigenvalues of the Jacobian of the flux function are (cf. for instance [11]) (1.9) As usual, a state (ρ, v), by definition, is • sonic if one of the eigenvalues vanishes, i.e. if |v| = k.
• supersonic if both eigenvalues have the same sign, i.e. if |v| > k.
• subsonic if the eigenvalues have different signs, i.e. if |v| < k.
We will need to distinguish between these regimes n order to design a robust scheme. An outline of this paper is as follows. We develop here high-order, well-balanced schemes for the relativistic Burgers equation (1.2) and next for the Euler equations (1.5) posed on a Schwarzschild background. The methodology is described first in Section 2, and the actual design of the schemes is presented in Section 3 (Burgers equation) and Section 5 (Euler equations). The efficiency, accuracy, and robustness of proposed algorithms are then investigated in Section Section 4 (Burgers equation) and Section 6 (Euler equations). We demonstrate that the global dynamics can be accurately determined by the proposed algorithms and we reach some definite conjectures in Sections 4.6 and 6.5, respectively.
In particular, by relying on extensive tests we demonstrate that the proposed schemes are numerically wellbalanced and we discuss the importance of this property in order to have reliable results. We study the late time behavior of the Burgers or Euler equations and discuss the roles of initial condition imposed at the boundary. We also describe how steady shocks behave under small or large perturbations. Further details are found in the corresponding sections.
We consider semi-discrete finite volume numerical methods of the form S(P t i (r), r) dr , (2.3) and the following notation is used: • I i = r i− 1 2 , r i+1/2 are the computational cells, whose length ∆r is supposed to be constant for simplicity. • V i (t) is the approximation of the average of the exact solution at the ith cell at time t, that is, • P t i (r) is the approximation of the solution at the ith cell given by a high-order reconstruction operator from the sequence of cell averages {V i (t)}: (2.5) where S i denotes the set of indices of the cells belonging to the stencil of the ith cell.
It can be then easily shown that, if the reconstruction operator is well-balanced for a continuous stationary solution V * of (2.2) then the numerical method is also well-balanced for V * . Given a stationary solution V * of (2.2), we use the following terminology.
• The numerical method (2.3) is said to be well-balanced for V * if the vector of cell averages of V * is an equilibrium of the ODE system (2.3).
• The reconstruction operator is said to be well-balanced for V * if P i (r) = V * (r), r ∈ [r i− 1 2 , r i+1/2 ], (2.7) where P i is the approximation of V * obtained by applying the reconstruction operator to the vector of cell averages of V * .
The following strategy to design a well-balanced reconstruction operator P i on the basis of a standard operator Q i was introduced in [2]. We are given a family of cell values {V i }, at every cell 1.
Seek the stationary solution V * i (x) such that 1 ∆r Apply the reconstruction operator to the cell values {W j } j∈Si given by in order to obtain Q i (r) = Q i (r; {W j } j∈Si ).

Define finally
(2.10) It can be then easily shown that the reconstruction operator P i in (2.10) is well-balanced for every stationary solution provided that the reconstruction operator Q i is exact for the zero function. Moreover, if Q i is conservative, then P i is conservative, in the sense that 1 ∆r (2.11) and P i has the same accuracy as Q i if the stationary solutions are sufficiently regular. The well-balanced property of the method can be lost if a quadrature formula is used to compute the integral appearing at the right-hand side of (2.3). In order to circumvent this difficulty, in [4] it is proposed to rewrite the methods as follows: where V t, * i is the stationary solution found in (2.8) at the ith cell at time t. In this equivalent form, a quadrature formula can be applied to the integral without losing the well-balanced property , and this leads to a numerical method of the form: (2.13) Here, α 0 , . . . , α q , r i,0 , . . . , r i,q represent the weights and the nodes of the chosen quadrature formula, whose order of accuracy must be greater or equal to the one of the reconstruction operator. If the quadrature formula is used to compute the averages of the initial condition as well: the reconstruction procedure has to be modified to preserve the well-balanced property: Steps 1 and 2 have to be replaced by the following: (2.14) 2. Apply the reconstruction operator to the cell values {W j } j∈Si given by Finally, for a first-order method we consider the trivial reconstruction operator, given the cell averages {V i } of a function V , and we consider the piecewise constant approximation of V It can be easily checked that the numerical method then reduces to Moreover, if the initial averages are computed with the midpoint rule, namely (2.17) then at the first step of the reconstruction procedure, the problem (2.14) reduces to finding the stationary solution satisfying For a second-order method the mid-point rule can also be selected, so that the problem (2.14) reduces again to (2.18).
3 Burgers-Schwarzschild model: designing the numerical algorithm

Preliminaries
For the Burgers-Schwarzschild equation (1.1), the steady state solutions are of the form In Figure 3.1.1 we plot the steady solutions for several values of K 2 . The domain of definition of these stationary solutions is It can be easily checked that, given a pair (K, r * ) such that r * ∈ D K , the discontinuous function defined in D K by is an entropy weak stationary solution of the Burgers-Schwarzschild model.

First-order method
If the midpoint rule is used to compute the initial averages, at the first step of the reconstruction procedure one has to search for K 2 i such that There is always a unique solution given by In order to apply the numerical method (2.12), this stationary solution has to be computed at r i± 1 2 and this is only possible if r i+1/2 ∈ D Ki , that is, provided If this condition is satisfied, then the numerical method (2.16) can be used. If this condition is not satisfied, then the standard trivial reconstruction is considered, i.e. Q i (r) = v i . The numerical method writes then as follows: Summing up, the expression of the semi-discrete first-order method reads The forward Euler method is used for the time discretization. We emphasize that, if (3.7) is not satisfied, then v i cannot be the point value of a stationary solution defined in the computational domain, so that the use of the standard reconstruction does not destroy the well-balanced property of the method, since in this case there is no stationary solution to preserve.

Second-order method
Let us suppose again that the midpoint rule is used to compute the cell averages and that the minmod reconstruction operator is considered: see [16]. The stationary solution v * i selected at the first stage of the reconstruction procedure is again (3.6) withK i given by (3.5). In order to compute the reconstructions, this stationary solution has to be computed at the points r i−1 , r i− 1 2 , r i+1/2 , r i+1 so that the following condition has to be satisfied r i+1 ∈ D Ki , i.e.
If this condition is satisfied, the following step of the reconstruction procedure consists in computing the fluctuations: (3.12) Then the reconstruction is defined as (3.14) Once the well-balanced reconstruction operator has been computed, the form (2.12) of the numerical method is considered and the midpoint rule is used again to approximate the integral. Observe however that, in this case: Therefore, the expression (2.13) reduces again to (2.16) with F i+1/2 = F(v t,− i+1/2 , v t,+ i+1/2 , r i+1/2 ). If (3.11) is not satisfied, then the standard MUSCL reconstruction is applied, namely The expression of the numerical method is given again by (3.9)-(3.10) with the difference that the second-order reconstructions are used now to compute the numerical fluxes. The TVDRK2 method [9] is used in order to discretize the equations in time.
Observe that, according to the well-balanced reconstruction procedure described in the previous section, the fluctuations to be reconstructed should be in this case but in (3.12) the sign of v i has been replaced by that of v j . This modification allows one to preserve not only the continuous stationary solutions solution but also the discontinuous stationary solutions of the family (3.3).

Third-order method
In order to design a third-order method the CWENO reconstruction of order 3 (see [5], [12]) will be considered and the two-point Gauss quadrature will be used to compute the initial averages and the integrals coming from the source term: Therefore, at the first step of the reconstruction procedure one has to look for K 2 i such that: If we define the function it can be easily verified that g is a positive decreasing function defined in the interval [0, Therefore there are two possibilities: has to be computed in the points {r i−1,0 , r i−1,1 , r i+1,0 , r i+1,1 } in the reconstruction procedure. Therefore, these points have to be in the interval of definition of v * i , and this happens if r i+1,1 ∈ D Ki . Therefore, the well-balanced reconstruction can be computed only if the following condition is satisfied: If this condition is satisfied, the fluctuations can be then computed: and the well-balanced reconstruction is given by where Q represents the CWENO approximation. If (3.17) is not satisfied, the standard CWENO reconstruction is applied: The expression of the semi-discrete method is finally (3.9) where the numerical fluxes are computed at the reconstructed states and (3.18) The TVDRK3 method of order 3 will be used for the time discretization.

Preserving the exact averages of the stationary solutions
The three methods presented above can be modified to preserve the exact averages of the stationary solutions instead of its approximation computed with a quadrature formula. To do this, the problem to be solved at the first stage of the well-balanced reconstruction procedure is the following one: find K 2 i such that: 1 ∆r If we define the function it can be easily checked that it is a decreasing function defined in [0, K 2 e,i ] where K 2 e,i = 1 − 2M/r i+1/2 −1 and g(0) = 1. Therefore, (3.19) has a unique solution K 2 i if |v i | ≤ g(K e,i ). (3.20) The explicit expression of g can be obtained: is a primitive of 1 − x 1 − 2M r . Therefore g(K e,i ) can be explicitly computed. The well-balanced reconstruction can be computed if (3.20) has a solution K i and the cells of the stencil S i are contained in D Ki . the domain of definition of the stationary solution corresponding to the constant K i solving (3.19) what gives restrictions similar to the ones found in the previous paragraphs. If these restrictions are not satisfied, the standard reconstruction is applied. The expression of the numerical methods is the same that the ones above.

.1 Preliminaries
In this section several numerical tests are applied to check the performance of the well-balanced numerical methods introduced in the previous section. We consider the spatial interval [2M, L] with M = 1 and L = 4, a 256-point uniform mesh and the CFL number is set to 0.5. At x = 2M we impose F − 1 2 = 0 as boundary condition because 1 − 2M r = 0. At x = L we will use a transmissive boundary condition using ghosth-cells if the initial condition is not a stationary solution; otherwise, the stationary solution is imposed in the ghost-cells. The following numerical flux will be used: where q(·; v L , v R ) is the self-similar solution of the Riemann problem for the standard Burgers equation with the initial condition In order to check the relevance of the well-balanced property, these methods will be compared with those based on the same numerical fluxes and the standard first-, second-, or third-order reconstructions.

Stationary solutions
We consider the initial condition corresponding to the positive stationary solution with K = 1/2. Table 1 shows the error in L 1 norm between the initial condition and the numerical solution at time t = 50.  Let us consider now as initial condition the negative stationary condition corresponding to K = 1/2: Table 2 shows the error in L 1 norm between the initial condition and the numerical solution at time t = 50.    Table 2 we need more time to see the differences between the well-balanced and non-well-balanced schemes of order 1,2 and 3 but the errors are again much smaller with the well-balanced schemes for this test.
Let us consider finally the discontinuous initial condition that is a stationary entropy weak solution of the family (3.3). Table 3 shows the error in L 1 norm between the initial condition and the numerical solution at time t = 50.

Moving shocks connecting two steady profiles
We consider now the initial condition The corresponding solution consists of a right-moving shock connecting two branches of stationary solutions. Figure 4.3.1 shows the numerical solutions obtained with the first-, second-, and third-order well-balanced methods and a reference solution computed with the first-order standard method using a mesh of 10000 cells. As it can be seen, the well-balanced methods capture correctly the shock with a resolution that increases with the order as expected. Similar conclusions can be drawn for the left-moving shock linking two branches of stationary solutions that   generates from the initial condition:

Perturbation of a steady shock solution
In this test case we consider the initial condition: where v 0 is the steady shock solution given by (4.3) and The first-, second-, and third-order well-balanced method have been applied to this problem. In Figure 4.4.1 it can be observed that, after the wave generated by the initial perturbation leaves the computational domain, the stationary solution (4.3) is not recovered: a different stationary solution of the family (3.3) is obtained whose shock is placed at a different location. Observe that all the three methods capture the same stationary solution.  Similar conclusions can be drawn if a perturbation at the right side of the shock is superposed to the stationary solution v 0 given by (4.3):ṽ see Figure 4.4.2. In this case we have used a 2000-point uniform mesh because the displacement of the shock is smaller in this case and more points in the mesh are needed in order to see that the steady shock is not recovered. Now we consider an initial condition of the form (4.6) with a perturbation p L such that:  In order to study the relation between the amplitude of the perturbation and the distance between the initial and the final stationary shocks, we consider the initial condition: v 0 (r) = v 0 (r) + p L (r) + p R (r), (4.12) where v 0 is the steady shock solution given by (4.3) and (p L (r) + p R (r))dr = 0. In particular we take p L (r) as in (4.7) and p R (r) as in (4.9). In Figure 4.4.5 it can be observed that, after the wave generated by the initial perturbation leaves the computational domain, the stationary solution (4.3) is not recovered: a different stationary solution of the family (3.3) with the shock is placed at a different location. This is a natural result because as we saw before the right perturbation creates a lower displacement than the left perturbation. Here we have used again a 2000-point uniform mesh. In order to study the relationship between the amplitude of the perturbation and the distance between the initial and the final shock locations, we consider the family of initial conditions: where v 0 is given again by (4.3) and with α > 0.                   Table 4 we can conclude the following:

Conclusion 1. If a perturbation δ v is added to a steady shock solution of the form
then the solution reaches at finite time another steady shock solution of the form: δ v = 0, then r 1 = r 0 , i.e. the initial stationary solution is recovered.

If
δ v , then r 1 = r 0 and a different stationary solution is obtained.
3. If δ v = 0, then r 1 = r 0 and a different stationary solution is obtained. In this case the distance between r 0 and r 1 depends linearly on the amplitude of the perturbation: see Table 4 and If v 0 (r) = 1 for r ∈ [2M, 2M + ), with > 0 and v 0 (L) = a, with a < 0, then in finite time the solution has the form of a right-moving shock that links the stationary solution v ≡ 1 and the negative steady solution v * such that v * (2M ) = −1 and v * 0 (L) = a, that is: (c) If v 0 (2M ) < 1 and v 0 (L) ≥ 0, then in finite time the solution coincides with the negative steady solution such that v * (2M ) = −1 and v * 0 (L) = 0, that is: (d) If v 0 (2M ) < 1 and v 0 (L) = a, with a < 0, then in finite time the solution coincides with the negative stationary solution v * such that v * (L) = a, that is: 2. For the unbounded domain [2M, ∞) the following conclusions can be drawn by passing to the limit when L → ∞: (b) If v 0 (r) = 1 for r ∈ [2M, 2M + ), with > 0 and lim r→∞ v 0 (r) = a, with a < 0, then in finite time t 0 the solution has the form of a rigt-moving shock that links the stationary solution v ≡ 1 and the negative stationary solution v * such that v * (2M ) = −1 and lim r→∞ v * 0 (r) = a, that is: (c) If v 0 (2M ) < 1 and lim r→∞ v 0 (r) ≥ 0, then the solution converges as t → ∞ to the negative stationary solution v * such that v * (2M ) = −1 and lim r→∞ v * 0 (r) = 0, that is: (d) If v 0 (2M ) < 1 and lim r→∞ v 0 (r) = a, with a < 0, then the solution converges as t → ∞ to the negative stationary solution v * such that v * (2M ) = −1 and lim r→∞ v * 0 (r) = a, that is: 5 Euler-Schwarzschild model: designing the numerical algorithm

Preliminaries
In the case of the Euler-Schwarzschild equations (1.4), the stationary solutions are implicitly given by the equations: The pair (v, ρ) of a stationary solution satisfies the ODE system: see [11].
that comes from solving dv dr = 0. In Figure 5.1.1 the red stationary solutions are those such that at r = r c the take the value v = ±k.
Given two constants C 1 and C 2 , in order to compute the stationary solution given by (5.1) in a point r = a, the following non linear system has to be solved: It is enough thus to solve, if it is possible, the nonlinear equation It can be easily checked that g satisfies: Moreover, g is strictly monotone in [−1, −k), (−k, k), and (k, 1]. As a consequence we can conclude: • If |K a | > g(k) equation (5.6) has no solution, i.e. a stationary solution given by C 1 and C 2 cannot be defined at r = a.
• If |K a | = g(k) then equation (5.6) has only one solution (k if K a > 0, −k if K a < 0). Therefore, (5.5) has only one solution that is a sonic state.
• Otherwise, (5.6) has two possible solutions. Therefore there are two states (ρ, v) that solve (5.5), one supersonic and one subsonic.
For the sake of clarity, together with the representation for the states, we will use V can be easily computed from V and, given V , V is also easily computed by (1.8) that comes from solving a second-degree equation.

First-order method
If the midpoint rule is used to compute the initial averages, given a family of cell values V i , in the first step of the well-balanced reconstruction procedure one has to find, if it is possible, a stationary solution V * i defined in Obviously such a stationary solution would correspond to the choice of constants: (5.10) According to the discussion above, the corresponding stationary solution is defined in r i±1/2 provided that: When |K i± 1 2 | < g(k) there are two possible values for V * i (r i± 1 2 ), one subsonic and one supersonic. Therefore, a criterion is needed to select one or the other. The following criterion wil be used here: • if V i is not sonic, then the state whose regime (sub or supersonic) is the same as V i is selected for V * i (r i± 1 2 ). • if V i is sonic, then the state whose regime is the same as V i+1 is selected for V * i (r i+ 1 2 ) and the state whose regime is the same as V i−1 is selected for V * i (r i− 1 2 ). Observe that this criterion aims to preserve the regime of the given cell values.
If condition (5.11) is satisfied, then the numerical method (2.16) is used. Otherwise the standard trivial reconstruction is considered.
The expression of the semi-discrete first-order method is then as follows: otherwise. (5.14) The forward Euler method is used again for the time discretization.

Second-order method
Let us use again the midpoint rule to compute cell averages and the minmod reconstruction operator. The stationary solution sought at the first stage of the well-balanced reconstruction procedure is again characterized by the constants (5.9)-(5.10). This time, this stationary solution has to be computed at the points r i−1 , r i− 1 2 , r i+1/2 , r i+1 so that the following condition has to be satisfied: where K i± 1 2 are given by (5.12) and If this condition is satisfied, the reconstruction is defined as follows: where the minmod function is applied component by component and Observe that the conserved variables V are used in the reconstruction procedure. If (5.15) is not satisfied, then the standard MUSCL reconstruction is applied: The expression of the numerical method is given again by (5.13)-(5.14) with the difference that the second-order reconstructions are used now to compute the numerical fluxes. The TVDRK2 method is used now to discretize the equations in time.

Third-order method
Although it will not be implemented in the present paper, let us briefly describe the first step of a third-order well-balanced reconstruction procedure based on the two-point Gauss quadrature i order to compute averages:. It consists on finding C 1 and C 2 such that where V * (r; C 1 , C 2 ) represents the value at r of a stationary solution characterized by the constants C 1 and C 2 .

.1 Preliminaries
In this section several tests are considered to check the performance of the first-and second-order well-balanced numerical methods for Euler-Schwarzschild model introduced in the previous section. We consider the spatial interval [2M, L] with M = 1 and L = 10, a 500-point uniform mesh, k = 0.3 and the CFL number is set to 0.5 again. At x = 2M we impose F − 1 2 = 0 as boundary condition because 1 − 2M r = 0. At x = L we will use a transmissive boundary condition in the case we are not in a stationary solution or we will expand the stationary solution if we are in one.
In order to test the dependency of the results on the numerical method, two different first-order numerical fluxes are considered: the Lax-Friedrichs numerical flux and a HLL-like numerical flux in PVM form (see [3]): where λ 1 and λ 2 are the eigenvalues of some intermediate matrix J i+1/2 of the form where v m is some intermediate value between v n i and v n i+1 . Given two states V L and V R , in order to choose an adequate intermediate value v m , we look for v such that the following Roe-type property is satisfied: where i.e. the factor (1 − 2M/r) is neglected for simplicity. Due to the form of the matrix, it is enough to find v such that This equality is equivalent to the second-order degree equation for v: Since the discriminant is always positive, there are always two real solutions: and it can be proven that: Finally, in the case α = 0 and V R = V L , we take v m = − γ β and in the case Once v m has been chosen, the expression of λ j , j = 1, 2 in (6.3) is as follows: Since for a 2-systems HLL and Roe methods are equivalent and the intermediate value chosen to compute the wave speeds satisfies a Roe-type property, this numerical flux will be called Roe-type numerical flux in what follows.
The numerical methods will be compared with those based on the same numerical flux and the standard firstand second-order reconstructions.

Stationary solutions
We take as initial condition the positive supersonic stationary solution satisfying ρ * (10) = 1, v * (10) = 0.6. (6.6) Table 5 shows the error in L 1 norm between the numerical solution at time t = 50 for the well-balanced and non-well-balanced methods using the Roe-type numerical flux.  Table 5: Well-balanced versus non-well-balanced schemes: L 1 errors at time t = 50 for the Euler-Schwarzschild model with the initial condition (6.6).
Let us consider now as initial condition the negative supersonic stationary solution V * that satifies ρ * (10) = 1, v * (10) = −0.8. (6.7) Table 6 shows the error in L 1 norm between the numerical solution at time t = 50. Figures 6.2.3, 6.2.4 show the difference between the numerical results given by well-balanced and non-well-balanced methods using the Roe-type numerical flux. Again the numerical solutions obtained with non-well-balanced methods depart from the initial steady state. Let us consider finally the initial condition V 0 (r) = V * − (r) if r ≤ 6, V * + (r) otherwise, (6.8) where V * − (r) is the supersonic stationary solution such that ρ * − (6) = 4, v * − (6) = 0.6 (6.9) and V * + (r) is the subsonic one such that . (6.10) V 0 is an entropy weak stationary solution of the system: see [11]. Table 7 shows the error in L 1 norm between the numerical solution at time t = 50 and Figures 6.2.5, 6.2.6 show the comparison of the numerical results obtained with well-balanced and non-well-balanced methods at selected times.
The numerical results of this section put on evidence, as for the Burgers-Schwarzschild system, the relevance of using well-balanced methods for the Euler-Schwarzschild model.

Perturbation of a regular stationary solution
In this test we consider the initial condition V 0 (r) = V * (r) + δ(r), (6.11) where V * is the supersonic stationary solution such that ρ * (10) = 1, v * (10) = 0.9 (6.12) and It can be observed in Figure 6.3.1 that the stationary solution V * is recovered once the perturbation has left the domain. In this Figure, the numerical results obtained with the first-and second-order well-balanced methods are compared with a reference solution computed using the first-order well-balanced method with a 5000-point mesh (the Roe-type numerical flux is used again). As expected, the second-order method is less diffusive.

Perturbation of a steady shock solution
We consider the initial condition: V 0 (r) = V * (r) + δ L (r), (6.14) where Let us consider now two different initial conditions: on the one hand where V * (r) is again the discontinuous stationay solution given by (6.8)-(6.10) and On the other hand, V 0,2 (r) = V * 2 (r) + δ R (r), (6.18) where δ R is given again by (6.16) and V * 2 (r) is the steady shock of the form (6.8) satifying ρ * − (6) = 5, v * − (6) = 0.6. (6.19) Observe that the definition of v is identical for both stationary solutions but ρ is different.   Table 7: Well-balanced versus non-well-balanced schemes: L 1 errors at time t = 50 for the Burgers-Schwarzschild model with the initial condition (6.8) After the passage of the perturbation, the shock starts moving leftward and, in both cases, the numerical solution converges to a smooth transonic stationary solution of the form: where r c is given by (5.4); V * − (r) and V * + (r) are respectively a subsonic and a supersonic stationary solution satisfying v * ± (r c ) = −k: see Figure 6.4.4. Nevertheless, the limits in time of the approximations of ρ are different: see Figure 6.4.5. Observe that, in the Euler-Schwarzschild model (5.1), there are infinitely many stationary solutions with the same function v and different ρ.
In order to study the relationship between the amplitude of the perturbation and the distance between the initial and the final shock locations, we consider the family of initial conditions: where V * is the steady shock solution given by (6.8)-(6.10) and with α > 0. In this case we will also use the Roe-type numerical flux and a 2000-point uniform mesh.  as we did for the Burgers-Schwarzschild model.  Let us finally consider a family of initial conditions that generate leftward movement of the initial steady shock : where V * is again the steady shock solution given by (6.8)-(6.10) and     1. If the perturbation moves the steady shock to the right, then a different stationary solution of the form with r 0 = r 1 , is obtained, and the distance between r 0 and r 1 seems to depend exponentially on the amplitude of the perturbation: see Table 8 and Figure 6.4.8.
2. If the perturbation moves the steady shock to the left, then a steady shock solution of the form (6.20) is obtained.   O1_WB_Roe_2000-mesh O1_WB_Roe_4000-mesh O1_WB_Lax_2000-mesh O1_WB_Lax_4000-mesh O1_WB_Roe_2000-mesh O1_WB_Roe_4000-mesh O1_WB_Lax_2000-mesh O1_WB_Lax_4000-mesh O1_WB_Roe_2000-mesh O1_WB_Roe_4000-mesh O1_WB_Lax_2000-mesh O1_WB_Lax_4000-mesh 8 (6.14): comparison between the first-order well-balanced method with different meshes using the Roe-type and the Lax numerical fluxes at selected times for the variable v. O1_WB_Roe_2000-mesh O1_WB_Roe_4000-mesh O1_WB_Lax_2000-mesh O1_WB_Lax_4000-mesh