Numerical solution of non-linear first order partial differential equation (HJB)

I am trying to solve a simple optimal control problem using the Hamilton-Jacobi-Bellman equation, numerically in Python. This is proving to be rather difficult as I end up having to solve the following: $$J_t - (J_x)^2 + x\cdot J_x = 0$$ I believe this to be a non-linear first order PDE. Being the HJB, we are given boundary condition at terminal time. $$\mbox{BC:} \quad J(x_f) = \dfrac{1}{4} x_f^2$$ I have attempted this problem myself by simply putting in central Euler approximations for $$J_x$$ and using backward difference for $$J_t$$. I then solve backwards in time by making $$J_{i}^{j-1}$$ the subject of the formula. $$J_{i}^{j-1} = J_{i}^{j}+\frac{k}{4h^2}\big[ J_{i+1}^{j}-J_{i-1}^{j} \big]^2 + \frac{k}{h}\big[ J_{i+1}^{j}-J_{i-1}^{j} \big]x_i$$ This only works for a few computations before the scheme becomes very unstable.

Would you please help guide me on how to solve this. I don't quite know whether my approach is correct. Please suggest resources or help me to develop a scheme that would work.

I haven't solved the HJB equations myself, but I can think of two things you can try based on my experience with other PDEs.

1. Use an "upwinded" approximation for $$J_x$$, instead of a central difference scheme. The $$x J_x$$ term in the equation suggests that information is flowing from left to right in your system as you move forward in time (assuming $$x > 0$$). Therefore it makes more sense to construct the $$J_x$$ approximation using $$J_i$$ and the nodes to its left, since that's where the information is coming from. The simplest such approximation is: $$J_x \approx \frac{1}{h}\left(J_i - J_{i-1}\right)$$ However, I'm not sure how this argument applies to the nonlinear $$(-J_x^2)$$ term. The negative sign on that term suggests a right-to-left direction, so you could probably experiment with a right-biased approximation for $$J_x$$ there.

2. Use an implicit time-marching scheme, instead of an explicit scheme. Explicit time-marching schemes, such as the one you are using, have stability conditions that severely restrict the size of the time-step you can take. In fact, there's a good chance your method will work (or survive for longer) if you reduce the time-step size by some large factor. But for many practical problems, this time-step constraint is simply too restrictive since it increases the cost/time required to obtain the final solution.

In contrast, implicit methods are unconditionally stable for many problems, meaning that the solutions remain bounded for any time-step size you choose. However, there's no free lunch; the extra stability comes at the cost of having to solve a system of equations at each time step. This is because all solution unknowns ($$J_i$$'s) involved in the spatial discretization are now evaluated at time $$j-1$$, instead of at time $$j$$. To make it more clear, here's how your discrete equation would look like with an implicit time-marching method, $$J_i^{j-1} - J_i^{j} = \frac{k}{4h^2} \left[J_{i+1}^{j-1} - J_{i-1}^{j-1} \right]^2 + \frac{k}{\textbf{2}h} \left[J_{i+1}^{j-1} - J_{i-1}^{j-1} \right] x_{i}$$

Note that all terms on the right-hand side are evaluated at time $$(j-1)$$, so your solution unknowns $$J_{i-1}^{j-1}, J_i^{j-1}, J_{i+1}^{j-1}..$$ etc appear on both sides of the equation. The above equation represents a nonlinear system of equations that needs to solved at each time-step, probably using something like the Newton's method.

P.S. I think you are missing a factor of 2 in the last term of your equation. The central difference for $$x J_x$$ needs to be divided by $$2h$$.

• This answer is amazing. Thank you very much. – Dylan Solms Dec 13 '18 at 10:39