Understanding NDSolve's implementation of the method of lines

Question

In order to better understand the MethodOfLines flag, I tried comparing the results of NDSolve for a driven heat equation with a manual implementation based on pdetoode. Despite my best efforts, the absolute discrepancy at collocation points is at the part per million level which seems rather high. Can anyone suggest a strategy to bring pdetoode into agreement with NDSolve in this case?

Edit: Closer examination suggests that the discrepancy may be related to the temporal grid, which contain 10018 and 10007 points for pdesol and odesol respectively. Is it possible to force the number of temporal grid points into agreement?

T = 10;
f[t_] := Sin[t]^2;
numpoints = 30;
grid = Array[# &, numpoints, {0, 1}];
removeredundant = #[[2 ;; -2]] &;
xdifforder = 4;

ptoofunc = pdetoode[u[x, t], t, grid, xdifforder];

pde = Derivative[0, 1][u][x, t] == Derivative[2, 0][u][x, t];
ic = u[x, 0] == 0;
bc = {u[0, t] == f[t], u[1, t] == 0};
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ode = removeredundant@ptoofunc@pde;
odeic = removeredundant@ptoofunc@ic;
odebc = ptoofunc@bc;

method = {"MethodOfLines", "TemporalVariable" -> t,
   "SpatialDiscretization" -> {"TensorProductGrid", 
     "MaxPoints" -> numpoints, "MinPoints" -> numpoints, 
     "DifferenceOrder" -> xdifforder}};

odesol = NDSolveValue[Join[ode, odeic, odebc], 
                      u /@ grid, {t, 0, T}, MaxStepSize -> 0.001];
pdesol = NDSolveValue[Join[{pde}, {ic}, bc], 
                      u, {x, 0, 1}, {t, 0, T}, Method -> method, 
                      MaxStepSize -> {Automatic, 0.001}];

index = Round[numpoints 0.5];
xcoord = grid[[index]];
Plot[{odesol[[index]][t], pdesol[xcoord, t]}, {t, 0, T}, PlotStyle -> {, Dashed}]
Timing@Plot[odesol[[index]][t] - pdesol[grid[[index]], t], {t, 0, T}, PlotRange -> All]

Answer 1

As you've noticed, the problem is related to the solver in t direction. The ODE/DAE solver is the most hard-to-understand part of NDSolve, so I cannot give a complete discussion for the topic (at least for now), see analysis in the following post for more info:

NDSolve very slow on 2D heat equation

But your specific problem is relatively easy because you just want to avoid the high absolute discrepancy. To achieve this, we simply need to choose exactly the same ODE solver for NDSolve.

We first change the following 2 lines

odeic = removeredundant@ptoofunc@ic;
odebc = ptoofunc@bc;

to

odeic = ptoofunc@ic;
scalefactor = 1;
odebc = diffbc[t, scalefactor]@ptoofunc@bc;

With this modification, the odebc will become ODEs rather than algebraic equations. (See discussion in this post for more info. ) Then we can choose a primary ODE solver (whose setting is fully controllable and easy to control) for the system. (There's no such primary DAE solver in NDSolve AFAIK. )

odemethod = {FixedStep, Method -> ExplicitEuler};
step = 10^-4;

method = {"MethodOfLines", "TemporalVariable" -> t, 
   "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> numpoints, 
     "MinPoints" -> numpoints, "DifferenceOrder" -> xdifforder}, 
          Method -> odemethod};

odesol = NDSolveValue[Join[ode, odeic, odebc], u /@ grid, {t, 0, T}, 
    StartingStepSize -> step, Method -> odemethod]; // AbsoluteTiming
(* {0.479614, Null} *)
pdesol = NDSolveValue[Join[{pde}, {ic}, bc], u, {x, 0, 1}, {t, 0, T}, 
    Method -> method, StartingStepSize -> {Automatic, step}]; // AbsoluteTiming
(* {2.74833, Null} *)

odesol[[index]]["ValuesOnGrid"] - pdesol["ValuesOnGrid"][[index]] // ListLinePlot

As we can see, now the error is around 10^-15, which is likely to be round-off error. If you still feel worried about this error, we can further check with a higher WorkingPrecision:

odesol = 
   NDSolveValue[Join[ode, odeic, odebc], u /@ grid, {t, 0, T}, 
    StartingStepSize -> step, Method -> odemethod, WorkingPrecision -> 16, 
    MaxSteps -> Infinity]; // AbsoluteTiming
(* {14.0822, Null} *)

pdesol = NDSolveValue[Join[{pde}, {ic}, bc], u, {x, 0, 1}, {t, 0, T}, 
    Method -> method, StartingStepSize -> {Automatic, step}, WorkingPrecision -> 16, 
    MaxSteps -> Infinity]; // AbsoluteTiming
(* {18.7217, Null} *)

odesol[[index]]["ValuesOnGrid"] - pdesol["ValuesOnGrid"][[index]] // 
  Abs // Max
(* 0.*10^-16 *)

As we can see, with WorkingPrecision -> 16, the largest absolute error is 0.*10^-16, so I think it's safe enough to say the 2 solutions coincide.