# NDSolve very slow on 2D heat equation

I am trying to solve the 2D heat equation

$$\begin{cases} u_{t}-u_{x x}-u_{y y}=f \\ u(0, x, y)=\sin (2 \pi x) \sin (2 \pi y) \\ u(t, 0, y)=0 \\ u(t, x, 1)=0 \\ u_{x}(t, 1, y)=2 \pi e^{-t} \sin (2 \pi y) \\ u_{y}(t, x, 0)=2 \pi e^{-t} \sin (2 \pi x) \end{cases} \quad x \in(0,1), y \in(0,1), t \in(0, T)$$ where $$f=8 \pi^{2} e^{-t} \sin (2 \pi x) \sin (2 \pi y)-e^{-t} \sin (2 \pi x) \sin (2 \pi y)$$ using NDSolve with the following code.

NDSolveValue[{Derivative[1, 0, 0][u][t, x, y] - Laplacian[u[t, x, y], {x, y}] ==
8 π^2 E^-t Sin[2 π x]*Sin[2 π y] - E^-t Sin[2 π x]*Sin[2 π y],
u[0, x, y] == Sin[2 π x]*Sin[2 π y],
u[t, 0, y] == 0, u[t, x, 1] == 0,
Derivative[0, 1, 0][u][t, 1, y] == 2 π E^-t Sin[2 π y],
Derivative[0, 0, 1][u][t, x, 0] == 2 π E^-t Sin[2 π x]},
u[t, x, y], {x, 0, 1}, {y, 0, 1}, {t, 0, 1}]


The exact solution is $$u(t, x, y)=e^{-t} \sin (2 \pi x) \sin (2 \pi y)$$.

I tried using low PrecisionGoal (3), AccuracyGoal (3), and MaxSteps (10^4), but the error is unacceptable, especially on the boundaries. Solutions from NDSolve are in the first row and the exact solution is in the second.

What would be the best method and options for solving such an equation?

• What is your question? What error were acceptable compared to which solution? Jan 17 at 9:10
• How long did it take to complete for you? few minutes? much more? can you be more exact. Its been running for me for 5 minutes but still not completed. Also it will be good to say which version you are using. Jan 17 at 9:21
• @Nasser, When I set 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. Jan 17 at 9:29
• You might want to report this to WRI ([email protected]). NDSolve should not hang. Jan 17 at 9:59
• @Nasser, the model does not hang, at least for me, you can easily see this if you add an evaluation monitor. Jan 17 at 10:15

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[FrontEndEvaluationNotebook[],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[], diffbc[t, sf2]@bcx[]},
ptoofunc2@{diffbc[t, sf1]@bcy[], diffbc[t, sf2]@bcy[]}}];
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} = NDSolveProcessEquations[seq2];
func2 = state2["NumericalFunction"];

{state} = NDSolveProcessEquations[seq];
func = state["NumericalFunction"];

iclst = N@Flatten@Table[ic[], {x, grid}, {y, grid}];

ListLinePlot[func[tend, iclst] - func2[tend, iclst], PlotRange -> All] Remark

1. The showStatus is from this post.

2. 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 == iclst}, u, {t, 0, tend},
jianshi[t]]; // AbsoluteTiming


So, there seems to be something mysterious in NDSolveStateData, but deeper analysis is beyond my reach.

• Thank you for the elaborate answer! Jan 18 at 2:35
• @xzczd Nice and interesting answer! I would like to know which Method instead of FiniteElement is used in your sol but sol["MethodInfomation"] gives no information. Jan 18 at 13:24
• @UlrichNeumann In this case TensorProductGrid is used for spatial discretization. (This can be checked with state["OptionValues"]. See also mathematica.stackexchange.com/q/140773/1871 ) As to method for time integration, when SolveDelayed->True is set, DAE solver (to be more specific, IDA method) has been used. Jan 18 at 13:43
• @xzczd Thank you for your helpful explanation! Jan 18 at 13:45

Try FiniteElement method with NeumannValue instead of derivative-bc:

U = NDSolveValue[{Derivative[1, 0, 0][u][t, x, y] -
Laplacian[u[t, x, y], {x, y}] ==
8 \[Pi]^2 E^-t Sin[2 \[Pi] x]*Sin[2 \[Pi] y] -
E^-t Sin[2 \[Pi] x]*Sin[2 \[Pi] y] +
NeumannValue[2 \[Pi] E^-t Sin[2 \[Pi] y], x == 1] -
NeumannValue[2 \[Pi] E^-t Sin[2 \[Pi] x], y == 0],
u[0, x, y] == Sin[2 \[Pi] x]*Sin[2 \[Pi] y], u[t, 0, y] == 0,
u[t, x, 1] == 0 }, u , {x, 0, 1}, {y, 0, 1}, {t, 0, 1}]

Manipulate[
Plot3D[{U[t, x, y], Exp[-t] Sin[2 Pi x] Sin[2 Pi y]}, {x, 0, 1},{y,0, 1}], {{t, 1/3}, 0, 1, Appearance -> "Labeled"}] Solution fits good. Hope it helps!

• Thank you, this is nice. It is weird that NDSolve with default options fails to solve this common equation in a few minutes. Jan 17 at 9:56
• @Charmbracelet The Method -> FiniteElement is not necessary. See this Jan 17 at 9:59
• @user21 Your link only shows example with DirichletConditions at the boundary. OP question concernes first order derivative boundaries Jan 17 at 10:04
• And removing Method -> "FiniteElement" will make the code more efficient, because currently you're using pure FEM to solve the problem (including in t direction). See this post for more info: mathematica.stackexchange.com/a/222373/1871 Jan 18 at 3:45
• I went ahead and removed the Method->"FiniteElement" options. There was a message embezzlement, that we should follow. People will see this, not understand it and make use of it when they should not. Time dependent PDEs should not be solved as a purely spacial problem - that's what you do when you specify Method->"FiniteElement" because this leads to a convection dominant PDE problem. Jan 19 at 15:04