The simplest solution I find is, add SolveDelayed -> True
to NDSolve
:
tend = 1; {xL, xR} = {yL, yR} = {0, 1};
sys = {eq, ic, bcx, bcy} =
With[{u = u[t, x, y]},
{D[u, t] - Laplacian[u,{x, y}] ==
8 π^2 E^-t Sin[2 π x] Sin[2 π y] - E^-t Sin[2 π x] Sin[2 π y],
u == Sin[2 π x] Sin[2 π y] /. t -> 0,
{u == 0 /. x -> xL,
D[u, x] == 2 π E^-t Sin[2 π y] /. x -> xR},
{u == 0 /. y -> yL,
D[u, y] == 2 π E^-t Sin[2 π x] /. y -> yR}}];
sol = NDSolveValue[sys, u, {t, 0, tend}, {x, xL, xR}, {y, yL, yR},
SolveDelayed -> True]; // AbsoluteTiming
(* {5.83636, Null} *)
The option SolveDelayed
is red, but don't worry. Alternatively, you can set Method -> {"EquationSimplification" -> "Residual"}
or Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> False}
.
If you just want to know a solution, you can stop here.
Sadly, I can't explain why the method above works. I thought the problem is related to following posts:
Boundary condition with spatial derivative is ignored by NDSolve
NDSolve uses different difference order for different spatial derivative when solving PDE
Nevertheless, even if I discretize the system in the same way (in principle) as that of NDSolve
using pdetoode
, NDSolve
behaves differently. The following is the code I use for test:
showStatus[status_]:=LinkWrite[$ParentLink,
SetNotebookStatusLine[FrontEnd`EvaluationNotebook[],ToString[status]]];
clearStatus[]:=showStatus[""];
clearStatus[]
jianshi[t_]:=EvaluationMonitor:>showStatus["t = "<>ToString[CForm[t]]]
mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}
points = 50; difforder = 4;
seq := Sequence[sys, u, {t, 0, tend}, {x, xL, xR}, {y, yL, yR},
Method -> mol[points, difforder], jianshi[t]];
(* Warning: the following stucks around t = 0.0035 *)
(* sol = NDSolveValue[seq]; // AbsoluteTiming *)
(* $Aborted *)
grid = Array[# &, points, {xL, xR}];
difforder2 = difforder + 1;
ptoofunc = pdetoode[u[t, x, y], t, {grid, grid}, difforder];
ptoofunc2 = pdetoode[u[t, x, y], t, {grid, grid}, difforder2];
del = #[[2 ;; -2]] &;
ode = del /@ del@ptoofunc@eq;
odeic = ptoofunc@ic;
odebc = With[{sf1 = 1,
sf2 = 0}, {del /@ ptoofunc2@{diffbc[t, sf1]@bcx[[1]], diffbc[t, sf2]@bcx[[2]]},
ptoofunc2@{diffbc[t, sf1]@bcy[[1]], diffbc[t, sf2]@bcy[[2]]}}];
vars = Outer[u, grid, grid];
seq2 :=
Sequence[{ode, odeic, odebc}, vars, {t, 0, 1}, jianshi[t],
Method -> {"EquationSimplification" -> "Solve"}];
sollst = NDSolveValue[seq2]; // AbsoluteTiming
(* {9.63662, Null} *)
sol2 = rebuild[sollst, {grid, grid}];
{state2} = NDSolve`ProcessEquations[seq2];
func2 = state2["NumericalFunction"];
{state} = NDSolve`ProcessEquations[seq];
func = state["NumericalFunction"];
iclst = N@Flatten@Table[ic[[2]], {x, grid}, {y, grid}];
ListLinePlot[func[tend, iclst] - func2[tend, iclst], PlotRange -> All]

Remark
The showStatus
is from this post.
Notice I've defined odebc
in an unusual way to mimic the behavior of PDE solver of NDSolve
. To understand why it's defined in this
way, please refer to the linked posts above.
As we can see, func[tend, iclst] - func2[tend, iclst]
(func
and func2
are NumericalFunction
generated by NDSolve
) is of 10^-12
order, which is likely to be caused by round-off error. Nevertheless, NDSolve
solves the system generated by pdetoode
in 10 seconds, but stucks when solving the PDE directly. At this point, the most reasonable explanation for the strange behavior seems to be, NDSolve
isn't doing a good job in choosing the ODE solver when dealing with the PDE, but:
sollst =
NDSolveValue[{ode, odeic, odebc}, vars, {t, 0, 1}, jianshi[t],
Method -> {"EquationSimplification" -> "Solve",
TimeIntegration -> StiffnessSwitching}]; // AbsoluteTiming
(* {15.6412, Null} *)
(* Warning: the following is very slow *)
sol = NDSolveValue[sys, u, {t, 0, tend}, {x, xL, xR}, {y, yL, yR},
Method -> Append[mol[points, difforder], Method -> StiffnessSwitching],
jianshi[t]]; // AbsoluteTiming
This isn't the end. Using the NumericalFunction
of state2
to build an equivalent (in principle) ODE system leads to trouble, too!:
(* The following stucks around t == 0.003: *)
NDSolveValue[{u'[t] == func2[t, u@t], u[0] == iclst}, u, {t, 0, tend},
jianshi[t]]; // AbsoluteTiming
So, there seems to be something mysterious in NDSolve`StateData
, but deeper analysis is beyond my reach.
PrecisionGoal -> 3, AccuracyGoal -> 3, MaxSteps -> 10^4
, it takes about 10 minutes to complete but it is not close to the exact solution as in the picture. $\endgroup$