2

I am working on an HJB recently and wish to get a numerical result. Can NDSolve on Mathematica 12.2 handle this?

$$ V_t+(e^{-\sqrt{t}}+0.4e^{-2.5t}+0.01x)V_{x}+1.12e^{-2t}V_{xx}=\frac{1}{2}e^{-3t}\frac{V_x^2}{V_{xx}} $$

with boundary condition $V(t,0)=1$ and $V(t,\infty)=0$ where $V_t=\frac{\partial}{\partial t}V(t,x)$, $V_x=\frac{\partial}{\partial x}V(t,x)$ and $V_{xx}=\frac{\partial^2}{\partial x^2}V(t,x)$

I tried it on Mathematica with NDSolve and it suggests me to use Inactive but in this case I don't know how to. So I wonder here if Mathematica can actually handle with it. If can, how should I rewrite the equation?

Thanks for any help.

NEW EDIT 5/27 1:40 GMT+8:

Code I used:

prob = {D[v[t, x], t] + (E^(-Sqrt[t]) + 0.4*E^(-2.5*t) + 0.01*x)*
       D[v[t, x], x] + 1.12*E^(-2*t)*D[v[t, x], x, x] - 
       0.5*E^(-3*t)*(D[v[t, x], x])^2/D[v[t, x], x, x] == 0, 
       v[t, 0] == 1, v[t, 1000] == 0}
NDSolve[prob, v, {t, 0, 20}, {x, 0, 20}]

And it turned out:


***Coefficient List: -((0.5 E^(-3 t) v$3459^2)/v$3460)+1.12 E^(-2 t) v$3460+v$3461+v$3459 \
(E^-Sqrt[t]+0.4 E^(-2.5 t)+0.01 x) is not a polynomial.

***NDSolve: The maximum derivative order of the nonlinear PDE coefficients 
for Finite Element Method is larger than 1. It may help to rewrite the PDE in 
inactive form.

That is all I got. The original text was not in English so I translated according to Mathematica documents.

Orsyke
  • 21
  • 2
  • 1
    Can you post the code you used, so we can check what went wrong, what the issues are, etc? –  May 26 '21 at 17:31
  • 2
    Do you have initial conditions? – user21 May 26 '21 at 17:32
  • Question edited with code used and resulted error. @DiSp0sablE_H3r0 – Orsyke May 26 '21 at 17:50
  • No initial conditions with t. In fact this is a simplified version with no restrictions on t. If this problem can be solved, then I may turn to the situation where another boundary condition $V(T,x)=1-sgn(x)$ where $T$ is a given moment and $sgn(.)$ returns the sign of x. @user21 – Orsyke May 26 '21 at 17:53
  • 1
    "No initial conditions with t. In fact this is a simplified version with no restrictions on t. " This doesn't simplify the problem at all, without the initial conditon, the particular solution cannot be determined, please add it to NDSolve. Also, please add some background info if possible. – xzczd Jun 26 '21 at 00:55

2 Answers2

4

It is not a classical quasilinear PDE, but it can be solved as optimization problem with using Euler wavelets and method, described in our paper as follows

u0[x_] := 
 Exp[-2 x]; L = 5; T = 1; prob = {D[v[t, x], 
     t] + (E^(-Sqrt[t]) + 0.4*E^(-2.5*t) + 0.01*x)*D[v[t, x], x] + 
    1.12*E^(-2*t)*D[v[t, x], x, x] - 
    0.5*E^(-3*t)*(D[v[t, x], x])^2/D[v[t, x], x, x] == 0, 
  v[t, 0] == 1, v[t, 1000] == 0};

UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] := 
 Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
      2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] := 
 Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 2; M0 = 4; nn = 
 Total[With[{k = k0, M = M0}, 
   Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn);  xl = Table[ l*dx, {l, 0, nn}]; tcol = 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1]; Psi[y_] := Psijk /. t1 -> y; 
int1[y_] := Int1 /. t1 -> y; int2[y_] := Int2 /. t1 -> y;
M = nn;


U1 = Array[a, {M, M}]; U2 = Array[b, {M, M}]; G1 = 
 Array[g1, {M}]; G2 = Array[g2, {M}]; G3 = Array[g3, {M}]; F1 = 
 Array[f1, {M}]; F2 = Array[f2, {M}]; F3 = Array[f3, {M}];


uxxt[t_, x_] := Psi[t] . U1 . Psi[x]; 
uxx[t_, x_] := int1[t] . U1 . Psi[x] + G1 . Psi[x]; 
ux[t_, x_] := int1[t] . U1 . int1[x] + G1 . int1[x] + G2 . Psi[t]; 
u1[t_, x_] := 
 int1[t] . U1 . int2[x] + G1 . int2[x] + x G2 . Psi[t] + G3 . Psi[t];
uxt[t_, x_] := Psi[t] . U1 . int1[x] + F1 . Psi[t];
ut[t_, x_] := Psi[t] . U1 . int2[x] + x F1 . Psi[t] + F2 . Psi[t]; 
u2[t_, x_] := 
 int1[t] . U1 . int2[x] + x F1 . int1[t] + F2 . int1[t] + F3 . Psi[x];


eq = Join[
  Flatten[Table[
    ut[tcol[[i]], xcol[[j]]]/
      T + (E^(-Sqrt[Ttcol[[i]]]) + 0.4E^(-2.5T tcol[[i]]) + 
        L 0.01xcol[[j]]) ux[tcol[[i]], xcol[[j]]]/L + 
     1.12E^(-2T tcol[[i]]) uxx[tcol[[i]], xcol[[j]]]/L^2 - 
     0.5E^(-3T tcol[[i]]) ux[tcol[[i]], xcol[[j]]]^2/
       uxx[tcol[[i]], xcol[[j]]], {i, M}, {j, M}]], 
  Flatten[Table[
    u1[tcol[[i]], xcol[[j]]] - u2[tcol[[i]], xcol[[j]]], {i, M}, {j, 
     M}]]]; bc = 
 Join[Table[u1[tcol[[i]], 0] == 1, {i, M}], 
  Table[u1[tcol[[i]], 1] == 0, {i, M}], 
  Table[u2[0, xcol[[i]]] == u0[L xcol[[i]]], {i, M}], 
  Table[u2[tcol[[i]], 0] == 1, {i, M}], 
  Table[u2[tcol[[i]], 1] == 0, {i, M}], 
  Table[u1[0, xcol[[i]]] == u0[L xcol[[i]]], {i, M}]];


var = Join[Flatten[U1], G1, G2, G3, F1, F2, F3];


sol1 = NMinimize[{eq . eq, bc}, var, MaxIterations -> 200];

Visualization

u = Interpolation[
  Flatten[Table[{t, x, u2[t/T, x/L] /. sol1[[2]]}, {t, 0, T, 
     T/10}, {x, 0, L, L/10}], 1]]

Plot3D[u[t, x], {t, 0, T}, {x, 0, L}, ColorFunction -> "Rainbow", 
 MeshStyle -> White, PlotRange -> All, AxesLabel -> Automatic]

Figure 1

Alex Trounev
  • 44,369
  • 3
  • 48
  • 106
1

Your equation is not written in a form that is solvable by Mathematica. If it can be rewritten in such form, please do so to get an answer. My personal guess is that it cannot, as it contains $V_{xx}^2$ term.

See NDSolve error: what does "It may help to rewrite the PDE in inactive form" mean?

Giorgio Busoni
  • 409
  • 2
  • 7
  • 1
    All right... Thanks a lot~~~ So if FEM cannot work, how about finite difference method in Mathematica 12.2? – Orsyke May 27 '21 at 02:12
  • I don't know, but I tend to doubt you can solve an equation in this form with it. On top of that, I think your boundary conditions are incompatible with the equation, at least to get a $\mathcal{C}^2$ solution. At $t\rightarrow\infty$, the solution should converge to the static solution of the equation, that is $xV_{x}V_{xx}=0$. Now, the only solution compatible with your initial conditions is $V=0$ if $x>0$, $V=1$ if $x=0$, that is not $\mathcal{C}^0$ in $x=0$. – Giorgio Busoni May 28 '21 at 10:43
  • ok here is what you can do. 1. Choose a starting function $f(x)$ such that $f(0)=1, f(\infty)=0, f_x(0)=0, f_{xx}(0)=0, f_x(\infty)=0, f_{xx}(\infty)=0$. 2 discretise the time evolution in steps $dt$, choose $dt$ 3. Find $V_t(x,0)$ setting $V(x,0)=f(x)$ 4. Find $V(x,dt)=V_(x,0)+dt*V_t(x,0)$ 4.Repeat unlimited times until you get to a $t$ as large as you need. t=25 should already be a good approximation 5. You can evolve the solution for larger $t$ by using the approximate equation $V_t + 0.01xV_x=0$ (you set exponentials to zero). This is linear, so it can be handled easily – Giorgio Busoni May 28 '21 at 11:24
  • This procedure seem to handle a solution for any starting $f(x)$ that satisfies the above conditions, so it seems tou need to set an additional condition there to get an unique solution – Giorgio Busoni May 28 '21 at 11:25
  • Thanks a lot! I'll try this iteration method and give feedback later. @Giorgio Busoni – Orsyke May 28 '21 at 11:35
  • I was trying to get a working example, but it seems that the boundary condition, evolving in $t$, stops being satisfied and departs from the chosen boundary solution after a few steps. The fact is that at each step one wants $V_x(t,x)=0, V_{xx}(t,x)=0$, and that this implies, after n steps $D[V,{x,2n}]=0$. So to keep the boundary condition satisfied you'd need an $f$ with all derivaties in 0 equal to 0. Seems a very irregular solution – Giorgio Busoni May 28 '21 at 12:39
  • Yes, I got into the same trap too. Maybe no solution indeed exist for this equation... Now I’m working on the existence of solution of this pde. – Orsyke May 29 '21 at 04:34
  • I think your best chance is to ask on the mathematics forum for a solution strategy, then you can try to implement it in mathematica using some iterative method like mine. Indeed, maybe for a special choice of the starting function the boundary coindition can remain satified, or maybe no solution exists. I am not expert of nonlinear-PDE, so can't help you with that. – Giorgio Busoni May 29 '21 at 07:28