# 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}} *)


The order is not 4! Similarly, we can verify that it's the same case for PDE system:

eqn =
With[{u = u[x, t], v = v[x, t]},
{D[u, t] == D[u, x] + D[u, {x, 2}] + D[u, {x, 3}] - D[u, {x, 4}],
D[v, t] == D[v, x] + D[v, {x, 2}]}];
ic = {u[x, 0] == 0, v[x, 0] == 0};
bc = {{u[0, t] == 0, u[1, t] == 0,
D[u[x, t], x] == 0 /. {{x -> 0}, {x -> 1}}},
{v[0, t] == 0, v[1, t] == 0}};

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

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


So, for a PDE or PDE system 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][] /. 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}]]] • two questions to consult: 1. is there a way to check that with your fix function, NDSolve indeed uses the same difference order for all the spatial derivatives? 2. Whether or not the fix function can be applied when solving 2 coupled equations with NDSolveValue? Btw, I found when solve 2 eqns with different highest spatial-derivative orders, your conclusion in OP would be modified as For a system of PDEs whose maximum spatial differential order is omax=max{eqn_1, eqn_2, ... eqn_i} ... Dec 8, 2021 at 7:16
• @Nobody 0. Do notice fix has been broken since v11.3. 1. Well, I think it's clear if one understands the source code of fix. Anyway, if you want to check, just add a Sow to my fix and use Reap to extract the modified NDSolveStateData[……] and do what I've done in the question. 2. Yes, it works for PDE system. 3. Thx for pointing out. The statement in the question is revised. Dec 8, 2021 at 7:41
• @zxczd, thanks, I am still a big fan of v.9, based on my personal experience it is often better than v11.3 at least in the aspect of numerical analysis :) Dec 8, 2021 at 9:23
• @zxczd, for the highest derivatives does NDSovle always use implicit representation of finite difference? If not, how could we know the scheme?
– jsxs
May 1 at 4:16
• @jsxs What do you mean by "implicit representation of finite difference"? May 1 at 5:53

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. Dec 11, 2016 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. Dec 11, 2016 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. Dec 15, 2016 at 13:22
• 教学相长 :D Dec 15, 2016 at 13:35