Solving an Hamilton Jacobi Bellman equation type /w nonlinear coefficients

Hi everyone,

I'm trying to solve numerically a Hamilton-Jacobi-Bellman PDE with nonlinear coefficients. Since I'm pretty new to the PDE toolbox of Matlab, I would like to share my first thoughts and tries so far, just to make sure I'm heading in the right direction.

The equation I'm trying to solve right now is in the form of

a1(x) + a2(t,x,uy) uy + a3(t,x) ux + a4() uxx + ut = 0

with :

• u a function of time t and two space variables x,y
• ux, uy, ut, uxx defined respectively as the first derivative wrt x, the first derivative wrt y, the first derivative wrt t, the second derivative wrt x
• x a space variable between 0.1 and 1 and y a space variable between 0 and 2 (geom will be a simple rectangular geom)
• t a space variable between 0 and 1.
• a1(x) = 1/x
• a2(t,x,uy) = [1-log2(x.uy)]
• a3(t,x) = x.t
• a4() = K1 (cst)
• an initial condition u0(x,y) = K2.y (K cst)
• boundary conditions quite unclear yet: the only thing I know for sure is that u(t,x,0) = 0, I will probably consider simple Neumann boundary conditions for everything else.

Here are a few thoughts I've had, in order to start this. Just let me know if you think I'm not heading in the right direction or missing something.

1. Defining a rectangular geom + mesh is a piece of cake.
2. Parabolic ( http://www.mathworks.fr/fr/help/pde/ug/parabolic.html ) seems like the function I will need to use to solve this PDE.
3. It sounds like I will have to define coefficients c,f,d,a for parabolic in function form, as described here ( http://www.mathworks.fr/fr/help/pde/ug/scalar-coefficients-in-function-form.html )
4. The coeffcicients a,d,f are simple to express. However, the way c coefficients (which will be used for a2, a3, a4) are defined in my case is quite unclear to me. I found this ( http://www.mathworks.fr/fr/help/pde/ug/c.html ), but I'm having a really hard time understanding what it really means and what is c is supposed to be (matrix ? with whose coefficients ?). Any simple answer on this would be really helpful.
5. I will have to express Dirichlet conditions on the only boundary I know (i.e. u(t,x,0) = 0), Neumann conditions can be added simply in the same custom boundary file, as suggested here ( http://www.mathworks.fr/fr/help/pde/ug/boundary-conditions-for-pde-systems.html )

So far, two questions:

- Am I doing this right ?

- Can anyone provide clarification on point 4 ? (a simple m-file illustrating how it's done or a simple explanation will do perfectly)

Thank you for your time, Best regards,

Matt

Answer 1

[Philip Caplan]

You have the right idea to use "parabolic" to solve your PDE, however, the terms cannot be manipulated into the format necessary for "parabolic". Considering you have a simple rectangular domain, this problem is well-suited to be solved with a finite difference scheme (as opposed to the finite-element discretization used by "parabolic").

I recommend using MATLAB's "ode45" to integrate in the temporal direction and discretize the differential terms with, say, second-order difference operators. Using the notation in your question, the equation is:

ut = -a1(x) -a2(t,x,uy)uy -a3(t,x)ux -a4()uxx

A simple nx by ny grid for the rectangular domain can be set up with

>> [x,y] = meshgrid(linspace(.1,1,nx),linspace(0,2,ny));

Assuming you have the solution at the previous time step, u, the second derivative (in x) at all interior points can be computed using

>> uxx(2:end-1,:) = ( u(3:end,:) -2*u(2:end-1,:) +u(1:end-2,:) )/dx^2;

where dx is the spacing between grid points in the x-direction (here, 0.9/(nx-1)).

You will need to be careful with points along the boundary and correctly impose boundary conditions, but the idea is the same. With all these difference operators (at all points in the grid), the right-hand side can be provided to "ode45" and you will obtain the time-dependent solution.

[Matthieu]

Thank you for your time.

I think I see what you are suggesting. I've never heard of this ode45 approach before, but I have been reading the documentation about it ( http://www.mathworks.fr/fr/help/matlab/ref/ode45.html ). I understand that I will have to define an odefun, modeling the right side of this equation

ut = -a1(x) -a2(t,x,uy)uy -a3(t,x)ux -a4()uxx = odefun(t,u)

and I will pass this odefun function to ode45.

Is the following code close to what you have in mind for odefun ?

     function dudt = odefun(t,u)

% How I should define the core of odefun

    dudt = zeros(1,1);
  dudt(1) = -a1(x) -a2(t,x,uy)uy -a3(t,x)ux -a4()uxx

% And obviously, uxx, ux and uy will have been previously defined for interior points as:

uxx = can be computed using this trick ( u(3:end,:) -2*u(2:end-1,:) +u(1:end-2,:) )/dx^2; 
ux = blablabla..
uy = blablabla..

% When not an interior point, I need to refer to the boundary conditions I have defined to express (But how ?)

uxx
ux
uy

    end

I could really use an example of code implementing ode45 for such a PDE problem. I admit I haven't been looking for any yet (I won't do it today, as it is quite late now and I really wanted to answer you before I go to bed), but if you have any good example/tutorial you could forward me on this topic, it would really help.

Thank you again for your time.

Answer 2

[Philip Caplan]

I have attached some sample code which solves the 2D heat equation (with unit diffusivity) using the "ode45" approach. The solution procedure is the same for your case but you will need to treat the boundary conditions appropriately.

Depending on the level of accuracy you seek, there are various ways of implementing the boundary conditions. You might want to try a first-order scheme for your Neumann conditions. In terms of implementation, I recommend updating the boundary points (using the scheme you choose) before computing du/dt for interior points. You would then leave du/dt=0 for the boundary points (but at least these get updated before computing the interior point derivatives). I've put a comment where this should be done.

Answer 3

[Bill Greene]

Hi Matt,

You can, in fact, express your equation in a form that the parabolic function in PDE Toolbox will accept.

The main "trick" is to recognize that the first derivative terms in x and y (ux and uy terms) have to be moved to the right-hand side and included in the parabolic f-coefficient-- i.e. the f-coefficient is a function of ux and uy. The page you reference: http://www.mathworks.fr/fr/help/pde/ug/scalar-coefficients-in-function-form.html shows how to calculate ux and uy (see section "Gradient or Derivatives of u") in the function you write for the f-coefficient.

Regarding your question about the c-coefficient, for your single-equation (single dependent variable) case, the c-coefficient can be thought of as a 2 x 2 matrix, [c11 c12; c21 c22]. The simplest way to see what the entries are for your equation is to expand the term div(c grad(u)) in the scalar parabolic pde. Because you have only a uxx second derivative term, only the c11 term is non-zero. I am unclear on what "K1(cst)" means but, in the most general case, you might want to define a function for the c-coefficient. It would look something like this:

function c = cfunc(p,t,u,time)
  number_of_elements = size(t, 2);
  cst = ???
  % define the diagonal terms of the c-coefficient matrix at
  % the centroid of each element (this code assumes the c-coefficient
  % doesn't vary as a function of x and y)
  c = repmat([K1(cst) 0]', 1, number_of_elements);
end

I should say one more thing about using the parabolic function to solve this equation. If the ux and uy terms are large compared with the uxx term, the parabolic function may run into numerical difficulties obtaining a solution.

Bill