Numerical method for hyperbolic conservation laws via forward backward SDEs

It is well known that for solutions of semi-linear parabolic PDEs, there are equivalent probabilistic interpretations, which yields the so called nonlinear Feymman-Kac formula. By adopting such formula, we consider in this work a novel numerical approach for solutions of hyperbolic conservation laws. Our numerical method consists in efficiently computing the viscosity solutions of conservation laws. However, instead of solving the viscosity problem directly (which is difficult), we find its equivalent probabilistic solution by adopting the Feymman-Kac formula, which relies on solving the equivalent forward backward stochastic differential equations. It is noticed that such framework possesses the following advantages: (i) the viscosity parameter can be chosen sufficiently small (say $10^{-10}$); (ii) the computational procedure on each discretized time level can be \textit{completely parallel}; (iii) the traditional CFL condition is dramatically weakened; (iv) one does not need to handle the transition layers and discertizations of derivatives. Thus, high accuracy viscosity solutions can be efficiently found. Several numerical examples are given to demonstrate the effectiveness of the proposed numerical method.


Introduction
Many problems in physics and engineering are modeled by hyperbolic systems of conservation laws. For example, the shallow water equations of hydrology, the Euler equations for inviscid and compressible flow, and the Magnetohydrodynamic equations of plasma physics, see [7,13]. The general initial value problems for hyperbolic conservation laws yield where U : R d → R m is the vector of unknowns and F j : R m → R m is the flux vector for the jth direction with m being a positive integer. It is well known that solutions of (1.1) develop discontinuities in finite time even when the initial data is smooth. This holds true even for the scalar case and solutions to (1.1) are sought in the weak sense. Furthermore, weak solutions are augmented by additional admissibility criteria or entropy conditions [7] in order to ensure uniqueness.
There has been a large amount of work on the classical numerical methods for solving hyperbolic conservation laws [13,19], starting from the Godunov scheme [14], to the Lax-Wendroff scheme [20] and the upwind scheme [31], to the total variation bounded high-order essentially non-oscillatory (ENO) scheme [15][16][17] and the weighted essentially non-oscillatory method [32], to some recent approaches such us the moving mesh method (cf. [34]) and the (local) discontinuous Galerkin methods (cf. [37,43]), to name a few.
Another classical approach for conservation laws is the viscosity method, in which one introduces a diffusion term into the original equation and obtains a semi-linear parabolic equation which admits a unique smooth solution, precisely, (1.2) Here, ǫ is a small parameter. It is shown that for certain conservation laws [4,6,9] the limit of the smooth solution u ǫ (ǫ → 0) is the viscosity solution of the original problem (1.1). The viscosity method plays an important role in theoretical analysis for hyperbolic conservation laws. However, the direct simulation for (1.2) is tough due to sharp transition layers, and thus, it is hard to compute the viscosity solution with high accuracy. For details about the viscosity method, please refer to [4,6,9,36] and references therein.
In this paper, we will essentially extend the Crank-Nicolson type scheme introduced in [39] for BSDEs to solve the corresponding weakly coupled nonlinear FBSDEs. It is noted that in such framework (compared to an direct solver for (1.2)), one does not need to handle the transition layers and discertizations of the derivatives. Furthermore, such approach possesses the following features: (i) the viscosity parameter can be chosen sufficiently small (say 10 10 ); (ii) the computational procedure on each discretized time level can be completely parallel ; (iii) the traditional CFL condition is dramatically weakened. To demonstrate the effectiveness of our method, several model problems such as the linear advection equation, the Burgers' equation, the Buckley-Leverett equation, and the Euler system are tested.
We summarize here the new features of our discussions: (i) we provide an alternative way to compute the viscosity solution of conservation laws; (ii) we extend the Crank-Nicolson type scheme introduced in [39] for BSDEs to solve the resulting weakly coupled nonlinear FBSDEs. Despite the large amount of the literatures, such problems have not been studied.
The rest of this paper is organized as follows. In Section 2, by introducing the nonlinear Feymman-Kac formula, we describe the equivalent probabilistic solution of the viscosity solution. Such probabilistic solution can be obtained by solving the corresponding FBSDEs. The numerical approach for FBSDEs will be discussed in Section 3. Several numerical tests are given in Section 4 to demonstrate the effectiveness of our numerical method. We finally given some conclusions in Section 5.

Conservation laws, viscosity solution, and the nonlinear Feymman-Kac formula
Without loss of generality, we consider here the following one dimensional hyperbolic system of conservation law where u = (u 1 , ..., u m ) ⊤ is the vector of unknowns. As we stated, the main purpose of this work is to find the following viscosity solution of (2.1) where ǫ is a small positive number. For scalar conservation laws, u ǫ converges to an entropy solution when ǫ goes to zero (see e.g. [4,6,9]). Furthermore, consider the following Hamilton-Jacobi problem with strictly convex hamiltonian H and compactly supported (or periodic) C 2 initial data u ǫ | t=0 = u 0 (x), there holds [21]: Theorem 2.1. Let {u ǫ } ǫ be the family of approximate viscosity solutions of (2.3), and u be the solution of the original problem without the viscosity term. For any fixed T , there exists C(T ) such that the following error estimate holds The viscosity method plays an important role in the theoretical analysis for conservation laws. However, solving equation (2.2) with high accuracy is a difficult work due to the sharp transition layers and the nonlinearity. Therefore, instead of solving (2.2) directly, we consider in this work solving its equivalent probabilistic solution. To do this, we first review some related results of the nonlinear Feymman-Kac formula.

Nonlinear Feymman-Kac formula
The nonlinear Feymman-Kac formula reveals that for solutions of certain semilinear parabolic PDEs, there are equivalent probabilistic interpretations, which yield the representations of the solutions of FBSDEs. To begin, let us first introduce the following coupled FBSDEs Under certain regularity conditions on the data b, σ, f and ϕ, the adapted solution (Y t , Z t ) of (2.4) can be represented in the following form [28,30]: where u(t, x) is the smooth solution of the following parabolic partial differential equation is the solution of the FBSDEs (2.4), then u defined by (2.5) is the viscosity solution of (2.6).
The representations in (2.5) are well known as the nonlinear Feynman-Kac formula.

Probability representation of the viscosity solution
Consider the viscosity equation in the general form of By the Feymman-Kac formula (2.5), the solution u ǫ of the viscosity equation (2.8) can be represented as where (Y τ , Z τ ) is the solution of the following FBSDE . Our aim of this paper is to numerically solve the associated FBSDEs (2.9) to obtain u ǫ (τ, X τ ) = Y τ (Note that we do not need to solve Z τ ). We remark that there has been a few works on using FBSDEs to study PDE problems (cf. [2,11] and references therein). However, the PDE models considered and the corresponding numerical methods derived in these works are totally different with the results presented in this paper.

Discretization of the FBSDEs
In this section, we discuss the discretizations of the FBSDEs (2.9). We remark that one can find related information on such numerical approach in e.g. [39][40][41].
We introduce the following equally distributed time partition, We define the density of the grids in D n h by 2. the number of elements in D n h,x is finite and uniformly bounded, that is, there exists a positive integer N e , such that, #D n h,x ≤ N e . We call D n h,x the neighbor grid set in D n h at x.

The reference equations
By the BSDE in (2.9), we have for 0 ≤ n ≤ N − 1 : Then, taking the conditional mathematical expectation E x tn [·] on both sides of (3.2) yields The integrand E x tn [f (s, X s , Y s )] on the right-hand side of (3.3) is a deterministic smooth function of time s under F tn , and we may use some numerical integration methods to approximate the integral in (3.3). In particular, we approximate this integral by the trapezoidal rule, namely, where R n y is the truncation error.
For the FSDE of X t , we introduce the following scheme where R n x is the error introduced by the scheme.
Remark 3.1. The equations of (3.5) and (3.6) are our reference equations of the FBSDEs, and our numerical method will be designed based on them.

The semi-discrete scheme
Based on the reference equations (3.5) and (3.6), let Y n and X n be the numerical approximations of Y tn and X tn respectively, we propose a Crank-Nicolson type numerical scheme to solve Y n : (3.7) From Scheme 1, it is clear that Y n is a function of X n . Note that Scheme 1 is nonlinear for Y n . To solve it, some iteration methods are needed. In this paper, we use the following iteration method: , and set a tolerance parameter δ.
2. For l=0, 1, ..., solve Y l+1 n . Note that we have used Y n+1 (X n ) as the initial gauss for Y 0 n , and it is noted that in such a framework, for a reasonable δ (say δ = 10 −5 ), about 2-3 iterations are enough for convergence.

Remark 3.2.
Although the main purpose of this paper is to provide an alternative way to compute the viscosity solutions of conservation laws, our algorithms Scheme 2 is, to the best of our knowledge, new for solving FBSDEs .

The fully discrete scheme
Note that Scheme 1 is a semi-discrete scheme. For any fixed t n , to solve Y (t n , X tn ), the space partition of X tn is needed, and the conditional expectations should be approximated at each grid point x ∈ D n h . We denote such an approximation by E x,h tn [·|X tn = x] ≈ E tn [·|X tn = x]. Generally speaking, after taking the conditional expectation E tn [·|X tn = x], X n+1 in Scheme 1 does not belong to D n+1 h for x ∈ D n h . Thus, interpolation methods are needed to approximate the value of Y n+1 at X n+1 by the values of Y n+1 on D n+1 h . Here, we will adopt a local interpolation operator I n h,x such that I n h,x g is the interpolation value of the function g at space point x ∈ R q by using the values of g only on D n+1 h,x . Note that any interpolation methods can be used here, however, care should be made if one wants to guarantee the stability and accuracy. Now, we are ready to introduce the following fully discrete scheme: Notice that the iteration method like Scheme 2 is still needed for the fully discrete scheme. Remark 3.3. It is worth to note that in each time level t n , the computations at the spatial grids are independent with each other, and thus, the computational procedure on each discretized time level can be completely parallel. This means the total computational time for each time level is decided by the computational time for one single grid point, provided that enough processors are available, and this is obviously promising for real world simulations.
Remark 3.4. For smooth data b, f and ϕ, it is known that the solution Y t = u(t, X t ) is smooth with respect to t and X t , and furthermore, the conditional expectation E x tn [f (s, X s , Y s )] is a smooth function of s. Thus, the local truncation error R n y has the local estimate which implies that the numerical method can be a high-order numerical scheme.

Gauss-Hermite quadrature rule for E x tn [·]
As mentioned in the last section, numerical approximations for the conditional expectations in Scheme 3 are needed. In this work, we will use the Gauss-Hermite integral to approximate the condition expectations. The Gauss-Hermite quadrature rule is an extension of Gaussian quadrature method for approximating the value of integrals of are the corresponding weights [1]: The truncation error R(g, L) of the Gauss-Hermite quadrature formula (3.10) is where η is a real number in R. The Gauss-Hermite quadrature formula (3.10) is exact for polynomial functions g of degree less than 2L − 1.
For a d-dimensional function g(x), x ∈ R d , the Gauss-Hermite quadrature formula becomes x 2 j , and .
It is well known that, for a standard d-dimensional normal random variable N (0, 1), it holds that Then by (3.12), we deduce where R GH E,L (g) is the truncation error of the Gauss-Hermite quadrature rule for g.
Recall that, in Scheme 3, the conditional expectation E xi is the approximation of E x tn [·], and I n h,x Y n+1 is the interpolation approximation of Y n+1 . By the nonlinear Feynman-Kac formula, Y n+1 has the following explicit representation.  (3.14). Other kinds of conditional expectations can be approximated in a similar way.
Remark 3.5. Note that the Gaussian Hermite approximation for the conditional expectations admits some local properties, namely, only some local informations of Y n+1 (function values around X n ) contribute much for the integral. Then, the scaling idea for the Gaussian Hermite approximation [35] can be used, and in such a way, 6-7 points are enough to guarantee a good approximation.

Numerical experiments
In this section, we provide some numerical examples to demonstrate the effectiveness of our numerical method for solving hyperbolic conservation laws. The numerical examples chosen here will involve all kinds of flow features, that is, shocks, contacts, expansion fans, sonic points, and smooth regions. As we discussed in the last section, we will use an equally spaced grid in space with step ∆x and a uniform time step ∆t. We denote by N x the number of grid points used in the spatial discretizations. In all our computations, 8 points will be used to evaluate conditional expectations.

The advection equation
We first consider the following scalar advection equation where φ(x) is a periodic function with period 2. Notice that the exact solution is In what follows, the computational domain will be fixed to [−1, 1].
To show the convergence rates of our numerical method, we consider here a smooth initial condition φ(x) = − sin(πx). The exact solution is smooth without sonic points.
In Table 1 we list the discrete L 1 -errors and convergence rates (r c ) at t = 2 with different parameters ǫ, ∆t/∆x. The notations Maxvalue and Minvalue represent the maximum and minimum values of the numerical solutions. The results clearly show that our scheme (FBSDEs) is a high-order scheme, and it conserves the Maxvalue principle. Moreover, the convergence rate is maintained even when a large CFL condition is considered (∆t/∆x = 1.6 the second row). Also, the convergence rate is independent of the viscosity parameter ǫ (the twice and third row, ǫ = 10 −6 , 10 −8 ).

The Burgers' equation
We next consider the following scaler nonlinear Burgers' equation The initial function φ(x) is given by (4.3) 9.9991e-1 -9.9991e-1 81 1.5330e-5 9.9999e-1 -9.9999e-1 161 1.9229e-6 1.0000 -1.0000 One can derive the following exact piecewise linear solution for (4.2) and (4.3): Notice that the shock and the expansion fan interact for t > 2/3, which complicates the solution. The exact solution admits at x = −1/3 a sonic expansion fan and at x = 1/3 a steady shock. We   Fig.1 (Left). Also, we provide the numerical results with the parameters ǫ = 10 −10 , ∆t ∆x = 2.0 and 2.8, and N x = 121 in Fig.1 (Right). Both cases show that our scheme behaves very well for such problem with the CFL condition numbers either bigger or smaller than 1.

The Buckley-Leverett equation
The Buckley-Leverett equation is where the non-convex flux function f (u) is given by

The Euler System
The Euler system of gas dynamics can be written as with u = (ρ, v, p) ⊤ , where ρ is the density, v is the velocity and p is the pressure. The matrix A is given by with a = γp ρ being the speed of sound. For details of the problem see e.g. [16]. The Euler system can be rewritten into the following equivalent characteristic  Test case 1. A well-known test problem is the shock tube problem proposed by Sod [33]. The initial conditions are given by We simulate the problem up to time t = 1.5 with ∆x = 0.05 and ∆t/∆x = 1/2.4. Test case 2. Another frequently used test problem is the Lax test case proposed by Lax [18], for which the initial states are (ρ L , v L , p L ) = (0.445, 0.698, 3.528), (ρ R , v R , p R ) = (0.5, 0, 0.571).
Compared to the Sod's problem, the contact discontinuity and the shock of such test are stronger. The problem is computed up to t = 1.5 with space step size ∆x = 0.05 and ∆t/∆x = 0.2. The numerical solutions for the density, the velocity and the presser are given in Fig.3, Fig.4 and Fig.5, respectively. The left plots are for the test case 1 while the right plots are for the test case 2. The plots show that our scheme admits no spurious overshoots or oscillations, and captures the expansion fan, the contact discontinuity and the shock very well.
We take ∆x = ∆y = 0.05 and solve the problem up to t = 0.5 with ∆t/∆x = 1/4 and ∆t/∆x = 1/3. The numerical solutions are shown in Fig.6 (a), and we also show the solutions that cut at some x and y, respectively in Fig.6 (b) and Fig.6 (c). It is noted that our scheme solves the problem very well.

Conclusions
In this work, we provide an alternative way to solve the viscosity solution of hyperbolic conservation laws. The method consists in solving the equivalent probability solution (FBSDEs) of the viscosity solution. Numerical approach for the resulting weakly coupled FBSDEs is discussed. It is noticed that in such framework, the viscosity parameter can be chosen sufficiently small (say 10 −10 ), and the computational procedure on each discretized time level can be completely parallel, furthermore, the traditional CFL condition is dramatically weakened. Several numerical tests are given to show the effectiveness of the numerical method. Although the original motivation of this work is to provide an alternative to compute the viscosity solutions of hyperbolic conservation laws, however, we would like to emphasize that our numerical method for solving FBSDEs is, to the best of our knowledge, novel.
We remark that for multi-dimensional problems, we have introduced the tensorized grids in space, which is computational inefficient. In our future studies, we will update this part into the sparse grid setting, which is efficient for high dimensional approaches. Numerical comparisons between our scheme and other methods (e.g., TVD, WENO, etc.) will also be part of our future work.
We address the issue here and want to open up the possibility of designing efficient stochastic approaches (vis FBSDEs) for more general semi-linear parabolic PDEs, such as general convection-dominated diffusion problems, the Schrödinger equations, etc.