Log In Sign Up

An adaptive sparse grid local discontinuous Galerkin method for Hamilton-Jacobi equations in high dimensions

We are interested in numerically solving the Hamilton-Jacobi (HJ) equations, which arise in optimal control and many other applications. Oftentimes, such equations are posed in high dimensions, and this poses great numerical challenges. This work proposes a class of adaptive sparse grid (also called adaptive multiresolution) local discontinuous Galerkin (DG) methods for solving Hamilton-Jacobi equations in high dimensions. By using the sparse grid techniques, we can treat moderately high dimensional cases. Adaptivity is incorporated to capture kinks and other local structures of the solutions. Two classes of multiwavelets are used to achieve multiresolution, which are the orthonormal Alpert's multiwavelets and the interpolatory multiwavelets. Numerical tests in up to four dimensions are provided to validate the performance of the method.


Adaptive sparse grid discontinuous Galerkin method: review and software implementation

This paper reviews the adaptive sparse grid discontinuous Galerkin (aSG-...

A sparse grid discontinuous Galerkin method for the high-dimensional Helmholtz equation with variable coefficients

The simulation of high-dimensional problems with manageable computationa...

Solving Stochastic Optimization by Newton-type methods with Dimension-Adaptive Sparse Grid Quadrature

Stochastic optimisation problems minimise expectations of random cost fu...

A Rotating-Grid Upwind Fast Sweeping Scheme for a Class of Hamilton-Jacobi Equations

We present a fast sweeping method for a class of Hamilton-Jacobi equatio...

Sparse grid time-discontinuous Galerkin method with streamline diffusion for transport equations

High-dimensional transport equations frequently occur in science and eng...

p-adaptive algorithms in Discontinuous Galerkin solutions to the time-domain Maxwell's equations

The Discontinuous Galerkin time-domain method is well suited for adaptiv...

Sparse Grid Discretizations based on a Discontinuous Galerkin Method

We examine and extend Sparse Grids as a discretization method for partia...

1 Introduction

In this paper, we consider the Hamilton-Jacobi (HJ) equation


on the bounded domain in arbitrary dimension, subject to initial condition and appropriate boundary conditions. The HJ equation has diverse applications in science and engineering, such as optimal control, seismic waves, crystal growth, robotic navigation, image processing, calculus of variations, among others [35]. In this paper, we develop a class of adaptive sparse grid (also called adaptive multiresolution) discontinuous Galerkin (DG) methods for approximating the viscosity solution of (1.1). The concept of viscosity solution was developed by Crandal and Lions in [12, 11] to single out the physically relevant weak solution. Under certain assumptions, the viscosity solution can be interpreted by the Hopf formula [18], and the numerical approximation to the viscosity solution is of practical interest [42]. It is well known that the viscosity solution of the HJ equation is only Lipschitz continuous and may develop discontinuous derivatives in finite time regardless of smoothness of the initial condition. Various numerical methods for solving (1.1) have been developed in the literature [42], such as the monotone methods [13, 1, 32], the essentially non-oscillatory (ENO) methods [38, 39], the weighted ENO (WENO) and the Hermite WENO (HWENO) methods [28, 47, 40] among many others. In this work, we choose to use the DG discretization due to their distinguished advantages in handling geometry, boundary conditions and accommodating adaptivity, which are highly desirable for efficiently solving the HJ equation. Several DG schemes have been proposed in the literature [24, 33, 5, 34, 46, 22, 6, 30]. Here, we use the local DG (LDG) method developed by Yan and Osher [46], which has provable property for the piecewise constant case and is easy to implement, although we remark that the extensions to other DG formulations are also possible.

Beyond the need to capture the viscosity solution, another major numerical challenge for HJ equation is that it is often posed in high dimensions, and any standard numerical discretization becomes inefficient due to the curse of dimensionality. Recent years have seen a surge of interests in designing numerical solutions of HJ equations in high dimensions. Various approaches have been proposed, including those using sparse grid

[3, 19, 29], model order reduction [31]

, tensor decomposition

[17], Hopf formula [16, 7, 8]

and machine learning

[23, 37, 14, 15], to name a few. Some of the work above is feasible for HJ equations in hundreds of dimensions for some special cases, and continued efforts to develop efficient numerical solvers for high-dimensional HJ PDEs constitutes a vibrant research area due to their wide applications in control and differential games.

In this paper, we take the sparse grid approach [4], which has been used in [3, 19, 29] for computations in moderately high dimension. The scheme we proposed relies on multiresolution analysis (MRA) [36] and is designed to be high order accurate. In a line of research, we have developed a family of adaptive sparse grid (or adaptive multiresolution) DG methods for linear transport equations with application to kinetic equations [21], hyperbolic conservation laws [25], and wave equations [26]. By incorporating MRA and the sparse grid ideas, our methods are able to efficiently capture small-scale structures, and more importantly, work very well in high dimensions. In particular, in [25, 26], we use two classes of multiwavelets to achieve MRA. The tensor-product Alpert’s multiwavelets are used as the DG function space, following the approach developed in [45, 20, 21] for linear equations. Besides, the interpolatory multiwavelets for MRA quadrature [44] are used for computations of nonlinear terms. Numerical experiments for benchmark tests in up to four dimension verify the efficiency and efficacy of the method in capturing the viscosity solution of the HJ equations.

The rest of the paper is organized as follows. In Section 2, we review the fundamentals of Alpert’s and interpolatory multiwavelets. In Section 3, we present the LDG method for the HJ equation with MRA. Some theoretical results and implementation details are discussed. Section 4 contains numerical examples. In Section 5, we include the conclusion of this paper.

2 Multiresolution Analysis and Multiwavelets

In this section, we first review the fundamentals of MRA of DG approximation spaces and the associated multiwavelets. Two classes of multiwavelets, namely the orthonormal Alpert’s multiwavelets [2] and the interpolatory multiwavelets [44], are considered. We also introduce a set of key notations used throughout the paper.

2.1 Alpert’s multiwavelets

We start with the construction of Alpert’s multiwavelets [2], which have been employed to develop a class of sparse grid DG methods for solving PDEs in high dimensions [45, 20]. For a unit sized interval , we define a set of nested grids , for which the -th level grid consists of uniform cells

Denote The piecewise polynomial space of degree at most on grid for is denoted by


Observing the nested structure

we can define the multiwavelet subspace , as the orthogonal complement of in with respect to the inner product on , i.e.,

By letting , we obtain a hierarchical decomposition , i.e., MRA of space . A set of orthonormal basis can be defined on as follows. When , the basis , are the normalized shifted Legendre polynomials in . When , the Alpert’s orthonormal multiwavelets [2] are employed as the bases and denoted by

We then follow a tensor-product approach to construct the hierarchical finite element space in multi-dimensional space. Denote as the mesh level in a multivariate sense, where denotes the set of nonnegative integers, we can define the tensor-product mesh grid and the corresponding mesh size Based on the grid , we denote as an elementary cell, and

as the tensor-product piecewise polynomial space, where represents the collection of polynomials of degree up to in each dimension on cell . If we use equal mesh refinement of size in each coordinate direction, the grid and space will be denoted by and , respectively. Based on a tensor-product construction, the multi-dimensional increment space can be defined as

The basis functions in multi-dimensions are defined as


for , and .

Using the notation of

and the same component-wise arithmetic operations and relations as defined in [45], we reach the decomposition


On the other hand, a standard choice of sparse grid space [45, 20] is


We skip the details about the property of the space, but refer the readers to [45, 20]. In Section 3, we will describe the adaptive scheme which adapts a subspace of according to the numerical solution, hence offering more flexibility and efficiency.

2.2 Interpolatory multiwavelets

Alpert’s multiwavelets described in Section 2.1 are associated with the projection operator. The interpolatory multiwavelets introduced in [44]

are constructed based on interpolation operators and also essential for efficient computation of integrals in the DG formulation, especially in high dimensions. In this work, only Lagrange interpolation is considered, while we note that Hermite interpolation can also be used but its implementation is more involved. The details are provided below.

We first define the set of interpolation points on the interval at zeroth mesh level by . Here, the number of points in is . We defer the discussion of the relations between and to Section 3.2.

The interpolation points at mesh level , can be obtained correspondingly as

We require the points to be nested, i.e.


This can be achieved by requiring .

Given the nodes, we define the basis functions on the zeroth level grid as Lagrange interpolation polynomials of degree which satisfy the property:

for . It is easy to see that With the basis function at mesh level zero, we can define the basis functions at mesh level :

which form a complete basis set for

We now introduce the hierarchical representations and the interpolatory multiwavelets. Define and for , then we have the decomposition

Denote the points in by . Then the points in for can be represented by

For notational convenience, we let The increment function space for is introduced as a function space that satisfies


and is defined through the multiwavelets that satisfies

for . Then is given by

The multi-dimensional construction follows similar lines as in Section 2.1. We let


while the sparse grid approximation space is

Note that the constructions by Alpert’s multiwavelets and the interpolatory multiwavelets deduce the same sparse grid space because of the same nested structure. Finally, the interpolation operator in multidimension is defined as :

where the multi-dimensional basis functions are defined in the same approach as (2.2) by tensor products:


For the sparse grid space or any adaptively chosen subspace of the interpolation operator, which is denoted by in later sections, can be defined accordingly, by taking only multiwavelet basis functions that belong to that space. The fast algorithms which transform point values at interpolation points to hierarchical coefficients are given in [44]. The detailed formulas of the interpolation points and the associated interpolatory multiwavelets used in this work, we refer readers to [25, 26].

3 Adaptive multiresolution LDG scheme

In this section, we present the adaptive multiresolution LDG method for simulating the HJ equation (1.1). We start with reviewing the LDG formulation by Yan and Osher in [46]. Then, by incorporating MRA and multiwavelets introduced in the previous section, we define our scheme.

3.1 LDG formulation

We consider periodic boundary conditions for simplicity, while the method can be adapted to other non-periodic boundary conditions. For illustrative purposes, we first introduce a set of shorthand notation. Denote by the union of the boundaries for all the elements in the partition . The jump and average of are defined as

where n is the unit normal. ‘

’ and ‘+’ represent that the directions of the vector point to interior and exterior at

, respectively. Note that , and we let denote the -th component of .

The key idea in [46] is to employ the standard LDG methodology, see e.g. [10], to reconstruct the first derivatives of , i.e., , . In particular, the LDG method computes two piecewise polynomials and , both approximating but using opposite one-sided numerical fluxes; that is, given we seek , , in such that for all


where the numerical fluxes are defined as

Note that and carry the information of from opposite directions. Hence, when the solution is smooth, and are almost identical, while if the solution involves nonsmooth corners, then and can be very different.

Then the semi-discrete scheme for solving (1.1) is defined as follows: seek such that, for all ,


where denotes a monotone numerical Hamiltonian that approximates , and , , are given in (3.1). In the simulations, we employ the following global Lax-Friedrichs Hamiltonian

where and

with the maximum being taken over the whole domain.

Depending on the choice of space , we obtain several LDG methods for (1.1) with distinct properties. If we recover the full grid LDG scheme in [46] on tensor-product meshes. If then we obtain the sparse grid LDG method. If is chosen adaptively, we have the adaptive sparse grid scheme. Noteworthy, besides the LDG formulation, we can employ other DG formulations as well, such as the direct DG method [6] and the indirect DG methods [24, 22]. The LDG formulation used is comparatively simpler to implement under the MRA framework.

If the Hamiltonian is linear, then the HJ equation (1.1) degenerates to a transport equation with constant coefficients, and the formulation (3.2) together with (3.1) is nothing but a standard upwind DG scheme. The results established in [20] can be adapted directly to the linear HJ equation that if the solutions is adequately smooth in terms of the mixed norm, then the sparse grid DG method using space is convergent of order with a polylogarithmic factor.

3.2 Semi-discrete scheme with multiresolution interpolation

For nonlinear problems, one major difficulty of implementation of formulation (3.2) is to compute the volume integral efficiently and accurately, especially in high dimensions. Naive implementation of numerical quadratures is inefficient due to the hierarchical structure of multiwavelets. To address the challenge, we follow the idea in [41, 25] and interpolate the numerical Hamiltonian by using the multiresolution Lagrange interpolation discussed in Section 2.2. In particular, we have the following modified formulation with interpolation. We find so that for all


By doing so, not only can we apply the unidirectional principle to facilitate the computation, but also fast algorithms can be utilized to improve efficiency. We omit the details regarding the fast algorithm and refer readers to [25].

Note that the interpolation procedure in the scheme formulation plays a role as a high order MRA numerical quadrature. There exist two types of points, namely inner points and interface points [44]. It is observed that the schemes using the interface points are more stable than those using the inner points with the same order accuracy (see, e.g. [25]), and hence we choose to use the interface points in the simulations. We also remark if the Hamiltonian is not smooth and , then to ensure stability, one must employ a very high order quadrature, i.e. large , for accurate computation of the volume integral. This is ascribed to the fact that the large quadrature error due to nonsmoothness of may pollute the numerical viscosity and lead to instability. This drawback is observed in [46], and the authors further coupled a nonlinear limiter to restore stability. In this paper, we propose to properly regularize the Hamiltonian so that the interpolation is adequately accurate for stability. The details will be presented in Section 4. Another possible approach is to add artificial viscosity, as done for solving conservation laws [25].

To preserve the accuracy of the original DG scheme, the interpolation operator needs to reach certain accuracy. Following [9], we can write the DG scheme with interpolation (3.3) into the semi-discrete form as


Here is an operator onto and is a discrete approximation of which satisfies


for all with , , determined by (3.1). Using similar techniques as in [9, 27], we have the following proposition on local truncation error of the sparse grid method with The proof is omitted for brevity.

Proposition 3.1 (Local truncation error analysis).

If the interpolation operator in (3.3) has the accuracy of order for sufficiently smooth functions, then the local truncation error of the semi-discrete DG scheme with interpolation (3.3) is of order . To be more precise, for sufficiently smooth Hamiltonian and function , the sparse grid DG method with interpolation (3.3) has the truncation error:


Here, we use to denote a generic constant that may depend on the solution , but does not depend on

The proposition indicates that, to preserve the order accuracy of the original scheme, we should use Hence, in the simulation we let . Meanwhile, many Hamiltonians are non-smooth functions, and indeed we observe numerically that we need for those cases. This will be further discussed in Section 4.

For time discretization, we employ the third order strong-stability-preserving Runge-Kutta (RK) scheme [43] to advance the semi-discrete scheme (3.4). The adaptive procedure follows the technique developed in [3, 21] to determine the space that dynamically evolves over time. The details are omitted for brevity. The main idea is that in light of the distinguished property of multiwavelets, we keep track of multiwavelet coefficients, i.e. norms of , as an error indicator for refining and coarsening, aiming to efficiently capture the viscosity solution of (1.1) which may develop discontinuous derivatives.

4 Numerical examples

In this section, we present a collection of numerical examples to demonstrate the performance of the proposed adaptive sparse grid LDG method for solving the HJ equation. We consider numerical examples up to with smooth and nonsmooth Hamiltonian, and with smooth and nonsmooth viscosity solutions. Noteworthy, we may need to tune for optimal performance. In particular, we observe that for some numerical tests, we can simply take to achieve satisfactory results and maintain the original accuracy of the DG method, while for the some other tests, we may need to take larger to ensure good performance. In all numerical simulations, the value of is taken between and

Example 4.1.

Consider the following Burgers’ equation in -dimension


with periodic boundary conditions.

At for and for , the solutions are still smooth, and we summarize the convergence study of the sparse grid method with in Tables 4.1-4.2, including the L errors and the associated order of accuracy, with various configurations of and . It is observed that larger leads to smaller error magnitude as expected. Slight order reduction is observed for , and it becomes more severe for . This is because, as time evolves, the viscosity solution of (4.1) develops larger and larger mixed derivatives especially in high dimensions. Hence, it may not be optimal to use the sparse grid space for approximating the viscosity solution. In Tables 4.3-4.4 we report the convergence study for the adaptive method for , respectively. In particular, by fixing the maximum mesh level , two rates of convergence are calculated [3]. The first one is with respect to the error threshold:

and the second is with respect to DoF:

For the full grid counterpart, we have for smooth solutions. It is observed that , which is similar to the Burgers’ equation [25]. Furthermore, using larger is beneficial, as the method with larger requires less DoF to attain a certain level of accuracy.

At for and for , the viscosity solutions have developed discontinuous derivatives. In Figure 4.1, we plot the solution profiles computed by the sparse grid method and the adaptive method for . We set , and , and the maximum mesh level for the adaptive method. It is observed the sparse grid method is able to capture the main structure of the solution, but severe oscillations appear due to lack of mesh resolution around the corners, while the adaptive method is able to effectively capture the viscosity solution by adding more DoF in the nonsmooth region. Hence, the sparse grid method with fixed space cannot reliably approximate the nonsmooth viscosity solution. Afterwards, we will only focus on the performance of adaptive method. For , we set , and , and the maximum mesh level , and plot the results generated by the adaptive method in Figure 4.2, including the 2D cuts of the solution and the associated active elements at final time. Similar results to are observed.

error order error order error order
3 2.63E-02 1.99E-02 1.99E-02
4 7.87E-03 1.74 5.88E-03 1.76 5.75E-03 1.79
5 3.82E-03 1.04 2.42E-03 1.28 2.24E-03 1.36
6 1.77E-03 1.11 8.27E-04 1.55 6.64E-04 1.76
7 7.95E-04 1.16 3.47E-04 1.25 2.10E-04 1.66
error order error order error order
3 5.68E-03 2.84E-03 2.81E-03
4 1.21E-03 2.23 3.75E-04 2.92 3.66E-04 2.94
5 3.10E-04 1.97 1.42E-04 1.40 1.41E-04 1.37
6 5.28E-05 2.56 1.97E-05 2.84 1.96E-05 2.85
7 9.58E-06 2.46 4.36E-06 2.18 4.32E-06 2.18
error order error order error order
3 1.14E-03 6.28E-04 6.10E-04
4 1.68E-04 2.76 6.86E-05 3.19 6.58E-05 3.21
5 2.59E-05 2.70 1.45E-05 2.24 1.44E-05 2.20
6 2.10E-06 3.63 7.94E-07 4.19 7.84E-07 4.20
7 2.77E-07 2.92 1.50E-07 2.41 1.49E-07 2.39
Table 4.1: Example 4.1, . Sparse grid. .
error order error order error order
3 3.34E-02 9.22E-03 7.46E-03
4 1.71E-02 0.96 3.69E-03 1.32 3.24E-03 1.20
5 6.93E-03 1.31 1.32E-03 1.49 1.21E-03 1.43
6 2.07E-03 1.75 5.49E-04 1.26 5.28E-04 1.19
7 8.07E-04 1.36 1.57E-04 1.80 1.54E-04 1.78
8 1.91E-04 2.08 4.09E-05 1.94 3.98E-05 1.95
error order error order error order
3 3.34E-02 9.22E-03 7.46E-03
4 1.71E-02 0.96 3.69E-03 1.32 3.24E-04 1.20
5 6.93E-03 1.31 1.32E-03 1.49 1.21E-03 1.43
6 2.07E-03 1.75 5.49E-04 1.26 5.28E-04 1.19
7 8.07E-04 1.36 1.57E-04 1.80 1.54E-04 1.78
8 1.91E-04 2.08 4.09E-05 1.94 3.98E-05 1.95
error order error order error order
3 8.86E-03 3.56E-03 2.48E-03
4 2.97E-03 1.58 1.10E-03 1.70 8.62E-04 1.53
5 9.97E-04 1.57 3.64E-04 1.59 2.93E-04 1.56
6 3.08E-04 1.70 9.78E-05 1.90 8.57E-05 1.77
7 6.49E-05 2.24 2.17E-05 2.17 2.04E-05 2.07
8 1.44E-05 2.17 5.02E-06 2.12 4.85E-06 2.07
Table 4.2: Example 4.1, . Sparse grid. .
DoF L-error
1.00E-03 448 1.56E-03
1.00E-04 1376 6.92E-04 0.35 0.73
1.00E-05 3520 2.55E-04 0.43 1.06
1.00E-06 10240 2.26E-05 1.05 2.27
1.00E-07 18688 1.05E-05 0.33 1.28
1.00E-03 270 6.67E-04
1.00E-04 720 3.10E-04 0.33 0.78
1.00E-05 1548 4.93E-05 0.80 2.40
1.00E-06 3492 1.34E-05 0.57 1.60
1.00E-07 7704 1.83E-06 0.86 2.51
1.00E-03 192 1.14E-03
1.00E-04 480 1.04E-04 1.04 2.62
1.00E-05 896 3.14E-05 0.52 1.92
1.00E-06 1856 7.14E-06 0.64 2.03
1.00E-07 3136 7.07E-07 1.00 4.41
Table 4.3: Example 4.1, . Adaptive sparse grid. . .
DoF L-error
1.00E-03 2432 7.87E-03
1.00E-04 14864 3.03E-03 0.41 0.53
1.00E-05 44656 1.17E-03 0.41 0.87
1.00E-06 152176 3.25E-04 0.56 1.04
1.00E-07 380976 9.05E-05 0.56 1.39
1.00E-03 2646 5.84E-03
1.00E-04 8208 9.84E-04 0.77 1.57
1.00E-05 21816 1.96E-04 0.70 1.65
1.00E-06 55404 6.11E-05 0.51 1.25
1.00E-07 133569 1.33E-05 0.66 1.74
1.00E-03 2048 2.01E-03
1.00E-04 6400 5.01E-04 0.60 1.22
1.00E-05 16384 9.26E-05 0.73 1.80
1.00E-06 35584 2.30E-05 0.60 1.79
1.00E-07 99840 2.82E-06 0.91 2.04
Table 4.4: Example 4.1, . Adaptive sparse grid. . .
Figure 4.1: Example 4.1, . , . . . =. (a) Numerical solution by sparse grids. (b) Numerical solutions by adaptive sparse grid. (c) Active elements.
Figure 4.2: Example 4.1, . , . . . =. (a) 2D-cuts of the numerical solution at . (b) Active elements.
Example 4.2.

Consider the following HJ equation with a nonconvex Hamiltonian


with periodic boundary conditions.

In Table 4.5-4.6, we report the convergence rates for the adaptive method for and at and , respectively. Similar results are observed to the previous example. In Figure 4.3, we report the solution profile together with the active elements used at when the viscosity solution has developed nonsmooth corners. In this simulations, we set and . Again, the adaptive method is able to efficiently and correctly capture the sharp corners. In Figure 4.4, we plot the results for at with configuration parameters , , maximum level , and . High resolution result is observed.

DoF L-error
1.00E-03 464 1.47E-03
1.00E-04 1616 4.60E-04 0.51 0.93
1.00E-05 3840 1.66E-04 0.44 1.18
1.00E-06 9056 2.37E-05 0.85 2.27
1.00E-07 17440 7.86E-06 0.48 1.68
1.00E-03 288 1.43E-03
1.00E-04 720 3.20E-04 0.65 1.64
1.00E-05 1656 9.24E-05 0.54 1.49
1.00E-06 3924 1.79E-05 0.71 1.90
1.00E-07 8406 4.00E-06 0.65 1.97
1.00E-03 192 1.53E-03
1.00E-04 512 1.85E-04 0.92 2.15
1.00E-05 960 4.56E-05 0.61 2.23
1.00E-06 2048 1.50E-05 0.48 1.47
1.00E-07 3968 1.64E-06 0.96 3.34
Table 4.5: Example 4.2, . Adaptive sparse grid. . .
DoF L-error
1.00E-03 2432 5.22E-03
1.00E-04 15680 2.70E-03 0.29 0.35
1.00E-05 46768 1.11E-03 0.39 0.82
1.00E-06 151480 4.62E-04 0.38 0.74
1.00E-07 391008 1.66E-04 0.45 1.08
1.00E-03 2646 5.61E-03
1.00E-04 8127 1.87E-03 0.48 0.98
1.00E-05 21492 6.16E-04 0.48 1.15
1.00E-06 55296 1.24E-04 0.70 1.70
1.00E-07 154926 3.75E-05 0.52 1.16
1.00E-03 1664 2.87E-03
1.00E-04 6400 1.09E-03 0.42 0.72
1.00E-05 17920 1.60E-04 0.83 1.87
1.00E-06 44416 4.19E-05 0.58 1.47
1.00E-07 186368 5.67E-06 0.87 1.40
Table 4.6: Example 4.2, . Adaptive sparse grid. . .
Figure 4.3: Example 4.2, . . , . . =. (a) Numerical solutions on adaptive grids. (b) Active elements.
Figure 4.4: Example 4.2, . , . . . =. (a) 2D-cuts of the numerical solution at . (b) Active elements.
Example 4.3.

We consider the following two-dimensional nonlinear problem


with periodic boundary conditions.

Note that unlike the previous two examples, the problem is genuinely nonlinear, and the Hamiltonian is smooth but nonconvex. When , the solution is still smooth, and we are able to test the convergence for the adaptive method. In the simulation, we set maximum level , , . It is observed in Table 4.7 that the method is able to achieve very accurate results by using a few DoFs. The convergence performance is similar to the previous examples. In Figure 4.5, we plot the solution at , when the viscosity solution becomes nonsmooth. It is observed that the adaptive method captures the corners correctly and efficiently, as compared with the results by other popular methods, see e.g. [34, 6].

DoF L-error
1.00E-03 180 9.78E-03
1.00E-04 448 2.05E-03 0.68 1.71
1.00E-05 952 1.26E-03 0.21 0.64
1.00E-06 1296 2.24E-04 0.75 5.60
1.00E-07 2952 2.74E-05 0.91 2.56
1.00E-03 135 2.73E-03
1.00E-04 306 3.95E-04 0.84 2.36
1.00E-05 594 1.97E-04 0.30 1.05
1.00E-06 1224 4.41E-05 0.65 2.07
1.00E-07 2565 1.31E-05 0.53 1.64
1.00E-03 112 5.38E-04
1.00E-04 256 1.67E-04 0.51 1.42
1.00E-05 560 4.29E-05 0.59 1.73
1.00E-06 832 1.21E-05 0.55 3.20
1.00E-07 1280 1.25E-06 0.99 5.27
Table 4.7: Example 4.3. Adaptive sparse grid. . .
Figure 4.5: Example 4.3. . . , . =. (a) Numerical solution profile. (b) Active elements.
Example 4.4.

We consider the classic nonlinear Eikonal equation


where and

An outflow boundary condition is imposed. The viscosity solution is

which is clearly smooth.

One additional challenge of this problem is that the Hamiltonian is not smooth, making the DG formulation unstable if the numerical quadrature is not sufficiently accurate, as mentioned in previous section. To circumvent the difficulty, we propose to employ a regularized Hamiltonian as follows.


It can be easily verified that is . In the simulation, we choose , where is the mesh size, hence the regularization will not affect the accuracy of the original method. We employ the regularized Hamiltonian for all tests, while we notice that it is only required for . In Table 4.8, we summarize the convergence study for the adaptive method for and . It is observed that the convergence rates and are similar for and , which is unsurprising, since the viscosity solution is only . Meanwhile, the error magnitude by is still much smaller than that by with the same number of DoF, demonstrating the efficiency of method with high order accuracy. In Figure 4.6, we report the contour plot of the numerical solution with , , , , . We observe that the rarefaction wave developed at the center of domain is correctly captured by the adaptive method. We also highlight the level set of .

DoF L-error DoF L-error
1.00E-03 236 2.25E-02 72 5.42E-03
1.00E-04 496 5.39E-03 0.62 1.92 108 3.59E-03 0.18 1.01
1.00E-05 1056 2.93E-03 0.26 0.81 324 1.09E-03 0.52 1.09
1.00E-06 1904 9.27E-04 0.50 1.95 900 4.41E-04 0.39 0.89
1.00E-07 5496 2.43E-04 0.58 1.26 2880 1.23E-04 0.56 1.10
1.00E-03 680 2.31E-02 108 6.46E-03
1.00E-04 1472 7.64E-03 0.48 1.43 351 3.18E-03 0.31 0.60
1.00E-05 2968 4.12E-03 0.27 0.88 1026 1.53E-03 0.32 0.68
1.00E-06 5080 1.72E-03 0.38 1.63 2970 5.76E-04 0.43 0.92
1.00E-07 23272 5.00E-04 0.54 0.81 11610 2.26E-04 0.41 0.69
1.00E-03 1872 2.30E-02 405 4.61E-03
1.00E-04 3792 2.52E-02 -0.04 -0.13 1053 2.90E-03 0.20 0.48
1.00E-05 8944 1.25E-02 0.30 0.81 3159 1.16E-03 0.40 0.84
1.00E-06 10624 2.59E-03 0.68 9.15 12312 5.71E-04 0.31 0.52
1.00E-07 - - - - 55080 2.03E-04 0.45 0.69
Table 4.8: Example 4.4, . Adaptive sparse grid. . .
Figure 4.6: Example 4.4. . , . . =. (a) Contour plot of the numerical solution. (b) Numerical error distribution. (c) Active elements.
Example 4.5.

In this example, we consider the following HJB equation [3]


where and is a set of vectors corresponding to possible controls. The function is the same as in the example 4.4. Note that this HJB equation is equivalent to the following HJ equation


The exact solution can be hence derived from (4.7):

Here, for a vector , in the component-wise sense. We apply the adaptive algorithm to simulate (4.7). The outflow boundary conditions are imposed. Note that the Hamiltonian is nonsmooth as with the Eikonal equation, and hence we regularize the absolute function using the technique (4.5) to ensure stability. In Figure 4.7, we plot the solution with configuration , , , . Note that the viscosity solution is , a rarefaction wave opens up at the center of the domain, which is well captured by the method. In Table 4.9, we summarize the convergence study for and . Note that when , the error does not decay anymore, since it has saturated already with the maximum level .

DoF L-error DoF L-error
1.00E-03 204 4.17E-03 63 8.62E-03
1.00E-04 444 1.62E-03 0.41 1.21 135 1.89E-03 0.66 1.99
1.00E-05 860 7.01E-04 0.36 1.27 207 6.26E-04 0.48 2.59
1.00E-06 876 6.43E-04 0.04 4.69 459 4.24E-04 0.17 0.49
1.00E-07 924 6.58E-04 -0.01 -0.43 855 3.97E-04 0.03 0.11
1.00E-03 608 5.44E-03 270 1.12E-02
1.00E-04 1328 2.12E-03 0.41 1.20 594 2.98E-03 0.58 1.68
1.00E-05 2576 9.15E-04 0.37 1.27 918 7.89E-04 0.58 3.05
1.00E-06 2624 8.39E-04 0.04 4.74 2052 4.36E-04 0.26 0.74
1.00E-07 2768 8.56E-04 -0.01 -0.38 3510 3.93E-04 0.05 0.19
1.00E-03 1616 6.65E-03 1053 1.35E-02
1.00E-04 3536 2.60E-03 0.41 1.20 1701 6.17E-03 0.34 1.63
1.00E-05 6864 1.12E-03 0.37 1.27 2997 1.24E-03 0.70 2.84
1.00E-06 6992 1.02E-03 0.04 4.69 6237 4.90E-04 0.40 1.26
1.00E-07 7376 1.05E-03 -0.01 -0.38 12069 4.20E-04 0.07 0.23
Table 4.9: Example 4.5, . Adaptive sparse grid. . .
Figure 4.7: Example 4.5. . , . . =. (a) Contour plot of the numerical solution. (b) Numerical error distribution. (c) Active elements.
Example 4.6.

In the last example, we consider the 2D problem related to controlling optimal cost determination [39]


Note that the Hamiltonian is not smooth. In Figure 4.8, we plot the solution profile, the optimal together with the active elements at final time . Again, we regularize the Hamiltonian as with previous examples. The adaptive method is able to capture the viscosity solution efficiently, and the numerical results agree with other methods in the literature, e.g. [24, 34, 6, 22, 30].

Figure 4.8: Example 4.6, . . , . . =. (a) Numerical solution profile. (b) Active elements. (c) Controls

5 Conclusion

In this work, we proposed an adaptive sparse grid LDG method for solving HJ equations in high dimensions. By incorporating the orthonormal Alpert’s multiwavelets as the DG finite element bases, and the interpolatory multiwavelets as efficient multiresolution numerical quadratures, we achieve efficient multiresolution schemes which is suitable for high dimensions. Benchmark numerical tests up to 4D are provided to validate the performance of the method. The code generating the results in this paper can be found at the GitHub link:, and it has the capability of computing higher dimensional problems.


We would like to thank Qi Tang and Kai Huang for the assistance and discussion in code implementation. Yingda Cheng would like to thank the support from IPAM to attend the workshop on “High-dimensional Hamilton-Jacobi PDEs”.


  • [1] R. Abgrall. Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes. Communications on Pure and Applied Mathematics, 49(12):1339–1373, 1996.
  • [2] B. Alpert. A class of bases in for the sparse representation of integral operators. SIAM Journal on Mathematical Analysis, 24(1):246–262, 1993.
  • [3] O. Bokanowski, J. Garcke, M. Griebel, and I. Klompmaker. An adaptive sparse grid semi-Lagrangian scheme for first order Hamilton-Jacobi Bellman equations. Journal of Scientific Computing, 55(3):575–605, 2013.
  • [4] H.-J. Bungartz and M. Griebel. Sparse Grids. Acta Numerica, 13:147–269, 2004.
  • [5] Y. Cheng and C.-W. Shu. A discontinuous Galerkin finite element method for directly solving the Hamilton–Jacobi equations. Journal of Computational Physics, 223(1):398–415, 2007.
  • [6] Y. Cheng and Z. Wang. A new discontinuous Galerkin finite element method for directly solving the Hamilton-Jacobi equations. Journal of Computational Physics, 268:134–153, 2014.
  • [7] Y. Chow, J. Darbon, S. Osher, and W. Yin. Algorithm for overcoming the curse of dimensionality for time-dependent non-convex Hamilton–Jacobi equations arising from optimal control and differential games problems. Journal of Scientific Computing, 73(2-3):617–643, 2017.
  • [8] Y. T. Chow, J. Darbon, S. Osher, and W. Yin. Algorithm for overcoming the curse of dimensionality for state-dependent Hamilton-Jacobi equations. Journal of Computational Physics, 387:376–409, 2019.
  • [9] B. Cockburn, S. Hou, and C.-W. Shu. The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. IV. The multidimensional case. Mathematics of Computation, 54(190):545–581, 1990.
  • [10] B. Cockburn and C.-W. Shu. The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM Journal on Numerical Analysis, 35(6):2440–2463, 1998.
  • [11] M. G. Crandall, L. C. Evans, and P.-L. Lions. Some properties of viscosity solutions of Hamilton-Jacobi equations. Transactions of the American Mathematical Society, 282(2):487–502, 1984.
  • [12] M. G. Crandall and P.-L. Lions. Viscosity solutions of Hamilton-Jacobi equations. Transactions of the American Mathematical Society, 277(1):1–42, 1983.
  • [13] M. G. Crandall and P. L. Lions. Two approximations of solutions of Hamilton–Jacobi equations. Mathematics of Computation, 43(167):1–19, 1984.
  • [14] J. Darbon, G. P. Langlois, and T. Meng. Overcoming the curse of dimensionality for some Hamilton–Jacobi partial differential equations via neural network architectures. arXiv preprint arXiv:1910.09045, 2019.
  • [15] J. Darbon and T. Meng. On some neural network architectures that can represent viscosity solutions of certain high dimensional Hamilton–Jacobi partial differential equations. arXiv preprint arXiv:2002.09750, 2020.
  • [16] J. Darbon and S. Osher. Algorithms for overcoming the curse of dimensionality for certain Hamilton–Jacobi equations arising in control theory and elsewhere. Research in the Mathematical Sciences, 3(1):19, 2016.
  • [17] S. Dolgov, D. Kalise, and K. Kunisch. Tensor decompositions for high-dimensional Hamilton-Jacobi-Bellman equations. arXiv preprint arXiv:1908.01533, 2019.
  • [18] L. Evans. Partial Differential Equations: Second Edition. American Mathematical Society, 2010.
  • [19] J. Garcke and A. Kröner. Suboptimal feedback control of pdes by solving hjb equations on adaptive sparse grids. Journal of Scientific Computing, 70(1):1–28, 2017.
  • [20] W. Guo and Y. Cheng. A sparse grid discontinuous galerkin method for high-dimensional transport equations and its application to kinetic simulations. SIAM Journal on Scientific Computing, 38(6):A3381–A3409, 2016.
  • [21] W. Guo and Y. Cheng. An adaptive multiresolution discontinuous Galerkin method for time-dependent transport equations in multidimensions. SIAM Journal on Scientific Computing, 39(6):A2962–A2992, 2017.
  • [22] W. Guo, F. Li, and J. Qiu. Local-structure-preserving discontinuous Galerkin methods with Lax-Wendroff type time discretizations for Hamilton-Jacobi equations. Journal of Scientific Computing, 47(2):239–257, 2011.
  • [23] J. Han, A. Jentzen, and E. Weinan.

    Solving high-dimensional partial differential equations using deep learning.

    Proceedings of the National Academy of Sciences, 115(34):8505–8510, 2018.
  • [24] C. Hu and C.-W. Shu. A discontinuous Galerkin finite element method for Hamilton–Jacobi equations. SIAM Journal on Scientific Computing, 21(2):666–690, 1999.
  • [25] J. Huang and Y. Cheng. An adaptive multiresolution discontinuous Galerkin method with artificial viscosity for scalar hyperbolic conservation laws in multidimensions. arXiv preprint arXiv:1906.00829, 2019.
  • [26] J. Huang, Y. Liu, W. Guo, Z. Tao, and Y. Cheng. An adaptive multiresolution interior penalty discontinuous Galerkin method for wave equations in second order form. arXiv preprint arXiv:2004.08525, 2020.
  • [27] J. Huang and C.-W. Shu.

    Error estimates to smooth solutions of semi-discrete discontinuous Galerkin methods with quadrature rules for scalar conservation laws.

    Numerical Methods for Partial Differential Equations, 33(2):467–488, 2017.
  • [28] G.-S. Jiang and D. Peng. Weighted ENO schemes for Hamilton–Jacobi equations. SIAM Journal on Scientific computing, 21(6):2126–2143, 2000.
  • [29] W. Kang and L. C. Wilcox. Mitigating the curse of dimensionality: sparse grid characteristics method for optimal feedback control and hjb equations. Computational Optimization and Applications, 68(2):289–315, 2017.
  • [30] G. Ke and W. Guo. An alternative formulation of discontinous Galerkin schemes for solving Hamilton–Jacobi equations. Journal of Scientific Computing, 78(2):1023–1044, 2019.
  • [31] K. Kunisch, S. Volkwein, and L. Xie. HJB-POD-based feedback design for the optimal control of evolution problems. SIAM Journal on Applied Dynamical Systems, 3(4):701–722, 2004.
  • [32] F. Lafon and S. Osher. High order two dimensional nonoscillatory methods for solving Hamilton–Jacobi scalar equations. Journal of Computational Physics, 123(2):235–253, 1996.
  • [33] O. Lepsky, C. Hu, and C.-W. Shu. Analysis of the discontinuous Galerkin method for Hamilton–Jacobi equations. Applied Numerical Mathematics, 33(1-4):423–434, 2000.
  • [34] F. Li and S. Yakovlev. A central discontinuous Galerkin method for Hamilton-Jacobi equations. Journal of Scientific Computing, 45(1-3):404–428, 2010.
  • [35] P.-L. Lions. Generalized solutions of Hamilton-Jacobi equations, volume 69. Pitman, London, 1982.
  • [36] S. Mallat. A Wavelet Tour of Signal Processing. Elsevier, 1999.
  • [37] T. Nakamura-Zimmerer, Q. Gong, and W. Kang. Adaptive deep learning for high dimensional Hamilton-Jacobi-Bellman equations. arXiv preprint arXiv:1907.05317, 2019.
  • [38] S. Osher and J. Sethian. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 79(1):12–49, 1988.
  • [39] S. Osher and C.-W. Shu. High-order essentially nonoscillatory schemes for Hamilton–Jacobi equations. SIAM Journal on Numerical Analysis, 28(4):907–922, 1991.
  • [40] J. Qiu and C.-W. Shu. Hermite WENO schemes for Hamilton–Jacobi equations. Journal of Computational Physics, 204(1):82–99, 2005.
  • [41] J. Shen and H. Yu. Efficient spectral sparse grid methods and applications to high-dimensional elliptic problems. SIAM Journal on Scientific Computing, 32(6):3228–3250, 2010.
  • [42] C.-W. Shu. High order numerical methods for time dependent Hamilton-Jacobi equations. Mathematics and Computation in Imaging Science and Information Processing, Lect. Notes Ser. Inst. Math. Sci. Natl. Univ. Singap, 11:47–91, 2007.
  • [43] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics, 77(2):439–471, 1988.
  • [44] Z. Tao, Y. Jiang, and Y. Cheng. An adaptive high-order piecewise polynomial based sparse grid collocation method with applications. arXiv preprint arXiv:1912.03982, 2019.
  • [45] Z. Wang, Q. Tang, W. Guo, and Y. Cheng. Sparse grid discontinuous Galerkin methods for high-dimensional elliptic equations. Journal of Computational Physics, 314:244–263, 2016.
  • [46] J. Yan and S. Osher. A local discontinuous Galerkin method for directly solving Hamilton–Jacobi equations. Journal of Computational Physics, 230(1):232–244, 2011.
  • [47] Y.-T. Zhang and C.-W. Shu. High-order WENO schemes for Hamilton–Jacobi equations on triangular meshes. SIAM Journal on Scientific Computing, 24(3):1005–1030, 2003.