In this short note we investigate the numerical performance of the method of artificial diffusion for second-order fully nonlinear Hamilton-Jacobi-Bellman equations. The method was proposed in (M. Jensen and I. Smears, arxiv:1111.5423); where a framework of finite element methods for Hamilton-Jacobi-Bellman equations was studied theoretically. The numerical examples in this note study how the artificial diffusion is activated in regions of degeneracy, the effect of a locally selected diffusion parameter on the observed numerical dissipation and the solution of second-order fully nonlinear equations on irregular geometries.
Hamilton-Jacobi-Bellman (HJB) equations, which describe the value function in the theory of optimal control, are fully nonlinear partial differential equations, which are of second-order if the underlying control problem is stochastic. One challenge arising in the numerical solution of these equations is the presence of spurious generalised solutions of the PDE which do not coincide with the value function. While these spurious solutions often possess the same regularity as the value function, they violate monotonicity properties exhibited by the value function. These properties lead to the concept of viscosity solutions.
Regarding the numerical solution of HJB equations, we would like to highlight three approaches within the finite element methodology, which have been employed to ensure convergence to the value function. For a review of discrete Markov chain approximations before application of the dynamic programming principle we refer to . For the method of vanishing moments we point to . For the approach by Barles and Souganidis we cite the original source  and the recent adaption  by the authors to a finite element framework. A more comprehensive outline of the literature is in .
2 Numerical Method
Let be a bounded Lipschitz domain in , . Let be a compact metric space and
be continuous. Consider the bounded linear operators of non-negative characteristic form
where belongs to . Then
We assume that the final-time boundary data is non-negative: on . For smooth let , where the supremum is applied pointwise. The HJB equation considered is
Let , , be a sequence of piecewise linear shape-regular finite element spaces, whose underlying meshes are strictly acute. Let , , denote the nodes of the mesh with associated hat functions . Set . The mesh size is denoted . The set of time steps is . Let the th entry of be .
For each and find an approximate splitting into linear operators
of the form , , and . To impose monotonicity, select the artificial diffusion parameters and as in  such that and .
Define, for , ,
Obtain the numerical solution by interpolation of . Then at time is defined, inductively, by interpolating the boundary data and by
where the supremum is understood to be applied component-wise to the collection of vectors
3 Selection of the Artificial Diffusion Parameter
That the diffusion coefficients and in (3) are placed outside of the scalar product originates from the non-divergence form of the linear operators of HJB operators . This structure makes it straightforward to implement local artificial diffusion parameters and , that is to implement a dependence on the nodal position . Indeed a change of corresponds in the assembly of to a scalar multiplication of the th row of the stiffness matrix, see (3).
Changes in the optimal artificial diffusion coefficients arise from variations of the local mesh quality and the local mesh Péclet number. The parameters and may be chosen by studying (35) in  or by examining the assembled matrices of the unstabilised operators. This is in particular simple for the , since only the signs of off-diagonal entries need be corrected to impose a local monotonicity property, thus leading to an algorithm which can be executed row-by-row.
Notice that in contrast, for problems in divergence form, the stabilised diffusion coefficients and need to be determined in the whole domain and not just at nodal positions.
4 An Exact Solution and Convergence Rates
We consider a triangular spatial domain with the vertices , and and the HJB equation
To see that equation (5) is a HJB equation, note that the Euclidean norm of the gradient satisfies
A calculation verifies that the function
is an exact solution of (5), where boundary and final-time conditions are determined by interpolation of the exact solution . The equation is to be solved on the time interval .
We consider a splitting which discretises the advection term explicitly with the minimum amount of diffusion needed for monotonicity. The remaining (linear) diffusion is incorporated in the implicit term, observing that this leads to a time-step restriction . Figure 2 shows the supremum norm error with a constant ratio , showing the uniform stability of the method in this setting.
Figure 1 illustrates a choice of artificial diffusion that is locally adapted to the Péclet number of a coarse mesh. It is seen that diffusion is only artificially introduced in a neighbourhood of the origin, which is where the operator becomes degenerate. Figure 3 illustrates the rates of convergence of the numerical solution in the -, - and -norms to the exact solution at the initial time, now using the constant ratio .
5 Eikonal Equation
The steady-state limit of the time-dependent eikonal equation,
equipped with homogeneous boundary and final-time conditions, measures the distance to the boundary. Due to (6) the eikonal equation belongs to the Hamilton-Jacobi-Bellman family. We consider the equation on a domain (Figure 4, top left) whose irregular shape leads to complicated curves on which is not differentiable. The height of the domain is equal to one.
Figure 4 compares locally-adapted choices of the artificial diffusion parameter with global choices, and illustrates the effect of mesh refinement on the quality of the solution. To compare the quality of the numerical solutions we compare their -norms—recalling that excessive numerical dissipation leads to a smearing out of extrema. The coarse grid solution, with internal nodes and with a global diffusion parameter of , only reaches a height of . In contrast, a computation on the same mesh with a local diffusion parameter, with mean value and standard deviation but maximal value , leads to a height of . With one (two) steps of uniform refinement the number of elements increases by a factor of (of ) and the artificial diffusion decreases to an average value of (of ), giving the improved value of (of ) for the -norm in the lower left (right) part of Figure 4. The -norm of a reference solution on a very fine mesh is .
6 A Second-Order Fully Nonlinear Equation
The final example, on the same domain as in the previous section, is concerned with the fully nonlinear case of second-order. We examine the equation
where is the unit ball in , , , and
The boundary and final-time conditions are homogeneous. As before, the advection term is discretised explicitly, with the locally minimal diffusion needed for monotonicity. The possibly remaining (fully nonlinear) diffusion is placed in the implicit term.
In this example the control of the second-order term is maximised independently of the first-order control ; in the sense that the optimal may be determined without knowledge of . Furthermore, as takes locally either a positive or negative value only the controls and are ever active in the HJB equation. This is an example of the Bang-bang principle. It is also reflected in the left plot in Figure 5, where the value of is plotted. Black colouring signifies maximises the operator, whereas white colouring corresponds to . Observe that no intermediate grey values can appear. The plot of the control mimics some of the features of the value function, which is plotted in the same figure on the right—in part because the Laplacian contains information about the curvature of the solution.
At each time-step of the method, a semi-smooth Newton method  was used to solve the nonlinear discrete equation in (4), where each iteration of the algorithm involves solving a linear system. To study the performance of the algorithm, the HJB equation was solved on a sequence of successively refined meshes, with a constant set of tolerances. The linear systems were solved to a tolerance of by GMRES. The stopping criterion for Newton’s method was a relative residual tolerance of in the maximum norm and a convergence of the iterations requirement of in the maximum norm. The sizes of the linear systems for the sequence of meshes were , , , and . The respective average number of Newton iterations for a single time step were , , , and , with respective standard deviations , , , and . This demonstrates a weak dependence of the number of iterations needed for an individual time step on the system size, thus showing that the total number of linear systems to be solved for a complete computation depends principally on the number of time-steps.
- G. Barles and P.E. Souganidis, Convergence of approximation schemes for fully nonlinear second order equations, J. Asymptotic Analysis, 4:271–283, 1991.
- X. Feng and R. Glowinsky and M. Neilan, Recent developments in numerical methods for fully nonlinear second order partial differential equations, submitted to SIAM Review.
- X. Feng and M. Neilan, The vanishing moment method for fully nonlinear second order partial differential equations: formulation, theory, and numerical analysis, arXiv:1109.1183, 2011.
- M. Jensen and I. Smears, On the convergence of finite element methods for Hamilton-Jacobi-Bellman equations, arXiv:1111.5423, 2011.
- H.J. Kushner and P. Dupuis, Numerical methods for stochastic control problems in continuous time, Applications of Mathematics 24, Springer-Verlag, 2001.