# NDSolve uses different difference order for different spatial derivative when solving PDE

I found something this tutorial for method of line doesn't tell us.

Consider the following toy example:

eqn = With[{u = u[x, t]},
D[u, t] == D[u, x] + D[u, {x, 2}] + D[u, {x, 3}] - D[u, {x, 4}]];

ic = u[x, 0] == 0;
bc = {u[0, t] == 0, u[1, t] == 0, D[u[x, t], x] == 0 /. {{x -> 0}, {x -> 1}}};

NDSolve[{eqn, ic, bc},
u, {x, 0, 1}, {t, 0, 2},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> 4}}]


Guess what difference order is chosen when those spatial derivatives (in this case $\frac{\partial u}{\partial x}$, $\frac{\partial ^2u}{\partial x^2}$, $\frac{\partial ^3u}{\partial x^3}$, $\frac{\partial ^4u}{\partial x^4}$) are discretized?

"What a needless question! The order is 4, as we set with "DifferenceOrder" -> 4! " About an hour ago, I thought so, too. But it's not true. Let's check the difference formula generated by NDSolve:

state = First@NDSolveProcessEquations[{eqn, ic, bc},
u, {x, 0, 1}, {t, 0, 2},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> 4}}];
funcexpr = state["NumericalFunction"]["FunctionExpression"]


Introduction for NDSolveProcessEquations can be found in tutorial/NDSolveStateData and tutorial/NDSolveDAE.

Then check the "DifferenceOrder" of these NDSolveFiniteDifferenceDerivativeFunction:

Head[#]@"DifferenceOrder" & /@ funcexpr[[2, 1]]
(* {{7}, {6}, {5}, {4}} *)


So, for a PDE whose maximum spatial differential order is omax, when "DifferenceOrder" -> n is set for "TensorProductGrid", the actual difference order for m-order spatial derivative is omax + n - m.

In certain cases, this design seems to cause trouble, here's an example.

To make this post a question, I'd like to ask:

1. Why NDSolve chooses this design?

2. If the 1st question is too hard, is there a easy way (e.g. a hidden option) to make NDSolve use the same difference order for every spatial derivative?

Note: fix is broken since v11.3, a new question has been started aiming at upgrade it.

Here's my approach for fixing the difference order. The key idea is modifying the NDSolveFiniteDifferenceDerivativeFunction inside NDSolveStateData directly:

Clear[tosameorder, fix]
tosameorder[state_NDSolveStateData, order_] :=
state /. a_NDSolveFiniteDifferenceDerivativeFunction :>
RuleCondition@NDSolveFiniteDifferenceDerivative[a@"DerivativeOrder", a@"Coordinates",
"DifferenceOrder" -> order, PeriodicInterpolation -> a@"PeriodicInterpolation"]

fix[endtime_, order_] :=
Function[{ndsolve},
Module[{state = First[NDSolveProcessEquations @@ Unevaluated@ndsolve], newstate},
newstate = tosameorder[state, order]; NDSolveIterate[newstate, endtime];
Unevaluated[ndsolve][[2]] /. NDSolveProcessSolutions@newstate], HoldAll]


Example:

bound = 0.25510204081632654;
upper = 99/100; lower = 1 - upper;
range = {L, R} = {-Pi/2, Pi/2};
endtime = 100;
xdifforder = 4;
eqn = With[{h = h[t, θ], ϵ = 5/10},
0 == -D[h, t] + D[h^3 (1 - h)^3 ϵ D[h, θ], θ]];
ic = h[0, θ] ==
SimplifyPWToUnitStep@Piecewise[{{upper, -bound < θ < bound}}, lower];
bc = {h[t, L] == lower, h[t, R] == lower};

mol[n_Integer, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}

With[{nd :=
NDSolveValue[{eqn, ic, bc}, h, {t, 0, endtime}, {θ, L, R},
Method -> mol[200, xdifforder], MaxSteps -> Infinity]},
With[{sol = nd, sold = fix[endtime, xdifforder]@nd},
Animate[Plot[{sol[t, th], sold[t, th]}, {th, L, R}, PlotRange -> {0, 1},
PlotLegends -> {"Before fix", "After fix"}], {t, 0, endtime}]]]


Complete control over the spatial decomposition of the PDE given in the answer by xzczd can be achieved by decomposing the PDE into a large set of ODEs, as described in the Introduction to the Numerical Method of Lines, provided in the Mathematica documentation. The following straightforward approach uses a uniform grid and second-order differencing.

Clear[u];
n = 200; d = (R - L)/n;
vars = Table[u[i, t], {i, 2, n}]; u[1, t] = lower; u[n + 1, t] = lower;
eq = Table[dup = (u[i + 1, t] - u[i, t])/d; dum = (u[i, t] - u[i - 1, t])/d;
up = (u[i + 1, t] + u[i, t])/2; um = (u[i, t] + u[i - 1, t])/2;
D[u[i, t], t] == (up^3 (1 - up)^3 dup - um^3 (1 - um)^3 dum) ϵ/d, {i, 2, n}];
init = Table[u[i, 0] == Piecewise[{{upper, -bound < L + (i - 1) d < bound}}, lower],
{i, 2, n}];
s = NDSolveValue[{eq, init}, vars, {t, 0, endtime}];
ListLinePlot[Evaluate@Table[Join[{lower},
Table[s[[i - 1]] /. t -> tt, {i, 2, n}], {lower}],
{tt, 0, endtime, endtime/10}], DataRange -> range, PlotRange -> 1]


A test of the accuracy of this result can be obtained by noting that the integral of D[h, t] (using the nomenclature in the answer by xzczd) over range is given by

h^3 (1 - h)^3 ϵ D[h, θ]


evaluated at R minus the same quantity evaluated at L. Moreover, numerical evaluation of this quantity at the two endpoints shows that it is very small. In other words, the integral of h over range should be essentially constant in time. The solution obtained here indeed is constant when integrated over range, as can be shown by evaluating

Table[Total@N@Table[s[[i - 1]] /. t -> tt, {i, 2, n}] d, {tt, 0, endtime, endtime/20}]
(* {0.539254, 0.539254, ..., 0.539254, 0.539254} *)


Consider, now, the "before fix" and "after fix" solutions obtained by xzczd and plotted here for t == endtime.

The "after fix" solution is similar but not identical to the solution t == endtime curve shown in the first plot of this answer. Moreover, the conserved quantity just described also varies in time.

ListPlot[Table[Quiet@NIntegrate[sold[t, th], {th, L, R},
Method -> {Automatic, "SymbolicProcessing" -> False}],
{t, 0, endtime, endtime/20}], DataRange -> {0, endtime}]


All this is not to suggest that xzczd's elegant answer (+1) is incorrect. In fact, merely increasing the number of grid points to 5000 reduces temporal variation of the conserved quantity in the "after fix" solution to within 0.5%,

and yields for t == endtime,

and the "after fix" curve is identical to the eye to the t == endtime curve in the first plot of this answer. Note that increasing the number of grid points does nothing to improve the accuracy of the "before fix" solution.

• Interesting! You may also give an answer to that question, I think your difference scheme is exactly what ojlm needs. – xzczd Dec 11 '16 at 6:15
• @xzczd As I look more at the ojlm's approach, I see that he analytically differentiates the right side of the PDE and then replaces the first and second spatial derivatives by their finite differences. My approach is to form the second-order finite difference without first analytically expanding the derivative, which I believe to be more accurate. I shall look in greater detail at this distinction tomorrow morning. – bbgodfrey Dec 11 '16 at 7:14
• @xzczd You ask interesting questions, and I have learned quite a bit from this one and others of yours that I have examined this past week. Thanks also for the bounty. – bbgodfrey Dec 15 '16 at 13:22
• 教学相长 :D – xzczd Dec 15 '16 at 13:35