3
$\begingroup$

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["NDSolve`FEM`"]
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.

Improve this question
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Xminer
    Commented Dec 30, 2019 at 18:29
  • $\begingroup$ 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))}$ $\endgroup$
    – NoobNoob
    Commented Dec 30, 2019 at 18:49
  • $\begingroup$ Was that what you meant? $\endgroup$
    – NoobNoob
    Commented Dec 30, 2019 at 19:15
  • $\begingroup$ yes, that forms make it easier to understand the problem :) $\endgroup$
    – Xminer
    Commented Dec 30, 2019 at 21:32

1 Answer 1

Reset to default
2
$\begingroup$

This is a nonlinear equation. To be able to use FEM you have to get the equation into the coefficient form:

enter image description here

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["NDSolve`FEM`"]
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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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!. $\endgroup$
    – NoobNoob
    Commented Jan 2, 2020 at 5:00
  • $\begingroup$ 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. $\endgroup$
    – user21
    Commented Jan 2, 2020 at 7: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.