# Why does smooth initial condition involving Piecewise cause mxsst warning?

This is a very simple one-dimensional heat-conduct equation, the only special part of it is the piecewise initial condition:

b = NDSolve[{D[tes[t, x], t] == D[tes[t, x], x, x] + Exp[-1/tes[t, x]],
tes[t, 0] == 1, tes[t, 1] == 1,
tes[0, x] == Piecewise[{{-100 (x - 0.1)^2 + 2, 0 <= x <= 0.1},
{2, 0.1 <= x <= 0.9},
{-100 (x - 0.9)^2 + 2, 0.9 <= x <= 1}}]},
{tes[t, x]}, {t, 0, 100}, {x, 0, 1}]


If you run the code, you will get this warning message:

NDSolve::mxsst: Using maximum number of grid points 10000 allowed by the MaxPoints or MinStepSize options for independent variable x

Why does this message come out? I read the help of mxsst, it states that

This warning is normally seen only if the initial condition for a partial differential equation has sharp or oscillating features…

but I don't think my initial condition has that kind of fault: it's piecewise but smooth, right?

I'd like to add another sample here since its behavior makes an interesting contrast to the sample above:

c = NDSolve[{D[tes[t, x], t] == D[tes[t, x], x, x] + Exp[-1/tes[t, x]],
tes[t, 0] == 1, tes[t, 1] == 1,
tes[0, x] == Piecewise[{{-Exp[-1000 x] + 2, 0 <= x <= 0.1},
{2, 0.1 <= x <= 0.9},
{-Exp[-1000 (1 - x)] + 2, 0.9 <= x <= 1}}]},
{tes[t, x]}, {t, 0, 100}, {x, 0, 1}]


In this example, Exp is chosen to be the transition between initial condition and boundary conditions. This initial condition is not even continuous at x=0.1 and x=0.9 in the view of math, and the change of Exp is more drastic than the polynomial function in my view, but it causes no warning message. Well, I should say, after all this time (notice the time I posted this question) I've already treated the warning as some kind of bug, but I still expect an in-depth explanation.

Aha, after playing with NDSolve for 10 years, I find the answer of my first question in this site!

I don't think my initial condition has that kind of fault: it's piecewise but smooth, right?

Yes, but prior error estimation of NDSolve doesn't only care about the smoothness of $$T(0,x)$$ (tes[0, x]) in this case. What's concerned is the smoothness of $$\frac{\partial ^2T(0, x)}{\partial x^2} + e^{-1/T(0,x)}$$ (D[tes[0, x], x, x] + Exp[-1/tes[0, x]]) i.e. right hand side (RHS) of the PDE at initial time.

Using higher order polynomial to build an initial condition (i.c.) with continuous second order derivative will resolve the problem, for example:

icexpr = Piecewise[{{2 - 1000 (0.1 - x)^3, 0 <= x <= 0.1},
{2 - 1000 (-0.9 + x)^3, 0.9 <= x <= 1}}, 2]

NDSolve[{D[tes[t, x], t] == D[tes[t, x], x, x] + Exp[-1/tes[t, x]],
tes[t, 0] == 1,
tes[t, 1] == 1, tes[0, x] == icexpr},
{tes[t, x]}, {t, 0, 100}, {x, 0, 1}]
(* Free of warning *)


First of all, let me untangle the sample a bit:

{xL = 0, xR = 1, dx = 0.1, T0 = 2, Tb = 1};

ic[n_][x_] :=
Piecewise[{{-dx^-n (dx - x)^n (T0 - Tb) + T0, xL <= x <= xL + dx},
{-dx^-n (x - (xR - dx))^n (T0 - Tb) + T0, xR - dx <= x <= xR}},
T0];

sample[n_, opt : OptionsPattern[]] :=
With[{T = T[t, x]},
NDSolve[{D[T, t] == D[T, x, x],
T == ic[n][x] /. t -> 0,
T == Tb /. {{x -> xL}, {x -> xR}}},
T, {t, 0, 100}, {x, xL, xR}, opt]]

sample[2]


NDSolve::mxsst

The $$e^{-1/T}$$ term (Exp[-1/tes[t, x]] in the question) has been removed because it turns out to be irrelevant to the problem.

To understand why the mxsst pops up, we need to know how the prior error estimation for PDE is done by NDSolve. The topic is discussed mostly in (but not limited to!) Spatial Error Estimates section of the obscure tutorial The Numerical Method of Lines. As a precondition, we also need to know how ODE(s) is pre-processed by NDSolve. I suggest reading this answer to understand the topic better………… OK, OK, I know it's non-trivial to read through all these, so let me try my best to retell the related parts in a easy-to-follow way.

When solving a time-dependent PDE (to be precise, PDE with independent variable that can be viewed as time), NDSolve will transform the PDE to a standard PDE, whose left hand side (LHS) is 1st order derivative respects to time. (Of course, for the specific problem discussed here, we don't need any further transform, because $$\frac{\partial T}{\partial t}=\frac{\partial ^2T}{\partial x^2}$$ is already in standard form. )

And the prior error estimation is built on the RHS of the obtained standard PDE.

The whole estimation process is a bit involved so I'd like not to talk too much about the details. We just need to remember: the RHS ($$\frac{\partial ^2T}{\partial x^2}$$ in our case) at initial time (t == 0 in our case) will be evaluated based on finite difference method (FDM), but the naive FDM won't be able to approximate a discontinuous function perfectly, because of the Gibbs phenomenon:

num = 101;
grid = If[$VersionNumber < 9, Table[x, {x, xL, xR, (xR - xL)/(num - 1)}], Array[# &, num, {xL, xR}]]; rhsfunc[func_, grid_] := NDSolveFiniteDifferenceDerivative[2, grid, func /@ grid] rhs1 = rhsfunc[ic[2], grid]; Plot[ic[2]''[x], {x, xL, xR}, PlotRange -> All]~Show~ ListPlot[{grid, rhs1}\[Transpose], PlotStyle -> {Red}]  The error caused by Gibbs phenomenon contributes quite a bit to the error estimation of NDSolve. For num = 101, the estimated error is: grid2 = If[$VersionNumber < 9, Table[x, {x, xL, xR, (xR - xL)/(2 num - 2)}],
Array[# &, 2 num - 1, {xL, xR}]];

rhs2 = rhsfunc[ic[2], grid2][[;; ;; 2]];

diff = Module[{lst = {rhs1, rhs2}},
lst[[All, ;; 3]] = 0.; lst[[All, -3 ;;]] = 0.;
Subtract @@ lst];

If[$VersionNumber < 11, ListPlot[{grid, diff}\[Transpose] /. {x_, y_} /; Abs@y > 1 :> Tooltip[{x, y}, TraditionalForm@"x" == N@x], PlotRange -> All] /. Tooltip[{_, pts_}, label_] :> {pts, Text[label, pts[[1, 1]] + {0.07, 0.5}]}, ListPlot[{grid, diff}\[Transpose] /. {x_, y_} /; Abs@y > 1 :> Callout[{x, y}, TraditionalForm@"x" == N@x], PlotRange -> All]]  p = If[$VersionNumber < 9, Infinity, 2];
NDSolveScaledVectorNorm[p, {0.00010264848819015053, 0.00010264848819015053},
NDSolve][diff, InternalMaxAbs[rhs2, ic[2] /@ grid]]
(* In v9 or later: 3808.44 *)
(* In v8 or earlier: 27061.1 *)


Remark

"Are you sure you've reproduced the error estimation of NDSolve?" Yes!:

norms = Trace[
sample[2, Method -> {"MethodOfLines",
"SpatialDiscretization" ->
{"TensorProductGrid",
"MinPoints" -> num}}],
HoldPattern@NDSolveScaledVectorNorm[__][__],
TraceInternal -> True] // Flatten // Union;

norms[[3]] // ReleaseHold
(* In v9 or later: 3808.44 *)
(* In v8 or earlier: 27061.1 *)


For any error estimation larger than 1, NDSolve will consider it as too large and increase number of grid aiming at a better estimation, but as we all know, the peak height of Gibbs phenomenon won't be relieved by denser grid, thus the mxsst warning finally comes up.

So, the problem won't show up when higher order polynomial is used to build the i.c. e.g. ic[3][x], because ic[3]''[x] is still continuous and is not influenced by Gibbs phenomenon:

Plot[ic[3]''[x], {x, xL, xR}, PlotRange -> All]~Show~
ListPlot[{grid, rhsfunc[ic[3], grid]}\[Transpose], PlotStyle -> {Red}]


sample[3]
(* Free of warning. *)

• Great. If there is no objection from your side, I'll use this on mxsst message ref page? Commented Oct 4, 2022 at 7:29
• @user21 Feel free to use it :) . Commented Oct 4, 2022 at 7:43
• This is really an excellent answer! Commented Oct 4, 2022 at 7:44

I still get a warning but I think it's less worrying if I use an alternative form of your initial condition :

b=NDSolve[{D[tes[t,x],t]==D[tes[t,x],x,x]+Exp[-1/(tes[t,x])],
tes[t,0]==1,tes[t,1]==1,
tes[0,x]==1+UnitStep[x-0] UnitStep[1-x]},{tes[t,x]},
{t,0,100},{x,0,1}]


• Haha, yeah, in fact this is the original problem I'm trying to solve, but at last I found that mathematica seems to have a high require for…er…consistency (it seems that the concept "one-sided limit" " One-Sided Derivative" etc. are not available in NDSolve), and this is just the reason why I turned to the problem now you see. Commented Jul 30, 2012 at 14:44