# Unstable result from numerical solution of second order partial differential equation using NDSolve

I want to solve following differential equation

$$$$\frac{\partial u(x,t)}{\partial t} -\frac{1}{2}\frac{\partial^2 u(x,t)}{\partial t^2} = M\left[\frac{\partial ^2 u[x,t]}{\partial x^2} + \frac{\partial }{\partial x} u[x,t](1-u[x,t])\right]$$$$

subject to initial condition $$u[x,0] = 1 ~~ if~~x and boundary condition $$u[0,t] =1, u[L,t] = 0$$. I have tried following

L = 50; M = 0.2;
g[x_] := Piecewise[{{1, x < L/2}, {0, x >= L/2}}];
pdeSec = D[uSec[x, t], t] - D[D[uSec[x, t], t], t]/2 -
M ( D[D[uSec[x, t], x], x] +
D[uSec[x, t] (1 - uSec[x, t]), x]) == 0;
solSec = NDSolve[{pdeSec, uSec[x, 0] == g[x], uSec[0, t] == 1,
uSec[L, t] == 0}, uSec[x, t], {x, 0, L}, {t, 0, 200},
Method -> StiffnessSwitching];
Plot[{uSec[x, t] /. solSec /. {t -> 20}}, {x, L/2 - 10, L/2 + 10}]


With this I get an error message

InitializePDECoefficients: The product of 2 times the Damkoehler number (0.19611613513818169) and the Peclet number (12.747548783981808) is 4.999999999999879 and is larger than the mesh order (2), and the computed result may not be stable. Adding artificial diffusion may help


Reading this message, I tried to change the coefficient $$M$$ to increase/decrease the diffusion term. The solution thus obtained, nevertheless, remains unstable. Any help to obtain a smooth, stable solution will be greatly appreciated.

• Add option to NDSolve: Method -> {"FiniteElement", "InterpolationOrder" -> {uSec -> 2}, "MeshOptions" -> {"MaxCellMeasure" -> 0.1}}.I tried on MMA 12.2.0 ,works fine. Commented Jan 31, 2021 at 9:26
• thanks a lot @MariuszIwaniuk, specifying the suggested method works as intented ..! If you could write it as an answer, I will be happy to accept it as solution. Commented Jan 31, 2021 at 10:04
• If you want to solve this as a time dependent problem at the very least you must specify a derivative for the initial condition because you have a second order in time equation. Something like Derivative[0, 1][uSec][x, 0] ==.... Commented Feb 1, 2021 at 6:35

Try:

L = 50; M = 0.2;
g[x_] := Piecewise[{{1, x < L/2}, {0, x >= L/2}}];
pdeSec = D[uSec[x, t], t] - D[D[uSec[x, t], t], t]/2 -
M ( D[D[uSec[x, t], x], x] +
D[uSec[x, t] (1 - uSec[x, t]), x]) == 0;
solSec = NDSolve[{pdeSec, uSec[x, 0] == g[x], uSec[0, t] == 1,
uSec[L, t] == 0}, uSec[x, t], {x, 0, L}, {t, 0, 200},
Method -> {"FiniteElement","MeshOptions" -> {"MaxCellMeasure" -> 0.1}}];
Plot[{uSec[x, t] /. solSec /. {t -> 20}}, {x, L/2 - 10, L/2 + 10}]
`

• You never need to specify "InterpolationOrder" for single equations. "InterpolationOrder" is only ever relevant for coupled equations. Where did you find the usage of this option? I'd like to update the documentation to make this explicit. I am going ahead and remove this option from yur code to prevent further spread of this. Commented Feb 1, 2021 at 6:31
• One other thing to mention is that you solve this as a spatial problem. (which in this case is fine +1) but the intend probably was to solve this as a time dependent problem. For that one would need to add a derivative of the initial condition. Commented Feb 1, 2021 at 6:33
• @user21 .Ok thanks for update. :) Commented Feb 1, 2021 at 7:51