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

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?


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 :> 
   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


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 Dec 11 '16 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 Dec 11 '16 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 Dec 15 '16 at 13:22
  • $\begingroup$ 教学相长 :D $\endgroup$ – xzczd Dec 15 '16 at 13:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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