# Problem with stiffness suspect and local spatial error when using NDSolve

I'm actually having a problem to solve this PDE, in $$0\le x \le 1$$ and $$t>0$$:

$$\frac{{\partial u}}{{\partial t}} + \tau \frac{{{\partial ^2}u}}{{\partial {t^2}}} + \frac{{\partial f(u)}}{{\partial x}} = \varepsilon \frac{\partial }{{\partial x}}\left( {v(u)\frac{{\partial u}}{{\partial x}}} \right)$$

with the boundary and initial conditions:

$$t = 0 \Rightarrow u(x,t) = g(x);$$

$$t = 0 \Rightarrow \frac{{\partial u}}{{\partial t}} = 0$$

$$x = 0 \Rightarrow u(x,t) = 0$$

$$x = 1 \Rightarrow u(x,t) = 1$$

Where:

$$v(u) := 4u(1 - u)$$

$$g(x) := \frac{1}{{1 + {e^{ - (x - a)/\varphi }}}}$$

$$f(u) := \frac{{{u^2}}}{{{u^2} + {{(1 - u)}^2}}}[1 - 5{(1 - u)^2}]$$

And:

$$\varepsilon=0.01$$;

$$a=1-1/2^{1/2}$$;

$$\varphi=0.01$$;

$$\tau=0.00001$$ and $$\tau=0.1$$

My following code are:

a = N[1 - 1/Sqrt[2]];
phi = 0.01;
epsilon = 0.01;
tau = {0.00001, 0.1};

vu = 4*u[x, t]*(1 - u[x, t]);
fu = u[x, t]^2/(u[x, t]^2 + (1 - u[x, t])^2)*(1 - 5*(1 - u[x, t]^2));
gx = 1/(1 + Exp[-(x - a)/phi]);

PDE1 = D[u[x, t], t] + tau[[1]]*D[u[x, t], t, t] + D[fu, x] ==
epsilon*D[(vu*D[u[x, t], x]), x];
PDE2 = D[u[x, t], t] + tau[[2]]*D[u[x, t], t, t] + D[fu, x] ==
epsilon*D[(vu*D[u[x, t], x]), x];
IC = {
u[x, 0] == gx,
(D[u[x, t], t] /. t -> 0) == 0
};
BC = {
u[0, t] == 0,
u[1, t] == 1
};

solND1 = NDSolve[{PDE1, IC, BC}, u, {x, 0, 1}, {t, 0, 10}];
solND2 = NDSolve[{PDE2, IC, BC}, u, {x, 0, 1}, {t, 0, 10}];


Which gives me the following errors:

If anyone can help me find the fix for this error, I would be very grateful.

Grateful for the attention.

• What do you try to simulate with this model? Commented Jul 5, 2022 at 9:34
• Look at the initial condition: u[x, 0]. u[0,0] gives: 1/(1 + Exp[-(- a)/phi]). However, the boundary condition specifies: u[0, 0] == 0. Similar for t==1. Commented Jul 5, 2022 at 10:38

First example with tau[[1]] can be solved with using method of lines and DAE solver as follows

a = 1 - 1/Sqrt[2];
phi = 1/100;
epsilon = 1/100;
tau = {10^-5, 10^-1};

vu = 4*u[x, t]*(1 - u[x, t]);
fu = u[x, t]^2/(u[x, t]^2 + (1 - u[x, t])^2)*(1 - 5*(1 - u[x, t]^2));
gx = 1/(1 + Exp[-(x - a)/phi]);
eqn[i_] :=
D[u[x, t], t] + tau[[i]]*D[u[x, t], t, t] + D[fu, x] ==
epsilon*D[(vu*D[u[x, t], x]), x]; ic = {u[x, 0] ==
gx, (D[u[x, t], t] /. t -> 0) == 0};
bc = {u[0, t] == 0, u[1, t] == 1};
Monitor[uu =
NDSolveValue[{eqn[1], ic, bc}, u, {x, 0, 1}, {t, 0, 1},
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 1370, "MaxPoints" -> 1370,
"DifferenceOrder" -> 2}}},
EvaluationMonitor :> (monitor =
Row[{"t=", CForm[t], " csol=", CForm[u[.5, t]]}])], monitor];


Visualization

{DensityPlot[uu[x, t], {x, 0, 1}, {t, 0, .5}, ColorFunction -> Hue,
PlotRange -> All, PlotPoints -> 100, PlotLegends -> Automatic],
Plot[Table[uu[x, t], {t, 0, .5, .05}] // Evaluate, {x, 0, 1},
PlotRange -> All]}


Second example with tau[[2]] is unstable one.