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@NDSolve`ProcessEquations[{eqn, ic, bc},
    u, {x, 0, 1}, {t, 0, 2}, 
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> 4}}];
funcexpr = state["NumericalFunction"]["FunctionExpression"]

Introduction for NDSolve`ProcessEquations can be found in tutorial/NDSolveStateData and tutorial/NDSolveDAE.

enter image description here

Then check the "DifferenceOrder" of these NDSolve`FiniteDifferenceDerivativeFunction:

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@NDSolve`ProcessEquations[{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?


2 Answers 2


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 NDSolve`FiniteDifferenceDerivativeFunction inside NDSolve`StateData directly:

Clear[tosameorder, fix]
tosameorder[state_NDSolve`StateData, order_] := 
 state /. a_NDSolve`FiniteDifferenceDerivativeFunction :> 
   RuleCondition@NDSolve`FiniteDifferenceDerivative[a@"DerivativeOrder", a@"Coordinates", 
    "DifferenceOrder" -> order, PeriodicInterpolation -> a@"PeriodicInterpolation"]

fix[endtime_, order_] := 
  Module[{state = First[NDSolve`ProcessEquations @@ Unevaluated@ndsolve], newstate}, 
    newstate = tosameorder[state, order]; NDSolve`Iterate[newstate, endtime]; 
   Unevaluated[ndsolve][[2]] /. NDSolve`ProcessSolutions@newstate], HoldAll]


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, θ] == 
   Simplify`PWToUnitStep@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}]]]

enter image description here

  • $\begingroup$ 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} ... $\endgroup$
    – Nobody
    Commented Dec 8, 2021 at 7:16
  • $\begingroup$ @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 NDSolve`StateData[……] 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. $\endgroup$
    – xzczd
    Commented Dec 8, 2021 at 7:41
  • $\begingroup$ @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 :) $\endgroup$
    – Nobody
    Commented Dec 8, 2021 at 9:23
  • $\begingroup$ @zxczd, for the highest derivatives does NDSovle always use implicit representation of finite difference? If not, how could we know the scheme? $\endgroup$
    – lxy
    Commented May 1, 2022 at 4:16
  • $\begingroup$ @jsxs What do you mean by "implicit representation of finite difference"? $\endgroup$
    – xzczd
    Commented May 1, 2022 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.

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}];
    Table[s[[i - 1]] /. t -> tt, {i, 2, n}], {lower}], 
    {tt, 0, endtime, endtime/10}], DataRange -> range, PlotRange -> 1]

enter image description here

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.

enter image description here

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}]

enter image description here

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%,

enter image description here

and yields for t == endtime,

enter image description here

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.

  • $\begingroup$ Interesting! You may also give an answer to that question, I think your difference scheme is exactly what ojlm needs. $\endgroup$
    – xzczd
    Commented Dec 11, 2016 at 6:15
  • 1
    $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Commented Dec 11, 2016 at 7:14
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Commented Dec 15, 2016 at 13:22
  • $\begingroup$ 教学相长 :D $\endgroup$
    – xzczd
    Commented Dec 15, 2016 at 13:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.