2
$\begingroup$

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_NDSolve`StateData, order_] := 
state /. a_NDSolve`FiniteDifferenceDerivativeFunction :> 
NDSolve`FiniteDifferenceDerivative[a@"DerivativeOrder", 
a@"Coordinates", "DifferenceOrder" -> order, 
PeriodicInterpolation -> a@"PeriodicInterpolation"]

fix[endtime_, order_] := 
Function[{ndsolve}, 
Module[{state = 
First[NDSolve`ProcessEquations @@ Unevaluated@ndsolve], 
newstate}, newstate = tosameorder[state, order]; 
NDSolve`Iterate[newstate, endtime];
Unevaluated[ndsolve][[2]] /. NDSolve`ProcessSolutions@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", 
  "Coefficients" -> "ImplicitRungeKuttaRadauIIACoefficients", 
  "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
$\endgroup$
  • $\begingroup$ @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 .. $\endgroup$ – Mike Jun 5 '17 at 1:12
1
$\begingroup$

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]] /. NDSolve`ProcessSolutions@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
$\endgroup$
  • $\begingroup$ 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 .. $\endgroup$ – Mike Jun 5 '17 at 1:25
  • $\begingroup$ 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. $\endgroup$ – xzczd Jun 5 '17 at 2:47

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.