# Solution of convection dominated PDE with conditions only on the interior of the domain

I'm Trying to solve the following PDE:

$$\frac{1}{20}\left(2 x_1^2+x_2^2\right) -25\left( \begin{bmatrix} 0\\ \frac{6 \cos(x_1)}{-37 + 3 \cos(2 x_1)} \end{bmatrix}^\intercal \nabla_{x_1,x_2}V(x_1,x_2) \right)^2 + \begin{bmatrix} x_2\\ \frac{6 (-98 + x_2^2 \cos(x_1)) \sin(x_1)}{-37 + 3 \cos(2 x_1)} \end{bmatrix}^\intercal \nabla_{x_1,x_2}V(x_1,x_2) = 0$$

The only condition I have is $$V(0,0)=0$$ and I am looking for a positive definite solution. To give some context, this comes from a Hamilton-Jacobi-Bellman PDE.

ClearAll[RHS]
RHS = 1/20 (2 x1^2 + x2^2) - 25({0, (6 Cos[x1])/(-37 + 3 Cos[2 x1])}.Inactive[Grad][V[x1, x2], {x1, x2}])^2 + {x2,( 6 (-98 + x2^2 Cos[x1]) Sin[x1])/(-37 + 3 Cos[2 x1])}.Inactive[Grad][V[x1, x2], {x1, x2}]

ClearAll[ufun, mesh]
Needs["NDSolveFEM"]
mesh = ToElementMesh[Rectangle[{0, 0}, {1, 1}], "MaxBoundaryCellMeasure" -> 0.005, "MeshElementType" -> TriangleElement];

ufun = NDSolveValue[{RHS == 0, DirichletCondition[V[x1, x2] == 0, x1 == 0 && x2 == 0]}, V, {x1, x2} \[Element] mesh]

Plot3D[ufun[x1, x2], {x1, 0, 1}, {x2, 0, 1}, ColorFunction -> "TemperatureMap", AxesLabel -> Automatic]


Which gives some interpolating function with the warning: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help. I am looking for a solution on a rather larger domain (e.g. Rectangle[{-10, -10}, {10, 10}), but I am finding several problems. First, I am not sure on how to enforce the condition $$V(0,0)=0$$ if (0,0) is not in the boundary of the region. Second, I am not sure on how to enforce the fact that V should be positive on the region and finally, I am not sure how to get rid of that warning about the PDE being convection dominated.

Any help would be greatly appreciated.

• Is it possible to rearrange the problem in the form of the objective functional and the transition function(or variational form)? it's much easier to understand the problem. – Xminer Dec 30 '19 at 18:29
• Sure, it is an inverted pendulum on a cart with objective: 1/10 x1^2 + 1/20 x2^2 + 1/100 u^2. The dynamics of the pendulum are standard, (two states, angles measured from the vertical, length 1/2, mass of the pendulum 2 and mass of the cart 8. The cart dynamics are neglected and the cart acceleration taken as the control. The right-hand side of the ODE can be reconstructed as $\dot x = f(x(t)) + g(x(t)) u(t)$ with $f(x) = {x_2,( 6 (-98 + x_2^2 \cos(x_2)) \sin(x_1))/(-37 + 3 \cos(2 x_1))}$ (the vectorfield multiplying the second gradient) and $g(x)={0, (6 \cos(x_1))/(-37 + 3 \cos(2 x_1))}$ – NoobNoob Dec 30 '19 at 18:49
• Was that what you meant? – NoobNoob Dec 30 '19 at 19:15
• yes, that forms make it easier to understand the problem :) – Xminer Dec 30 '19 at 21:32

This is a nonlinear equation. To be able to use FEM you have to get the equation into the coefficient form: Since the form you present is not in this coefficient form strange things happen. The easiest way around this is probably to use Activate and have the FEM solver figure it out - which seems to work in this case but is not a general cure.

ClearAll[ufun, mesh]
Needs["NDSolveFEM"]
mesh = ToElementMesh[Rectangle[{0, 0}, {1, 1}],
"MaxBoundaryCellMeasure" -> 0.005,
"MeshElementType" -> TriangleElement];
Block[{e = 1, s = 1, j = 1},
ufun = NDSolveValue[{Activate[RHS] == 0,
DirichletCondition[V[x1, x2] == 0, x1 == 0 && x2 == 0]},
V, {x1, x2} \[Element] mesh,
Method -> {"PDEDiscretization" -> {"FiniteElement",
"PDESolveOptions" -> {"FindRootOptions" -> {EvaluationMonitor \
:> e++, StepMonitor :> s++,
Jacobian -> {Automatic, EvaluationMonitor :> j++}}}}}];
{"Steps" -> s, "Function Evaluations" -> e,
"Jacobian Evaluations" -> j}]


This gives messages

InitializePDECoefficients::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

InitializePDECoefficients::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

FindRoot::dfmin: The minimal damping factor of 1/10000 has been reached.

{"Steps" -> 6, "Function Evaluations" -> 9,
"Jacobian Evaluations" -> 5}


Yes, the PDE is convection dominated. This is a warning. Not all convection dominated PDEs are unstable. But the bigger problem is that the nonlinear solver can not find a solution. You could try to help it with specifying an InitialSeeding. The ones I tried did not work, but you may know more about the problem.

• Thanks, I will try this. I was wondering tho if this would help with the fact that I need a solution on Rectangle[{a,b},{c,d}] with zero in the interior of the domain. I don't have my mathematica at hand, but I remember that the solver was complaining that V[0,0]=0 needed to be specified on the boundary. Also, If you could quickly comment on what sort of Initial seeding is needed. Thanks!. – NoobNoob Jan 2 at 5:00
• Do you impose zero with a DirichletConditon? Make sure the ElementMesh has a node there by using the "IncludePoints" option to ToElementMesh. To see how InitialsSeeding is to be used have a look at the ref page. – user21 Jan 2 at 7:39