# Well-balanced high-order finite difference methods for systems of balance laws

In this paper, high order well-balanced finite difference weighted essentially non-oscillatory methods to solve general systems of balance laws are presented. Two different families are introduced: while the methods in the first one preserve every stationary solution, those in the second family only preserve a particular stationary solution. The accuracy, well-balancedness, and conservation properties of the methods are discussed, as well as their application to systems with singular source terms. The strategy is applied to derive third and fifth order well-balanced methods for a linear scalar balance law, Burgers' equation with a nonlinear source term, and for the shallow water model.

## Authors

• 4 publications
• 3 publications
07/02/2020

### An order-adaptive compact approximation Taylor method for systems of conservation laws

We present a new family of high-order shock-capturing finite difference ...
03/23/2020

### High-order accurate entropy stable finite difference schemes for the shallow water magnetohydrodynamics

This paper develops the high-order accurate entropy stable (ES) finite d...
10/20/2020

### A universal centred high-order method based on implicit Taylor series expansion with fast second order evolution of spatial derivatives

In this paper, a centred universal high-order finite volume method for s...
07/18/2020

### An iterative scaling function procedure for solving scalar non-linear hyperbolic balance laws

The scaling of the exact solution of a hyperbolic balance law generates ...
05/18/2020

### Entropy conservation property and entropy stabilization of high-order continuous Galerkin approximations to scalar conservation laws

This paper addresses the design of linear and nonlinear stabilization pr...
05/04/2020

### High order discretely well-balanced finite volume methods for Euler equations with gravity – without any à priori information about the hydrostatic solution

We introduce novel high order well-balanced finite volume methods for th...
04/11/2020

### High-order well-balanced finite volume schemes for hydrodynamic equations with nonlocal free energy

We propose high-order well-balanced finite-volume schemes for a broad cl...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

We consider 1d systems of balance laws of the form

 Ut+F(U)x=S(U)Hx, (1)

where takes value in , is the flux function; ; and is a known function from (possibly the identity function ). The system is supposed to be hyperbolic, i.e. the Jacobian of the flux function is assumed to have

different real eigenvalues. PDE systems of this form appear in many fluid models in different contexts: shallow water models, multiphase flow models, gas dynamic, elastic wave equations, etc.

Systems of the form (1) have non trivial stationary solutions that satisfy the ODE system:

 F(U)x=S(U)Hx. (2)

The objective of well balanced schemes is to preserve exactly some of these steady state solutions. In the context of shallow water equations Bermúdez and Vázquez-Cendón introduced in [2] the condition called C-property: a scheme is said to satisfy this condition if it solves correctly the steady-state solutions corresponding to water at rest. Since then, the design of high-order well-balanced numerical methods for systems of balance laws has been a very active front of research: see, for instance, [1], [4], [5], [8], [9], [10], [11], [12], [13], [16] [17], [19], [22], [24], [28], [29], [30], [31], [34], [39] [40], [41], [42]

We focus here on finite difference methods. Previous well-balanced high-order finite difference types have been based in techniques like the one proposed in [20] consisting on a formal transformation of (1) into a conservative system through the definition of a ‘combined flux’ formed by the flux and a primitive of the source term (see, for instance, [4], [16]). In [39]-[40], finite difference WENO schemes where applied to problems in which can be written as a function of variables that are constants for the stationary solutions to be preserved.

The technique introduced here is different and it is based on a simple idea: if is the numerical approximation of the solution at the node at time , then if one can find a stationary solution satisfying

 U∗i(xi)=Ui,

one has

 S(Ui)Hx(xi)=S(U∗i(xi))Hx(xi)=F(U∗i(xi))x.

Therefore, the source term can be numerically computed with high order accuracy by applying the finite difference reconstruction operator to the point values . We will show that this simple strategy leads to high-order numerical methods that preserve either one particular stationary solution or all of them. Of course, the main difficulty comes from the computation of the stationary solution at every node at every time step.

This strategy has been inspired by the concept of well-balanced reconstruction introduced in [5] to develop well-balanced high-order finite volume numerical methods: see [11] for a recent follow-up. In that case, a stationary solution whose average is the numerical approximation at every cell has to be computed at every time step. From this point of view, the technique introduced here for finite difference methods is easier, since the local problems to be solved are standard Cauchy problems for ODE system (2). In all the numerical tests considered here, explicit or implicit expressions of the stationary solutions are available, which will allow us to find easily. For cases where this is not possible, a numerical method for ODE can be used to compute at the stencil of in the spirit of [6].

While the finite volume methods based on well-balanced reconstructions have the property of being conservative for any conservation law contained in (1) this is only the case for the finite difference methods that preserve one particular stationary solution to be described here but not for those that preserve every stationary solutions.

The case in which

has jump discontinuities will be also considered. At a discontinuity of

, a solution is expected to be discontinuous too and the source term cannot be defined within the distributional framework: it becomes a nonconservative product whose meaning has to be specified. There are different mathematical theories that allow one to give a sense to nonconservative products. In the theory developed in [15], nonconservative products are interpreted as Borel measures whose definition depends on the choice of a family of paths that, in principle, is arbitrary. As the Rankine-Hugoniot conditions, and thus the definition of weak solution, depend on the selected family of paths, its choice has to be consistent with the physics of the problem. Although for general nonconservative systems the adequate selection of paths may be difficult, in the case of systems of balance laws with singular source term there is a natural choice in which the paths are related to the stationary solutions of a regularized system: the interested reader is addressed to [11] for a detailed discussion. Although in general finite volume methods can be more easily adapted to deal with nonconservative products than finite difference schemes (see [32]), it will be shown that the numerical methods that preserve every stationary solution can be easily adapted to deal properly with singular source terms.

The organization of the article is as follows: finite difference high order methods based on a reconstruction operator are recalled in Section 2 where the particular example of WENO methods is highlighted. In Section 2 the case in which is continuous and a.e. differentiable is considered. First, the numerical methods that preserve any stationary solution are introduced, together with the proofs of their accuracy and their well-balanced property. Next, methods that preserve only one stationary solution are introduced and the conservation property is discussed. The implementation of well-balanced WENO methods is also discussed. Section 3 is devoted to the more difficult case in which is a.e. differentiable and piecewise continuous with isolated jump discontinuities. The definition of the nonconservative product is briefly discussed and the numerical methods introduced in Section 2 are adapted to this case. The numerical methods are applied to the linear transport equation with linear source term, Burgers’ equation with a nonlinear source term, and the shallow water equations. In particular, it will be shown that the methods introduced here preserve all the stationary solutions for the shallow water, including cases in which the bottom has a step. Finally, some conclusions are drawn and further developments are discussed.

## 2 Numerical methods

### 2.1 General case

We will consider here high-order finite difference numerical methods for (1) based on a reconstruction operator, i.e. an operator which, given the point values of a function in the nodes of a mesh

 vi=v(xi),∀i,

provides a numerical flux function

 ^vi+1/2=^v(vi−r,…,vi+s),∀i,

such that the flux difference approximates the derivative to -th order accuracy

 1Δx(^vi+1/2−^vi−1/2)=v′(xi)+O(Δxk),∀i.

when is smooth enough. Finite-difference ENO or WENO reconstructions are examples of such reconstruction operators: see [23], [35], [36]. Uniform meshes of constant step will be considered and the following notation will be also used for the intercells

 xi+1/2=xi+Δx2.

Once the reconstruction operator has been chosen, a possible semi-discrete numerical method for solving (1) can be written as follows:

 dUidt+1Δx(ˆFi+1/2−ˆFi−1/2)=S(Ui)Hx(xi), (3)

where are the flux reconstructions:

 ˆFi+1/2=ˆF(F(Ui−r),…,F(Ui+s)), (4)

The semidiscrete method (3) will be discretized in time by using TVD-RK methods: see [21].

### 2.2 WENO reconstructions

In the particular case of the WENO reconstruction of order , two flux reconstructions are computed using the values at the points :

 ˆFLi+1/2 = ˆFL(F(Ui−k),…,F(Ui+k)), (5) ˆFRi−1/2 = ˆFR(F(Ui−k),…,F(Ui+k)). (6)

Once these reconstructions have been computed, an upwind criterion can be chosen to define . For instance, for scalar problems

 ut+f(u)x=s(u)Hx, (7)

an approximate value of is chosen and the numerical flux is defined then as follows:

 ˆfi+1/2=a+i+1/2fLi+1/2+a−i+1/2fRi+1/2,

where

 a±=12(a±|a|).

For systems, a matrix with real different eigenvalues

 λi+1/2,1,…,λi+1/2,N

that approximates the Jacobian of the flux function has to be chosen and then the numerical flux can be defined by

 ˆFi+1/2=P+i+1/2ˆFLi+1/2+P−i+1/2ˆFRi+1/2, (8)

where

 P±i+1/2=Ki+1/2D±i+1/2K−1i+1/2. (9)

Here is the diagonal matrix whose coefficients are

 12(1±sign(λi+1/2,j)),j=1,…,N,

and

is a matrix whose columns are eigenvectors.

An alternative approach is to split the flux

 F(U)=F+(U)+F−(U)

in such a way that the eigenvalues of the Jacobian (resp. ) of (resp. ) are positive (resp. negative). Then, the reconstruction operator is applied to :

 ˆF+i+1/2 = ˆFL(F+(Ui−k),…,F+(Ui+k)), (10) ˆF−i−1/2 = ˆFR(F−(Ui−k),…,F−(Ui+k)), (11)

and finally,

 ˆFi+1/2=ˆF+i+1/2+ˆF−i+1/2. (12)

A standard choice is the Lax-Friedrichs flux-splitting:

 F±(U)=12(F(U)±αU),

where is the maximum of the absolute value of the eigenvalues of , this maximum being taken over either local (WENO-LLF) or global (WENO-LF): see [23], [36].

In both cases (the upwind or the splitting implementations) the values used to compute are those at the points

 xi−k,…,xi+k+1

so that, in the general notation, and .

## 3 Well-balanced high-order finite difference methods: H continuous

### 3.1 Definition of the method: general case

In order to tackle the difficulties gradually, let us suppose first that is continuous and a.e. differentiable. In this case, the stationary solutions of (1) are the solutions of the ODE system:

 d dxF(U)=S(U)Hx. (13)

Once the reconstruction operator has been chosen, the following semi-discrete numerical method is proposed:

 dUidt+1Δx(ˆFi,i+1/2−ˆFi,i−1/2)=0, (14)

where are computed as follows:

1. Look for the stationary solution such that:

 U∗i(xi)=Ui. (15)
2. Define

 Fj=F(Uj)−F(U∗i(xj)),j=i−1−r,…,i+s
3. Compute

 ˆFi,i+1/2 = ˆF(Fi−r,…,Fi+s), ˆFi,i−1/2 = ˆF(Fi−1−r,…,Fi−1+s).

The following result holds:

###### Proposition 1.

If the numerical method (3) is well-defined and the stationary solutions of (1) are smooth, then it has order of accuracy .

###### Proof.

Let be a smooth solution of (1). Given a time and an index , the reconstruction procedure is applied to to obtain , where represents the solution of (17) that satisfies

 U∗,ti(xi)=U(xi,t). (16)

One has:

 ∂tU(xi,t)+1Δx(ˆFi,i+1/2−ˆFi,i−1/2) =∂tU(xi,t)+∂xF(U)(xi,t)−∂xF(U∗,ti)(xi)+O(Δxk) =∂tU(xi,t)+∂xF(U)(xi,t)−S(U∗,ti(xi))∂xH(xi)+O(Δxk) =∂tU(xi,t)+∂xF(U)(xi,t)−S(U(xi,t))∂xH(xi)+O(Δxk) =O(Δxk),

where the facts that is a solution of (1) and a stationary solution satisfying (16) have been used. ∎

Observe that the method is well-defined if the first step of the reconstruction procedure can be always performed, i.e. if for every the Cauchy problem

 ⎧⎪⎨⎪⎩d dxF(U)=S(U)Hx,U(xi)=Ui, (17)

has a unique solution whose interval of definition contains the stencil . When (17) has no solution, or when it has a unique one but the stencil is not contained in its interval of definition, the cell values , to which the reconstruction operator is applied, cannot be the point values of a stationary solution of (1). In this case, there is no stationary solution to preserve and the numerical method (3) is applied to compute . When (17) has more than one solution, an additional criterion is needed to select one of them in the first stage of the reconstruction procedure.

###### Remark 1.

In the notation the index corresponds to the intercell and the index to the center of the cell where the initial condition of (17) is imposed. Therefore, in general

 ˆFi,i+1/2≠ˆFi+1,i+1/2 (18)

as one can expect due to the non conservative nature of the system of equations. Notice that two reconstructions have to be computed at every stencil: and .

###### Remark 2.

If the Jacobian of the flux function is regular, the ODE system can be written in normal form:

 ⎧⎨⎩dUdx=J(U)−1S(U)Hx,U(xi)=Ui, (19)

and it is expected to have a unique maximal solution . If one of the eigenvalues of vanishes in the stencil (the problem is said to be resonant in this case), (17) may have no solution or have more than one.

### 3.2 Well-balanced property

The numerical method (14) is well-balanced in the sense given by the following

###### Proposition 2.

Given a stationary solution of (1

), the vector of its point values

is an equilibrium of the ODE system given by the semi-discrete method (14).

###### Proof.

Observe that, in this case

 U∗,ti≡U∗,∀i,t

and thus

 F(U∗(xj))−F(U∗,ti(xj))=0,∀i,j=i−1−r,…,i+s.

The proof is trivial from this equality. ∎

### 3.3 Numerical method with WENO reconstructions

Let us discuss the implementation of the numerical method (14) in the particular case of WENO reconstructions with the upwind or the flux-splitting approach.

#### 3.3.1 Upwind approach

The implementation in this case is as follows: once the solution has been computed:

• Define

 Fj=F(Uj)−F(U∗i(xj)),j=i−k−1,…,i+k+1
• Compute

 ˆFLi,i+1/2 = ˆFL(Fi−k,…,Fi+k), ˆFRi,i−1/2 = ˆFR(Fi−k,…,Fi+k), ˆFRi,i+1/2 = ˆFR(Fi−k+1,…,Fi+k+1), ˆFLi,i−1/2 = ˆFL(Fi−k−1,…,Fi+k−1).
• Choose intermediate matrices .

• Define

 ˆFi,i+1/2 = P+i+1/2ˆFLi,i+1/2+P−i+1/2ˆFRi,i+1/2, ˆFi,i−1/2 = P+i−1/2ˆFLi,i−1/2+P−i−1/2ˆFRi,i−1/2,

where the projection matrices are given by (9).

#### 3.3.2 Flux-splitting approach

The implementation of WENO with splitting approach will be as follows: once the solution has been computed:

• Define

 F+j=F+(Uj)−F+(U∗i(xj)),j=i−k−1,…,i+k F−j=F−(Uj)−F−(U∗i(xj)),j=i−k,…,i+k+1
• Compute

 ˆF+i,i+1/2 = ˆF+L(F+i−k,…,F+i+k), ˆF−i,i−1/2 = ˆF−R(F−i−k,…,F−i+k), ˆF−i,i+1/2 = ˆF−R(F−i−k+1,…,F−i+k+1), ˆF+i,i−1/2 = ˆF+L(F+i−k−1,…,F+i+k−1).
• Define

 ˆFi,i±1/2=ˆF−i,i±1/2+ˆF+i,i±1/2.

In the particular case of the Lax-Friedrichs splitting, the reconstruction operators and will be applied to the values

 F(Uj)−F(U∗i(xj))±α(Uj−U∗i(xj))

in the corresponding stencil. Again is the local or global maximum eigenvalue of : although is not explicitly taken into account in the numerical viscosity, the numerical method has been shown to be stable under a CFL number of 1/2 in all the test cases considered here.

###### Remark 3.

Observe that, while for a conservative system 2 reconstructions are computed at every stencil , 4 reconstructions have to be computed at every stencil: 2 using , 1 using , 1 using .

### 3.4 Numerical methods that preserve only one stationary solution

If there is only one known stationary solution to preserve (which is the case if, for instance, the initial condition is a perturbation of a given steady state), the following semidiscrete numerical method can be used:

 dUidt+1Δx(ˆFi+1/2−ˆFi−1/2)=(S(Ui)−S(U∗(xi))Hx(xi), (20)

where are computed now as follows:

1. Apply the reconstruction operator to the point values

 Fj=F(Uj)−F(U∗(xj)),j=i−r,…,i+s

to obtain

 ˆFi+1/2=ˆF(Fi−r,…,Fi+s).
2. Apply the reconstruction operator to the point values

 Fj=F(Uj)−F(U∗(xj)),j=i−1−r,…,i−1+s

to obtain

 ˆFi−1/2=ˆF(Fi−r,…,Fi+s).

It can be easily checked that this method preserves the stationary solution and, if this solution is smooth enough, the order of the method is equal to that of the reconstruction operator. Observe that the first stage of the reconstruction procedure, that consisted of the resolution of a Cauchy problem for a ODE system, is avoided; therefore the computational cost is expected to be lower than that of the methods that preserve any stationary solution.

Another important difference between numerical methods (14) and (20) lies in their conservation properties: let us assume that there exists such that:

 S(U)=[0,…,0,SI+1(U),…,SN(U)]T,

i.e. the first equations of (1) are conservation laws. Then, (20) is conservative for these equations: in effect, the first components reduce to the conservative expressions:

 duj,idt+1Δx(ˆfj,i+1/2−ˆfj,i−1/2)=0,j=1,…,I,

where and represent respectively the th component of and . This is not the case for since inequality (18) is expected to hold at every component for (14) due to the fact that, in general, .

## 4 Well-balanced high-order finite difference methods: H discontinuous

Let us suppose now that is a.e. differentiable with finitely many isolated jump discontinuities. In this case, the definition of weak solutions (and, in particular, of stationary solutions) of (1) becomes more difficult: a solution is expected to be discontinuous at the discontinuities of and, in this case, the source term cannot be defined within the distributional framework. The source term becomes then a nonconservative product that can be defined in infinitely many different forms: see [15]. Like in [11], the following criterion is assumed here to define the weak solutions of (1):

• A pair of states can be the limits of an admissible weak solution of (1) to the left and to the right of a discontinuity point of , , if and only if there exists a solution of the ODE system

 d dσF(V)=S(V) (21)

satisfying

 V(H±)=U±. (22)

This admissibility criterion can be understood as a particular choice of paths within the theory developed by DalMasso, LeFloch, and Murat in [15] that allows one to give a sense to the nonconservative product as a Borel measure. In [11] it has been shown that the admissible jumps according to this criterion can be regularized through a smooth stationary solution in the following sense: is an admissible jump at a discontinuity point of if and only if for any given it is possible to find two smooth functions , defined in such that is monotone and:

 limx→x∗−ϵUϵ(x)=U−,limx→x∗+ϵUϵ(x)=U+,
 limx→x∗−ϵHϵ(x)=H(x∗−),limx→x∗+ϵHϵ(x)=H(x∗+),
 ∂xF(Uϵ)=S(Uϵ)∂xHϵ,%in(x∗−ϵ,x∗+ϵ).

Moreover, the jumps at the discontinuities of can be interpreted as stationary contact discontinuities of an extended system and the admissibility criterion (AC) is equivalent to assuming that Riemann invariants are preserved through these discontinuities (see [11] for details).

Let us remark that, if is a solution of (21), then is a stationary solution of (1): observe that (13) is trivially satisfied in smooth regions and (AC) is trivially satisfied in discontinuities. The well-balanced numerical method (14) can be easily adapted to this case by looking for stationary solutions of this form. More precisely, let us assume that the mesh has been designed so that the discontinuities of lie in intercells. Then, the numerical method (14) can be still applied with the only difference that now the first stage of the reconstruction step (15) is now as follows:

1. Look for the solution of (21) such that:

 V∗i(Hi)=Ui, (23)

where .

Every discontinuity point will be placed at an intercell . The reconstruction so defined leads to the following approximation of the Dirac delta produced by the discontinuity of :

 (F(V(H+(xI−1/2)))−F(V(H−(xI−1/2))))δ|x=x∗, (24)

where is a solution of (21): therefore, it is consistent with the assumed definition of the nonconservative products. Although consistency is not enough to guarantee the convergence to the right weak solution of nonconservative systems (see [33], [7]) the numerical tests in Section 5 show that the numerical methods capture the correct weak solutions.

In order to compare the behaviour of the different methods in the presence of a discontinuity of , for WENO methods (1) the Dirac delta issued from a discontinuity of will be approached by

 SI−1/2H+(xI−1/2)−H−(xI−1/2)Δxδ|x=x∗, (25)

where is some intermediate value of at the discontinuity. An upwind treatment of the singular source term is then used, so that the numerical method for the neighbor nodes writes as follows:

 dUI−1dt+1Δx(ˆFI−1/2−ˆFI−3/2)=S(UI−1)Hx(xI−1)+S−I−1/2, (26) dUIdt+1Δx(ˆFI+1/2−ˆFI−1/2)=S(UI)Hx(xI)+S+I−1/2,

where

 S±I−1/2=P±i+1/2SI−1/2H(xI)−H(xI−1)Δx.

Here, are the projection matrices given by (9). It will be seen in Section 5 that these methods do not converge to the assumed weak solutions. Moreover, the behaviour of the numerical methods at the discontinuity of depends both on and the chosen intermediate state.

In the case of the methods that only preserve one stationary solution , the source term in the neighbor nodes of the discontinuity will be computed as follows:

 dUI−1dt+1Δx(ˆFI−1/2−ˆFI−3/2)=(S(UI−1)−S(U∗(xI−1)))Hx(xI−1)+S−I−1/2−S∗−I−1/2, (27) dUIdt+1Δx(ˆFI+1/2−ˆFI−1/2)=(S(UI)−S(U∗(xI)))Hx(xI)+S+I−1/2−S∗+I−1/2,

with

 S∗±I−1/2=P±i+1/2S∗I−1/2H(xI)−H(xI−1)Δx.

Clearly, if for the right-hand sides vanish and the approximation of the Dirac mass is given again by (24) but, if it is not the case, a combination of (24) and (25) is used. Therefore, this method is only consistent with the definition of weak solution when the stationary solution is not perturbed at the neighbour nodes of the discontinuity, as it will be seen in Section 5

## 5 Numerical tests

In this section we apply the numerical methods introduced in Sections 2 and 3 to a number of test cases with continuous or discontinuous. Three families of methods based on WENO reconstructions of order will be compared:

• WENO: methods of the form (3).

• WBWENO: methods of the form (14) that preserve any stationary state.

• WB1WENO: methods of the form (20) that preserve only one given stationary state.

Nevertheless, in many test cases the results obtained with WBWENO and WB1WENO are indistinguishable: in those cases, only the results corresponding to WBWENO will be shown. In all cases, the global Lax-Friedrichs flux-splitting approach is used for WENO implementation and the third order TVD-RK3 method is applied for the time discretization: see [21]. The CFL parameter is set to 0.5.

### 5.1 A linear problem

We consider the linear scalar problem

 ut+ux=uHx. (28)

In this case, (21) reduces to

 dvdσ=v, (29)

whose solutions are

 V(σ)=Ceσ,C∈R.

The stationary solutions of (28) for any given are thus given by:

 u∗(x)=CeH(x),C∈R. (30)

The solution of (23) is thus

 v∗i(σ)=uieσ−Hi.

Therefore, well-balanced methods are based on the reconstructions of

 fj=uj−uie(Hj−Hi)Δx,j=i−r,…,i+s.

#### 5.1.1 Order test

Let us consider (28) with

 H(x)=x.

It can be easily checked that the solution of (28) with initial condition:

 u(x,0)=u0(x),x∈R

is given by

 u(x,t)=etu0(x−t),x∈R.

Let us consider the initial condition:

 u0(x)=⎧⎨⎩0if x<0,p(x)if 0≤x≤1,1otherwise, (31)

where is the 11th degree polynomial

 p(x)=x6(5∑k=0(−1)k(5+kk)(1−x)k)

such that

 p(0)=0,p(1)=1,pk(0)=pk(1)=0,k=1,…,5

see Figure 1.

We solve (28) with initial condition (31) with well-balanced and non well-balanced third order methods in the interval . Free boundary conditions based on the use of ghost cells are used at both extremes. Table 1 shows the -errors and the empirical order of convergence corresponding to WENO and WBWENO, . As can be seen, both methods are of the expected order and the errors corresponding to methods of the same order are almost identical. In order to capture the expected order, the smooth indicators of the WENO reconstruction have been set to 0 and has been chosen for the fifth order methods.

#### 5.1.2 A moving discontinuity linking two stationary solutions

Next, we consider (28) with again , and initial condition

 u0(x)={4exif x<0,exotherwise.

The solution consists of a discontinuity linking two steady states that travels at speed 1:

 u(x,t)={4exif x

Figure 2 shows the exact and the numerical solutions at time obtained with WBWENO3 and WENO3 (left) and a zoom of the differences of the numerical and the exact solution at the same time (right). It can be observed that the stationary states at both sides of the discontinuity are better captured with the well-balanced method. For fifth order methods the differences are lower, but still noticeable: see Figure 3. The results obtained with WBWENO, and the WB1WENO, that only preserve the stationary solution , are indistinguishable.

### 5.2 Burgers’ equation with source term

We consider next the scalar equation

 ut+f(u)x=s(u)Hx, (32)

with

 f(u)=12u2,s(u)=u2.

The stationary solutions are also given by (30). The reconstruction operator has to be applied in this case to

 fj=u2j2−u2ie2(Hj−Hi)Δx2,j=i−r,…,i+s.

#### 5.2.1 Preservation of a stationary solution with smooth H

In this test case we consider and the stationary solution

 u(x)=ex. (33)

Let us solve (32) taking this stationary solution as initial condition in the interval . As boundary conditions, the value of the stationary solution is imposed at the ghost cells. Figures 4 and 5 show the differences between the stationary solution and the numerical solutions obtained at time with WENO and WBWENO, using a 200-cell mesh: the well-balanced methods capture the stationary solution with machine accuracy. This is confirmed by Tables 2 and 3 that show the -errors and the empirical order of convergence corresponding to WENO and WBWENO, .

#### 5.2.2 Preservation of a stationary solution with oscillatory smooth H

Let us consider now (32) with a function that has an oscillatory behavior:

 H(x)=x+0.1sin(100x), (34)

(see Figure 6). We consider again the interval and we take as initial condition the stationary solution

 u(x)=ex+0.1sin(100x),

(see Figure 6) and, as boundary conditions, the value of this stationary solution is also imposed at the ghost cells.

We consider a 100-cell mesh, so that the period of the oscillations is close to . Figure 7 shows the numerical solutions at time corresponding to WBWENO, WENO, : while the well-balanced methods preserve the stationary solution with machine precision, the non well-balanced methods give a wrong numerical solution. Of course, they give more accurate solutions if the mesh is refined: see next paragraph, where the reference solution is computed using WENO3.

#### 5.2.3 Perturbation of a stationary solution with oscillatory smooth H

We consider again Burgers’ equation (32) with as in (34), and an initial condition that is the stationary solution approximated in the previous test with a small perturbation

 u0(x)=ex+0.1sin(100x)+0.1e−200(x+5)2,

(see Figure 8). If a mesh with 100 cells is used, the non well-balanced methods are unable to follow the evolution of the perturbation, since the numerical errors observed in the previous test are much larger than the perturbation. Let us see what happens when well-balanced methods are used: Figure 9 shows the numerical solutions obtained with WBWENO3 and WBWENO5 and a reference solution computed with WENO3 using a mesh of 5000 cells at time . As it can be seen, both methods are able to follow the evolution of the perturbation.

#### 5.2.4 Preservation of a stationary solution with piecewise continuous H

Let us consider now (32) with a piecewise continuous function :

 H(x)={0.1x if x≤0;0.9+xotherwise; (35)

(see Figure 10). We consider again the interval and we take as initial condition the stationary solution

 u(x)={e0.1x if x≤0;e0.9+xotherwise; (36)

(see Figure 10) and, as boundary conditions, the value of this stationary solution is also imposed at the ghost cells. WBWENO, preserve the stationary solution again with machine accuracy: see Table 4.

For WENO methods, the non-conservative product appearing at the source term is discretized by (26) with two different definitions of : a centered one

 sI−1/2=s(0.5(uI−1+uI)) (37)

or an upwind one

 sI−1/2=(1+sign(uI−1/2)2)s(uI−1)+(1−sign(uI−1/2)2)s(uI). (38)

Here, is the index of the intercell is located and is the arithmetic mean of and .

In Figure 11 we compare the exact solution with the numerical solutions obtained at time obtained with WBWENO3 using a 300-cell mesh (its graph and the one of the exact solution are identical at the scale of the figure) and with WENO3 or WENO5 using different implementations:

• WENO3-UPW1: WENO3 with upwind implementation and (38);

• WENO3-UPW2: WENO3 with upwind implementation and (37);

• WENO3-LF: WENO3 with LF implementation and (38);

• WENO5-LF: WENO5 with LF implementation and (38).

As it can be observed, the results of the WENO methods depends on the chosen implementation, on the numerical definition of the source term, and on the order. Moreover, these differences remain as tends to 0. This is in good agreement with the difficulties of convergence of finite difference methods to nonconservative systems: see [7].

#### 5.2.5 Perturbation of a stationary solution with piecewise continuous H

We consider again (32) with (35) and an initial condition that is the stationary solution approximated in the previous test with a small perturbation

 ~u0(x)=