# Fixing the difference order for space derivatives in a PDE when using NDSolve

I am trying to implement the fix to make the difference order uniform for all space derivatives in a PDE, as described in this post: NDSolve uses different difference order for different spatial derivative when solving PDE

It does not work, and I am out of my wits as to why. (Actually, it is the system of 2 coupled PDEs that I am solving.) After the system is set up I included the standard call to NDSolve to show that the problem was set up correctly, and then in the end there is a minimal implementation of the fix. As I wrote it, it does not work, although the the output indicates that time stepping was performed. Can anyone please help?

Clear[tosameorder, fix]
tosameorder[state_NDSolveStateData, order_] :=
state /. a_NDSolveFiniteDifferenceDerivativeFunction :>
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]

RHSCeq = D[H[t, x], x] D[Subscript[C, Bb][t, x], x, x];
RHSHeq = Subscript[C, Bb][t, x] D[H[t, x], x, x];
Nkmax = 6.482250811894023;
Nλmax = 2 Pi/Nkmax;

domainlength = 20 Nλmax; center1 = domainlength/8; w1 = 0.05;
IniShapeH = 1 - 0.02 Exp[-(x - center1)^2/w1^2];
IniShapeC = 1/2;

Δt = 0.1

solution =
NDSolve[{D[Subscript[C, Bb][t, x], t] == RHSCeq,
D[H[t, x], t] == RHSHeq,
Subscript[C, Bb][t, 0] == Subscript[C, Bb][t, domainlength],
H[t, 0] == H[t, domainlength], Subscript[C, Bb][0, x] == IniShapeC,
H[0, x] == IniShapeH,
WhenEvent[Mod[t, Δt] == 0, Print[t]]}, {Subscript[C,
Bb], H}, {x, 0, domainlength}, {t, 0, 1},
AccuracyGoal -> MachinePrecision/2,
PrecisionGoal -> MachinePrecision/2, InterpolationOrder -> All,
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> 2, MinPoints -> 1000, MaxPoints -> 1000,
AccuracyGoal -> MachinePrecision/2,
PrecisionGoal -> MachinePrecision/2},
Method -> {"ImplicitRungeKutta",
"DifferenceOrder" -> 2}}]

Plot[{IniShapeH, First[H[1, x] /. solution]}, {x, 0, domainlength},
PlotRange -> All,
PlotStyle -> {{Thickness[0.0055], Black}, {Thickness[0.0055], Blue}}]

Plot[{IniShapeC, First[Subscript[C, Bb][1, x] /. solution]}, {x, 0,
domainlength}, PlotRange -> All,
PlotStyle -> {{Thickness[0.0055], Black}, {Thickness[0.0055], Blue}}]

endtime = 1; nspacepoints = 1000; xdifforder = 4;

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

nd := NDSolveValue[{D[Subscript[C, Bb][t, x], t] == RHSCeq,
D[H[t, x], t] == RHSHeq,
Subscript[C, Bb][t, 0] == Subscript[C, Bb][t, domainlength],
H[t, 0] == H[t, domainlength], Subscript[C, Bb][0, x] == IniShapeC,
H[0, x] == IniShapeH,
WhenEvent[Mod[t, Δt] == 0, Print[t]]}, {Subscript[C,
Bb], H}, {x, 0, domainlength}, {t, 0, endtime},
Method -> mol[nspacepoints, xdifforder]]

sold = fix[endtime, xdifforder]@nd

• @xzczd Thank you! But how to extract solution from sold ? Plot[{IniShapeH, First[H[x, 0.1] /. sold]}, {x, 0, domainlength}, PlotRange -> All, PlotStyle -> {{Thickness[0.0055], Black}, {Thickness[0.0055], Blue}}] results in error ..
– Mike
Commented Jun 5, 2017 at 1:12

It's actually a matter of evaluation order. Notice that in my answer, I've defined nd with a With[{nd:= ……}, ……] structure. Why? Because the last "step" inside fix is

Unevaluated[ndsolve][[2]] /. NDSolveProcessSolutions@newstate


The Unevaluated[ndsolve][[2]] part is for extracting the dependent variable from ndsolve, which requires the dependent variables to explictly exist. If you simply define nd:=…… without With, nd will still be nd inside Unevaluated. In other words, what happened here can be reproduced with the following example:

lst := {a, b}
Unevaluated[lst][[2]]


Part::partd

In this example, though the OwnValues of lst is {a, b}, when it's not evaluated, it's no more than a Symbol in Part's eyes, so the code doesn't work as expected.

Then why it'll work if I define nd with a With? Because the… Er… scoping mechanism of With is quite similar to rule substitution. When e.g.

With[{lst := {a, b}}, Unevaluated[lst][[2]]]


is executed, lst in Unevaluated[lst][[2]] is simply replaced by {a, b}.

So, to fix your code, we just need to replace nd with the code inside nd e.g. modify the last line to:

sold = Unevaluated@fix[endtime, xdifforder][nd] /. OwnValues@nd

• Thank you ! But how to extract the solution from sold ? Plot[{IniShapeH, First[H[0.1,x] /. sold]}, {x, 0, domainlength}, PlotRange -> All, PlotStyle -> {{Thickness[0.0055], Black}, {Thickness[0.0055], Blue}}] results in error ..
– Mike
Commented Jun 5, 2017 at 1:25
• Notice sold is already a List of InterpolatingFunction so sold[[1]] is the solution for Subscript[C, Bb] and sold[[2]] is the solution for H`. Commented Jun 5, 2017 at 2:47