1 Introduction
In this paper, we consider the HamiltonJacobi (HJ) equation
(1.1) 
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 nonoscillatory (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 highdimensional 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 smallscale 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 tensorproduct 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
(2.1) 
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 tensorproduct approach to construct the hierarchical finite element space in multidimensional space. Denote as the mesh level in a multivariate sense, where denotes the set of nonnegative integers, we can define the tensorproduct mesh grid and the corresponding mesh size Based on the grid , we denote as an elementary cell, and
as the tensorproduct 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 tensorproduct construction, the multidimensional increment space can be defined as
The basis functions in multidimensions are defined as
(2.2) 
for , and .
Using the notation of
and the same componentwise arithmetic operations and relations as defined in [45], we reach the decomposition
(2.3) 
On the other hand, a standard choice of sparse grid space [45, 20] is
(2.4) 
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.
(2.5) 
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
(2.6) 
and is defined through the multiwavelets that satisfies
for . Then is given by
The multidimensional construction follows similar lines as in Section 2.1. We let
then
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 multidimensional basis functions are defined in the same approach as (2.2) by tensor products:
(2.7) 
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 nonperiodic 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 onesided numerical fluxes; that is, given we seek , , in such that for all
(3.1) 
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 semidiscrete scheme for solving (1.1) is defined as follows: seek such that, for all ,
(3.2) 
where denotes a monotone numerical Hamiltonian that approximates , and , , are given in (3.1). In the simulations, we employ the following global LaxFriedrichs 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 tensorproduct 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 Semidiscrete 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
(3.3) 
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 semidiscrete form as
(3.4) 
Here is an operator onto and is a discrete approximation of which satisfies
(3.5) 
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 semidiscrete 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:
(3.6) 
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 nonsmooth 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 strongstabilitypreserving RungeKutta (RK) scheme [43] to advance the semidiscrete 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
(4.1) 
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.14.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.34.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.63E02  –  1.99E02  –  1.99E02  –  
4  7.87E03  1.74  5.88E03  1.76  5.75E03  1.79  
5  3.82E03  1.04  2.42E03  1.28  2.24E03  1.36  
6  1.77E03  1.11  8.27E04  1.55  6.64E04  1.76  
7  7.95E04  1.16  3.47E04  1.25  2.10E04  1.66  
error  order  error  order  error  order  
3  5.68E03  –  2.84E03  –  2.81E03  –  
4  1.21E03  2.23  3.75E04  2.92  3.66E04  2.94  
5  3.10E04  1.97  1.42E04  1.40  1.41E04  1.37  
6  5.28E05  2.56  1.97E05  2.84  1.96E05  2.85  
7  9.58E06  2.46  4.36E06  2.18  4.32E06  2.18  
error  order  error  order  error  order  
3  1.14E03  –  6.28E04  –  6.10E04  –  
4  1.68E04  2.76  6.86E05  3.19  6.58E05  3.21  
5  2.59E05  2.70  1.45E05  2.24  1.44E05  2.20  
6  2.10E06  3.63  7.94E07  4.19  7.84E07  4.20  
7  2.77E07  2.92  1.50E07  2.41  1.49E07  2.39  
error  order  error  order  error  order  

3  3.34E02  –  9.22E03  –  7.46E03  –  
4  1.71E02  0.96  3.69E03  1.32  3.24E03  1.20  
5  6.93E03  1.31  1.32E03  1.49  1.21E03  1.43  
6  2.07E03  1.75  5.49E04  1.26  5.28E04  1.19  
7  8.07E04  1.36  1.57E04  1.80  1.54E04  1.78  
8  1.91E04  2.08  4.09E05  1.94  3.98E05  1.95  
error  order  error  order  error  order  
3  3.34E02  –  9.22E03  –  7.46E03  –  
4  1.71E02  0.96  3.69E03  1.32  3.24E04  1.20  
5  6.93E03  1.31  1.32E03  1.49  1.21E03  1.43  
6  2.07E03  1.75  5.49E04  1.26  5.28E04  1.19  
7  8.07E04  1.36  1.57E04  1.80  1.54E04  1.78  
8  1.91E04  2.08  4.09E05  1.94  3.98E05  1.95  
error  order  error  order  error  order  
3  8.86E03  –  3.56E03  –  2.48E03  –  
4  2.97E03  1.58  1.10E03  1.70  8.62E04  1.53  
5  9.97E04  1.57  3.64E04  1.59  2.93E04  1.56  
6  3.08E04  1.70  9.78E05  1.90  8.57E05  1.77  
7  6.49E05  2.24  2.17E05  2.17  2.04E05  2.07  
8  1.44E05  2.17  5.02E06  2.12  4.85E06  2.07 
DoF  Lerror  

1.00E03  448  1.56E03  –  –  
1.00E04  1376  6.92E04  0.35  0.73  
1.00E05  3520  2.55E04  0.43  1.06  
1.00E06  10240  2.26E05  1.05  2.27  
1.00E07  18688  1.05E05  0.33  1.28  
1.00E03  270  6.67E04  –  –  
1.00E04  720  3.10E04  0.33  0.78  
1.00E05  1548  4.93E05  0.80  2.40  
1.00E06  3492  1.34E05  0.57  1.60  
1.00E07  7704  1.83E06  0.86  2.51  
1.00E03  192  1.14E03  –  –  
1.00E04  480  1.04E04  1.04  2.62  
1.00E05  896  3.14E05  0.52  1.92  
1.00E06  1856  7.14E06  0.64  2.03  
1.00E07  3136  7.07E07  1.00  4.41 
DoF  Lerror  

1.00E03  2432  7.87E03  –  –  
1.00E04  14864  3.03E03  0.41  0.53  
1.00E05  44656  1.17E03  0.41  0.87  
1.00E06  152176  3.25E04  0.56  1.04  
1.00E07  380976  9.05E05  0.56  1.39  
1.00E03  2646  5.84E03  –  –  
1.00E04  8208  9.84E04  0.77  1.57  
1.00E05  21816  1.96E04  0.70  1.65  
1.00E06  55404  6.11E05  0.51  1.25  
1.00E07  133569  1.33E05  0.66  1.74  
1.00E03  2048  2.01E03  –  –  
1.00E04  6400  5.01E04  0.60  1.22  
1.00E05  16384  9.26E05  0.73  1.80  
1.00E06  35584  2.30E05  0.60  1.79  
1.00E07  99840  2.82E06  0.91  2.04 
Example 4.2.
Consider the following HJ equation with a nonconvex Hamiltonian
(4.2) 
with periodic boundary conditions.
In Table 4.54.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  Lerror  

1.00E03  464  1.47E03  
1.00E04  1616  4.60E04  0.51  0.93  
1.00E05  3840  1.66E04  0.44  1.18  
1.00E06  9056  2.37E05  0.85  2.27  
1.00E07  17440  7.86E06  0.48  1.68  
1.00E03  288  1.43E03  
1.00E04  720  3.20E04  0.65  1.64  
1.00E05  1656  9.24E05  0.54  1.49  
1.00E06  3924  1.79E05  0.71  1.90  
1.00E07  8406  4.00E06  0.65  1.97  
1.00E03  192  1.53E03  
1.00E04  512  1.85E04  0.92  2.15  
1.00E05  960  4.56E05  0.61  2.23  
1.00E06  2048  1.50E05  0.48  1.47  
1.00E07  3968  1.64E06  0.96  3.34 
DoF  Lerror  

1.00E03  2432  5.22E03  –  –  
1.00E04  15680  2.70E03  0.29  0.35  
1.00E05  46768  1.11E03  0.39  0.82  
1.00E06  151480  4.62E04  0.38  0.74  
1.00E07  391008  1.66E04  0.45  1.08  
1.00E03  2646  5.61E03  –  –  
1.00E04  8127  1.87E03  0.48  0.98  
1.00E05  21492  6.16E04  0.48  1.15  
1.00E06  55296  1.24E04  0.70  1.70  
1.00E07  154926  3.75E05  0.52  1.16  
1.00E03  1664  2.87E03  –  –  
1.00E04  6400  1.09E03  0.42  0.72  
1.00E05  17920  1.60E04  0.83  1.87  
1.00E06  44416  4.19E05  0.58  1.47  
1.00E07  186368  5.67E06  0.87  1.40 
Example 4.3.
We consider the following twodimensional nonlinear problem
(4.3) 
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  Lerror  

1.00E03  180  9.78E03  –  –  
1.00E04  448  2.05E03  0.68  1.71  
1.00E05  952  1.26E03  0.21  0.64  
1.00E06  1296  2.24E04  0.75  5.60  
1.00E07  2952  2.74E05  0.91  2.56  
1.00E03  135  2.73E03  –  –  
1.00E04  306  3.95E04  0.84  2.36  
1.00E05  594  1.97E04  0.30  1.05  
1.00E06  1224  4.41E05  0.65  2.07  
1.00E07  2565  1.31E05  0.53  1.64  
1.00E03  112  5.38E04  –  –  
1.00E04  256  1.67E04  0.51  1.42  
1.00E05  560  4.29E05  0.59  1.73  
1.00E06  832  1.21E05  0.55  3.20  
1.00E07  1280  1.25E06  0.99  5.27  
Example 4.4.
We consider the classic nonlinear Eikonal equation
(4.4) 
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.
(4.5) 
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  Lerror  DoF  Lerror  
1.00E03  236  2.25E02  72  5.42E03  
1.00E04  496  5.39E03  0.62  1.92  108  3.59E03  0.18  1.01  
1.00E05  1056  2.93E03  0.26  0.81  324  1.09E03  0.52  1.09  
1.00E06  1904  9.27E04  0.50  1.95  900  4.41E04  0.39  0.89  
1.00E07  5496  2.43E04  0.58  1.26  2880  1.23E04  0.56  1.10  
1.00E03  680  2.31E02  108  6.46E03  
1.00E04  1472  7.64E03  0.48  1.43  351  3.18E03  0.31  0.60  
1.00E05  2968  4.12E03  0.27  0.88  1026  1.53E03  0.32  0.68  
1.00E06  5080  1.72E03  0.38  1.63  2970  5.76E04  0.43  0.92  
1.00E07  23272  5.00E04  0.54  0.81  11610  2.26E04  0.41  0.69  
1.00E03  1872  2.30E02  405  4.61E03  
1.00E04  3792  2.52E02  0.04  0.13  1053  2.90E03  0.20  0.48  
1.00E05  8944  1.25E02  0.30  0.81  3159  1.16E03  0.40  0.84  
1.00E06  10624  2.59E03  0.68  9.15  12312  5.71E04  0.31  0.52  
1.00E07          55080  2.03E04  0.45  0.69  
Example 4.5.
In this example, we consider the following HJB equation [3]
(4.6) 
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
(4.7) 
The exact solution can be hence derived from (4.7):
Here, for a vector , in the componentwise 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  Lerror  DoF  Lerror  
1.00E03  204  4.17E03  63  8.62E03  
1.00E04  444  1.62E03  0.41  1.21  135  1.89E03  0.66  1.99  
1.00E05  860  7.01E04  0.36  1.27  207  6.26E04  0.48  2.59  
1.00E06  876  6.43E04  0.04  4.69  459  4.24E04  0.17  0.49  
1.00E07  924  6.58E04  0.01  0.43  855  3.97E04  0.03  0.11  
1.00E03  608  5.44E03  270  1.12E02  
1.00E04  1328  2.12E03  0.41  1.20  594  2.98E03  0.58  1.68  
1.00E05  2576  9.15E04  0.37  1.27  918  7.89E04  0.58  3.05  
1.00E06  2624  8.39E04  0.04  4.74  2052  4.36E04  0.26  0.74  
1.00E07  2768  8.56E04  0.01  0.38  3510  3.93E04  0.05  0.19  
1.00E03  1616  6.65E03  1053  1.35E02  
1.00E04  3536  2.60E03  0.41  1.20  1701  6.17E03  0.34  1.63  
1.00E05  6864  1.12E03  0.37  1.27  2997  1.24E03  0.70  2.84  
1.00E06  6992  1.02E03  0.04  4.69  6237  4.90E04  0.40  1.26  
1.00E07  7376  1.05E03  0.01  0.38  12069  4.20E04  0.07  0.23  
Example 4.6.
In the last example, we consider the 2D problem related to controlling optimal cost determination [39]
(4.8) 
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].
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: https://github.com/JuntaoHuang/adaptivemultiresolutionDG, and it has the capability of computing higher dimensional problems.
Acknowledgment
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 “Highdimensional HamiltonJacobi PDEs”.
References
 [1] R. Abgrall. Numerical discretization of the firstorder HamiltonJacobi 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 semiLagrangian scheme for first order HamiltonJacobi 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 HamiltonJacobi 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 timedependent nonconvex Hamilton–Jacobi equations arising from optimal control and differential games problems. Journal of Scientific Computing, 73(23):617–643, 2017.
 [8] Y. T. Chow, J. Darbon, S. Osher, and W. Yin. Algorithm for overcoming the curse of dimensionality for statedependent HamiltonJacobi equations. Journal of Computational Physics, 387:376–409, 2019.
 [9] B. Cockburn, S. Hou, and C.W. Shu. The RungeKutta 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 timedependent convectiondiffusion 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 HamiltonJacobi equations. Transactions of the American Mathematical Society, 282(2):487–502, 1984.
 [12] M. G. Crandall and P.L. Lions. Viscosity solutions of HamiltonJacobi 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 highdimensional HamiltonJacobiBellman 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 highdimensional 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 timedependent transport equations in multidimensions. SIAM Journal on Scientific Computing, 39(6):A2962–A2992, 2017.
 [22] W. Guo, F. Li, and J. Qiu. Localstructurepreserving discontinuous Galerkin methods with LaxWendroff type time discretizations for HamiltonJacobi equations. Journal of Scientific Computing, 47(2):239–257, 2011.

[23]
J. Han, A. Jentzen, and E. Weinan.
Solving highdimensional 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 semidiscrete 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. HJBPODbased 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(14):423–434, 2000.
 [34] F. Li and S. Yakovlev. A central discontinuous Galerkin method for HamiltonJacobi equations. Journal of Scientific Computing, 45(13):404–428, 2010.
 [35] P.L. Lions. Generalized solutions of HamiltonJacobi equations, volume 69. Pitman, London, 1982.
 [36] S. Mallat. A Wavelet Tour of Signal Processing. Elsevier, 1999.
 [37] T. NakamuraZimmerer, Q. Gong, and W. Kang. Adaptive deep learning for high dimensional HamiltonJacobiBellman equations. arXiv preprint arXiv:1907.05317, 2019.
 [38] S. Osher and J. Sethian. Fronts propagating with curvaturedependent speed: algorithms based on HamiltonJacobi formulations. Journal of Computational Physics, 79(1):12–49, 1988.
 [39] S. Osher and C.W. Shu. Highorder 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 highdimensional elliptic problems. SIAM Journal on Scientific Computing, 32(6):3228–3250, 2010.
 [42] C.W. Shu. High order numerical methods for time dependent HamiltonJacobi 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 nonoscillatory shockcapturing schemes. Journal of Computational Physics, 77(2):439–471, 1988.
 [44] Z. Tao, Y. Jiang, and Y. Cheng. An adaptive highorder 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 highdimensional 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. Highorder WENO schemes for Hamilton–Jacobi equations on triangular meshes. SIAM Journal on Scientific Computing, 24(3):1005–1030, 2003.