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

Question

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.

Answer 1

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$.