# Numerical solution of a nonlinear PDE that develops a growing piecewise linear region

I am trying to improve the numerical solution of some PDEs that develop a piecewise behavior during their evolution. The simplest example of one such PDE is for a function $$u(t,x)$$ with $$t \in [-T,T]$$ and $$x \in \mathbb{R}$$: $$u_t(t,x) = \frac{\sqrt{T^2-t^2}}{2\pi T^2} \frac{t^2 f(x) u_{xx}(t,x)}{1- t u_x(t,x)},$$ with initial condition $$u(-T,x) = f(x)$$ and boundary conditions $$u(t,x) \sim f(x)$$ for large $$|x|$$, implemented numerically by choosing a value $$L$$ and imposing $$u(t,\pm L) = f(\pm L)$$. For the examples here I will focus on $$f(x) = {\rm LogisticSigmoid}(x) = 1/(1+\exp(-x))$$.

The behavior of the solution depends on the value of the parameter $$T$$, which sets the initial and final "time". For sufficiently small $$T$$ the solution remains smooth, but for sufficiently large $$T$$ I expect there will be a time $$t_\ast$$ and point $$x_\ast$$ such that $$1 - t_\ast u_x(t_\ast,x_\ast) = 0$$--i.e., the denominator vanishes, but I also expect this will be compensated by $$u_{xx}(t_\ast,x_\ast) = 0$$. The solution can continue for $$t_\ast < t$$ if it becomes piecewise linear, such that there is a range $$x^-(t) < x < x^+(t)$$ for which $$1 - t u_x(t,x) = u_{xx}(t,x) = 0$$, and otherwise obeys the above PDE when $$x$$ is outside this range.

I have solved the above PDE numerically using two methods. The first is NDSolve using MethodOfLines with a TensorProductGrid discretization:

mineg[T_, L_] :=
NDSolve[{D[u[t, x], t] ==
Sqrt[T^2 - t^2]/(2*Pi*T^2)*t^2*LogisticSigmoid[x]*
D[u[t, x], x, x]/(1 - t*D[u[t, x], x]),
u[-T, x] == LogisticSigmoid[x], u[t, -L] == LogisticSigmoid[-L],
u[t, L] == LogisticSigmoid[L]}, u, {t, -T, T}, {x, -L, L},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> {600}}}}]

egT10L40 = mineg[10, 40]

(* Plot solution versus initial condition *)
Manipulate[
Plot[{LogisticSigmoid[x], Evaluate[u[t, x] /. egT10L40]}, {x, -10,
10}, PlotRange -> {0, 1}], {t, -10, 10}]

(* Plot denominator to see that it runs up against zero *)
Manipulate[
Plot[{1 - t*LogisticSigmoid'[x],
Evaluate[1 - t*Derivative[0, 1][u][t, x] /. egT10L40]}, {x, -10,
10}, PlotRange -> {-4, 4}], {t, -10, 10}]

(* Plot second derivative *)
Manipulate[
Plot[{LogisticSigmoid''[x],
Evaluate[Derivative[0, 2][u][t, x] /. egT10L40]}, {x, -10, 10},
PlotRange -> {-0.25, 0.25}], {t, -10, 10}]


You can see that as $$t$$ runs from $$-T$$ to $$T$$ the solution eventually reaches a point where the denominator appears to touch zero, and then continues to be pressed against zero instead of becoming negative, and in this region the second derivative is also close to zero. The behavior gets a bit noisy as $$t$$ approaches $$T$$, and while this qualitatively agrees with the picture I claimed above, I am not sure if there is any artifact coming from the numerator and denominator both approaching zero numerically.

Solution at $$t = T = 10$$, compared to initial condition.

Denominator at $$t = 10$$ (compared to a fictitious denominator with the initial condition plugged in)

Second derivative at $$t = 10$$. Note the sharp peak near $$x = -5$$.

I attempted to improve this numerical solution by using the pdetoode tool developed by @xzczd:

Tp = 10; Lp = 40;
eqn = D[u[t, x], t] ==
Sqrt[Tp^2 - t^2]/(2*Pi*Tp^2)*t^2*LogisticSigmoid[x]*
D[u[t, x], x, x]/(1 - t*D[u[t, x], x]);
ic = {u[-Tp, x] == LogisticSigmoid[x]};
bc = {u[t, -Lp] == LogisticSigmoid[-Lp],
u[t, Lp] == LogisticSigmoid[Lp]};

(*Difference order of x:*)
xdifforder = 4;

points = 500;
grid = Array[# &, points, {-Lp, Lp}];

(*There're 2 b.c.s,so we need to remove 2 equations from every \
PDE/i.c.,usually the difference equations that are the "closest" ones \
to the b.c.s are to be removed:*)
removeredundant = #[[2 ;; -2]] &;
(*Use pdetoode to generate a "function" that discretizes the spatial \
derivatives of PDE(s) and corresponding i.c.(s) and b.c.(s):*)

ptoofunc = pdetoode[u[t, x], t, grid, xdifforder];

odeqn = eqn // ptoofunc // removeredundant;
odeic = removeredundant /@ ptoofunc@ic;
odebc = bc // ptoofunc;

sollst = NDSolveValue[{odebc, odeic, odeqn}, u /@ grid, {t, -Tp, Tp},
MaxSteps -> Infinity];
(*Rebuild the solution for the PDE from the solution for the ODE set:*)

sol = rebuild[sollst, grid];

(* Solution compared to initial condition *)
Manipulate[
Plot[{LogisticSigmoid[x], sol[t, x]}, {x, -10, 10},
PlotRange -> All], {t, -Tp, Tp}]

(* Denominators *)
Manipulate[
Plot[{1 - t*LogisticSigmoid'[x],
1 - t*Derivative[0, 1][sol][t, x]}, {x, -10, 10},
PlotRange -> {-4, 4}], {t, -Tp, Tp}]

(* Second derivatives *)
Manipulate[
Plot[{LogisticSigmoid''[x], Derivative[0, 2][sol][t, x]}, {x, -10,
10}, PlotRange -> {-0.25, 0.25}], {t, -Tp, Tp}]


Solution at $$t = 10$$

Denominator:

Second derivative:

This appears to improve the stability of the solution, can go to larger values of $$T$$, and maintains the qualitative picture I described where it appears that the solution is becoming close to piecewise linear in the center. However, it does not appear to have as flat a region as I would expect, and attempts to find zero crossings numerically only return a single value. This value also does not coincide with the point where the solution appears to abruptly become exactly equal to the initial condition for negative $$x$$ (roughly near -4 or -5).

Perhaps this is a limitation of the method of solution, but I am wondering if the solution can be improved further to actually implement the piecewise linear behavior by identifying when the denominator vanishes and then imposing that at later times the solution is linear between edges $$x^-(t)$$ and $$x^+(t)$$ that move in time and must also be tracked. This sounds superficially similar to problems on this site solving PDEs with Stefan boundary conditions, though these boundaries appear inside the solution at some intermediate time $$t_\ast$$.

Specific issues/questions:

1. I imagine that to do this I would need to first identify when ($$t_\ast$$) and where ($$x_\ast$$) the denominator vanishes, and then track the evolution of the boundary points for times $$t_\ast < t$$. I know WhenEvent does not work on bivariate functions; I imagine something can be done with the pdetoode discretization, but I am not sure how to implement it.

2. I am not sure how to derive equations that dictate the evolution of the piecewise linear region boundary points $$x^\pm(t)$$. I would have guessed that I would define these points by $$1 - t u_x(t,x^\pm(t)) = 0$$ and differentiate this relationship, but it did not appear to give me something useful. This also does not seem to be how the boundary conditions are derived in the Stefan boundary problems I saw in other questions (e.g., here.). I am not clear on whether the moving boundaries in those questions are just given as part of the model/problem statement, or if there is a way to derive differential equations for moving boundaries from the governing differential equation for the full solution and the definition of the boundary.

Any help figuring out how to improve the numerical solution will be much appreciated.

• Why not to post this question on math.stackexchange.com ? Commented Jul 29, 2023 at 1:52

With the following modification, I find a solution:

1. Densify the grid to points = 1000;

2. Modify definition of odeic and odebc to

odeic = ptoofunc@ic;
odebc = diffbc[t, 1]@bc // ptoofunc;


I make this modification because empirically this makes the system less stiff.

3. The most important part: move the denominator (1 - t D[u[t, x], x]) to left hand side of the PDE, and add a tiny extra term -D[t D[u[t, x], x], t] eps to the system. I find this by accident so cannot give a theoretical explanation, but this just relieve the stiffness of the system.

The following is the complete code sample:

showStatus[status_]:=LinkWrite[$ParentLink, SystemSetNotebookStatusLine[FrontEndEvaluationNotebook[], ToString[status]]]; clearStatus[]:=showStatus[""]; clearStatus[] jianshi[t_]:=EvaluationMonitor:>showStatus["t = "<>ToString[CForm[t]]] mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, "MinPoints" -> n, "DifferenceOrder" -> o}} Tp = 10; Lp = 40; eps = 10^-4; eqn = -D[t D[u[t, x], x], t] eps + (1 - t D[u[t, x], x]) D[u[t, x], t] == Sqrt[Tp^2 - t^2] t^2 LogisticSigmoid[x] D[u[t, x], x, x]/(2 Pi Tp^2); ic = {u[-Tp, x] == LogisticSigmoid[x]}; bc = {u[t, -Lp] == LogisticSigmoid[-Lp], u[t, Lp] == LogisticSigmoid[Lp]}; xdifforder = 4; points = 1000; grid = Array[# &, points, {-Lp, Lp}]; removeredundant = #[[2 ;; -2]] &; ptoofunc = pdetoode[u[t, x], t, grid, xdifforder]; odeqn = eqn // ptoofunc // removeredundant; odeic = ptoofunc@ic; odebc = diffbc[t, 1]@bc // ptoofunc; sollst = NDSolveValue[{odebc, odeic, odeqn}, u /@ grid, {t, -Tp, Tp}, MaxSteps -> Infinity, jianshi[t], SolveDelayed -> True]; sol = rebuild[sollst, grid]; Manipulate[Plot[sol[t, x], {x, -40, 40}], {t, -Tp, Tp}]  Adjusting eps to 10^-3, the solution is visually the same, which seems to suggest the solution is reliable. • Thank you for your experiments (and the pdetoode function itself!). Your change does appear to lead to a flatter region like I was hoping for, albeit with some oscillations in the flat region of the derivatives (that get worse if I try to increase the number of points further). I was originally avoiding multiplying through by the denominator because in some related systems of PDEs the denominators involve higher powers of this same denominator (meaning I would end up with LHSs$(1-t u_x) u_t$,$(1-tu_x)v_t\$, etc, so I will see how this sort of trick works in these more complicated examples. Commented Aug 4, 2023 at 17:56

For theory please see, for example, Parabolic Equations with Changing Direction of Time

Solution can be extended to t=T with small modification of the code by multiplying both parts on (1 - (t*D[u[t, x], x])^2) then we have

mineg[T_, L_] :=
NDSolve[{(1 - (t*D[u[t, x], x])^2) D[u[t, x], t] ==
Sqrt[T^2 - t^2]/(2*Pi*T^2)*t^2*LogisticSigmoid[x]*
D[u[t, x], x, x] (1 + t*D[u[t, x], x]),
u[-T, x] == LogisticSigmoid[x], u[t, -L] == LogisticSigmoid[-L],
u[t, L] == LogisticSigmoid[L]}, u, {t, -T, T}, {x, -L, L},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> {400}}}}]

egT10L40 = mineg[10, 40]

(*Plot solution versus initial condition*)
Manipulate[
Plot[{LogisticSigmoid[x], Evaluate[u[t, x] /. egT10L40]}, {x, -10,
10}, PlotRange -> {0, 1}], {t, -10, 10}]

(*Plot denominator to see that it runs up against zero*)
Manipulate[
Plot[{1 - t*LogisticSigmoid'[x],
Evaluate[1 - t*Derivative[0, 1][u][t, x] /. egT10L40]}, {x, -10,
10}, PlotRange -> {-4, 4}], {t, -10, 10}]

(*Plot second derivative*)
Manipulate[
Plot[{LogisticSigmoid''[x],
Evaluate[Derivative[0, 2][u][t, x] /. egT10L40]}, {x, -10, 10},
PlotRange -> {-0.25, 0.25}], {t, -10, 10}]


Visualization at t=T

• Wow, this is surprising. (Of course +1) Commented Jul 29, 2023 at 5:13
• @xzczd Thank you. I first derived this type of equation in the paper link.springer.com/article/10.1007/BF00919627 :) Commented Jul 29, 2023 at 6:27
• Thanks for the references, and the suggestions to increase points and multiply by (1- (t u_x)^2). It seems that increasing MinPoints to 600 in my original code is sufficient to extend the solution to T = 10 -- the difference between that and the proposed modification here appear only in the 2nd derivative, and vary only by 0.01 or so. But, I'll explore the idea on some of my more complicated examples involving systems of PDEs to see if that helps at all there. Commented Aug 4, 2023 at 17:46
• If you have any other references on these sorts of singular parabolic PDEs, I would also appreciate them! Commented Aug 4, 2023 at 17:48
• See, for example, M. Gevrey, Sur les equations aux derivees partielles du type parabolique,J. Math. Pures Appl.,9,10, 305–475 105–148 (1913, 1914). M. S. Baouendi and P. Grisvord, Sur une equation d'evolution changeant de type,J. Funct. Anal.,2, 352–367 (1968). C. D. Pagani and G. Talenti, On a forward-backward parabolic equation, Ann. Math. Pura Appl.,90, 1–58 (1971). Commented Aug 5, 2023 at 7:13