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.


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

  • $\begingroup$ This answer is amazing. Thank you very much. $\endgroup$ Commented Dec 13, 2018 at 10:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.