# Boundary condition for heat equation in polar coordinates deduced with L'Hôpital's rule fails for method of lines, but works well for FDM

## Update

I manage to find a way that resolves the problem. I'd like not to make this solution public for the moment so other answerers will have more chance to get the bounty. Here's a hint: it's a documented but seldom used Method combination.

I've hesitated for a while about whether this should be posted to scicomp.SE or not, but still decide to discuss it here because I believe this will make the discussion more convenient.

All the code in this question has been tested in v14.1, v12.3 and v8.0.4 (with minimal modification to make it v8 compatible).

Consider the following initial-boundary value problem for 1+1D axisymmetric heat equation in polar coordinates: $$\newcommand{\p}{\partial}\newcommand{\f}{\frac}$$ $$\left\{ \begin{array}\\ \f{\p u}{\p t} = \f{1}{r} \f{\p u}{\p r}+ \f{\p^2 u}{\p r^2}\\ \left.u\right|_{t=0}=1-r^2\\ \left.u\right|_{r=1}=0\\ \end{array} \right.$$

It's well known that the problem has a removable singularity at $$r=0$$. If one wants to solve the problem numerically in $$r\in[0,1]$$, the singularity should be handled properly. There exist many possible methods for handling the singularity, of course. For example, use FiniteElement method of NDSolve:

R = 1;
tend = 1;

sys = With[{u = u[t, r]}, {D[u, t] == Laplacian[u, {r, θ}, "Polar"],
u == 1 - r^2 /. t -> 0, u == 0 /. r -> R}];

solfem = NDSolveValue[sys, u, {t, 0, tend}, {r, 0, R},
Method -> {"MethodOfLines", "SpatialDiscretization" -> "FiniteElement"}];


But this question will focus on one method only i.e. the method based on L'Hôpital's rule, which has been mentioned in various sources e.g. this and this post, and the paper A computational fluid dynamic technique valid at the centerline for non-axisymmetric problems in cylindrical coordinates by Griffin. The idea is simple: given that the $$\f{1}{r} \f{\p u}{\p r}+ \f{\p^2 u}{\p r^2}$$ term is finite at $$r=0$$, we can transform it based on L'Hôpital's rule:

$$\lim_{r\to 0}\left(\f{1}{r} \f{\p u}{\p r}+ \f{\p^2 u}{\p r^2}\right)=\lim_{r\to 0}\left(\f{\left(\f{\p\f{\p u}{\p r}}{\p r}\right) }{\left(\f{\p r}{\p r}\right)} + \f{\p^2 u}{\p r^2}\right)=\left.2\f{\p^2 u}{\p r^2}\right|_{r=0}$$

For this specific problem, $$\left.u\right|_{r=1}=0$$ is equivalent to $$\left.\f{\p u}{\p t}\right|_{r=1}=0$$. So we have:

$$\left\{ \begin{array}\\ \f{\p u}{\p t} = \f{1}{r} \f{\p u}{\p r}+ \f{\p^2 u}{\p r^2}\\ \left.u\right|_{t=0}=1-r^2\\ \left.\f{\p u}{\p t}\right|_{r=1}=0\\ \left.\left(\f{\p u}{\p t} -2\f{\p^2 u}{\p r^2}\right)\right|_{r=0}=0 \end{array} \right.$$

OK, here comes the trouble. At this point, method of lines (the method internally used by NDSolve, which discretizes the problem in $$r$$ direction with finite difference formula and solve the obtained ODE system with any ODE solver) seems to be a reasonable choice for solving the problem, but it's not true. Just try the following implementation (though my pdetoode can be used to facilitate coding, I've avoided it to make the code more transparent):

n = 100;
Δr = R/(n - 1);

rhs[u_List] :=
Table[Which[i == 1,
2 {2, -5, 4, -1} . {u[[1]], u[[2]], u[[3]], u[[4]]}/Δr^2,

i == n,
0.,

True,
With[{r = (i - 1) Δr}, 1/r (u[[i + 1]] - u[[i - 1]])/(2 Δr)
+ (u[[i + 1]] - 2 u[[i]] + u[[i - 1]])/Δr^2]], {i, n}]

iclst = Table[1. - r^2, {r, 0, R, Δr}];

sol = NDSolveValue[{u'[t] == rhs@u[t], u[0] == iclst}, u, {t, 0, tend/500}];

Plot[sol[t][[1]], {t, 0, tend/500}, PlotRange -> All]


Remark

At $$r=0$$, the one-sided difference formula $$f''\left(x_i\right)\simeq \frac{2 f\left(x_i\right)-5 f\left(x_{i+1}\right)+4 f\left(x_{i+2}\right)-f\left(x_{i+3}\right)}{h^2}$$ has been used.

Clearly, the solution is unstable.

"Well, perhaps the method based on L'Hôpital's rule is simply wrong?" It doesn't seem to be, because if we further discretize in $$t$$ direction with finite difference formula and solve the resulting linear algebraic equation system, the desired solution is obtained!:

m = 100;
Δt = tend/(m - 1);

eqlst = Flatten@Table[With[{r = (i - 1) Δr},
With[{lhs =
Which[j == 1 || i == 1 || i == n,
u[j, i],
j == m,
{1, -4, 3} . {u[j - 2, i], u[j - 1, i], u[j, i]}/(2 Δt),
True,
(u[j + 1, i] - u[j - 1, i])/(2 Δt)],

rhs =
Which[j == 1, 1 - r^2,
i == 1, (2 {2, -5, 4, -1}.{ u[j, 1], u[j, 2], u[j, 3], u[j, 4]})/Δr^2,
i == n,
0,

True,
1/r (-u[j, -1 + i] + u[j, 1 + i])/(2  Δr) +
(u[j, -1 + i] - 2 u[j, i] + u[j, 1 + i])/Δr^2]},

lhs == rhs]], {j, m}, {i, n}];

vars = Array[u, {m, n}] // Flatten;
sollst = LinearSolve[#2, -N@#] & @@
CoefficientArrays[eqlst, vars]; // AbsoluteTiming

solfdm = ListInterpolation[ArrayReshape[sollst, {m, n}], {{0, tend}, {0, R}}];

Manipulate[
Plot[{solfem[t, r], solfdm[t, r]}, {r, 0, R}, PlotStyle -> {Automatic, Dashed},
PlotRange -> {0, 1.1}], {t, 0, tend}, ControlPlacement -> Top]


Remark

At $$t=1$$, 2nd order one-sided difference formula $$f' (x_n)\simeq \frac{f (x_{n}-2h)-4 f (x_{n}-h)+3 f (x_n)}{2 h}$$ has been used. This is not necessary, the simple $$f'\left(x_i\right)\simeq \frac{f\left(x_{i}\right)-f\left(x_{i-1}\right)}{h}$$ with a dense enough grid will result in the desired solution, too.

Why does the pure finite difference method (FDM) succeed, but the method of lines fail?

Just found something surprising. If we solve the problem with NDSolve only i.e. make use of the automatic discretization of NDSolve, the problem can be solved correctly!:

newsys = With[{u = u[t, r]}, {D[u, t] ==
PiecewiseExpand@If[r == 0, 2 D[u, r, r], Laplacian[u, {r, θ}, "Polar"]],
u == 1 - r^2 /. t -> 0, u == 0 /. r -> R}];

argu = {newsys, u, {t, 0, tend}, {r, 0, R},
Method -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 0},
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 100,
"MinPoints" -> 100, "DifferenceOrder" -> 2}},
MaxStepSize -> {tend/50, Automatic}};

tst = NDSolveValue@@argu;

Plot[{solfem[tend, r], tst[tend, r]}, {r, 0, R}, PlotStyle -> {Automatic, Dashed}]


Remark

1. The bcart warning will show up, but it's merely a warning in this case. It shows up because there's no explicit b.c. at $$r=0$$ in the code. (We have to hide the b.c. at $$r=0$$ into the PDE via Piecewise, otherwise bdord warning will show up and NDSolve will fail in parsing the problem. )

2. The Method option isn't actually necessary, I set them just to make the underlying setting closer to rhs above.

But the difference scheme used by NDSolve is almost equivalent to the rhs above:

state = First[NDSolveProcessEquations @@ argu];

ndsolverhs = state["NumericalFunction"];

grid = Array[# &, n, {0, R}];

ndsolverhs[12345, iclst] - rhs@iclst // Abs // Max
(* 3.18368*10^-12 *)


Remark

What's more surprising is, NDSolve will fail on the ODE system re-built with the ExperimentalNumericalFunction[…] generated by NDSolve!:

tst2 = NDSolveValue[{u'[t] == ndsolverhs[23456, u[t]],
u[0] == iclst}, u, {t, 0, tend/500}];

Plot[tst2[t][[1]], {t, 0, tend/500}]


This might be related to the observation here.

If it's an ODE solver issue, which method for ODE solving should we choose to solve the problem correctly? If not, what magic has been used by NDSolve in definition of tst for stablizing the solution?

Just in case, here's the implementation of method of lines via pdetoode:

R = 1; tend = 1; n = 100; difforder = 2;

With[{u = u[t, r]},
{D[u, t] == Laplacian[u, {r, θ}, "Polar"],
u == 1 - r^2 /. t -> 0,
u == 0 /. r -> R,
D[u, t] == 2 D[u, r, r] /. r -> 0}];

grid = Array[# &, n, {0, R}];
(* Definition of pdetoode isn't included in this post,
ptoofunc = pdetoode[u[t, r], t, grid, difforder];

ode = #[[2 ;; -2]] &@ptoofunc@eq; // Quiet
odeic = ptoofunc@ic;

solfunclst = NDSolveValue[{ode, odeic, odebc}, u /@ grid, {t, 0, tend/500}];
solfunc = rebuild[solfunclst, grid];
Plot[solfunc[t, 0], {t, 0, tend/500}, PlotRange -> All]


As we can see, the solution is unstable, just similar to sol. (Yeah it's not exactly the same, and the reason isn't fully understood, either. I suspect it's essentially the same question as asked above. )

• This is a question about the stability of the numerical method. Since the system of ordinary differential equations in this case is linear, then answering this question does not present any difficulty. It is possible to build a code to test the convergence. Don't pay attention to the NDSolve "magic". It seems that FEM and MOL ignoring L'Hôpital's rule and simply used D[u[t,r],r]=0 at r=0. Commented Aug 6 at 6:39
• @AlexTrounev "It seems that FEM and MOL ignoring L'Hôpital's rule and simply used D[u[t, r], r] == 0 at r == 0." I haven't looked into how the singularity is handled by FiniteElement, but for TensorProduct sub-method, by checking state["NumericalFunction"]["FunctionExpresion"] we can see the b.c. deduced from L'Hôpital's rule is indeed used, so I don't think they've just used D[u[t, r], r] == 0 under the hood. (At least I can't find any evidence for this guess. ) Commented Aug 6 at 7:53
• Let check Table[Derivative[1, 0][#][t, 0] - 2 Derivative[0, 2][#][t, 0] & /@ {solfem, tst}, {t, 0, 1, .01}] and Table[Derivative[0, 1][#][t, 0] & /@ {solfem, tst}, {t, 0, 1, .01}]. Max abs values in first table is about 5.4 10^-3, while second one is about 1.5 10^-5. :) Commented Aug 6 at 9:44
• @alex I'd argue the comparison is a bit unfair, because the interpolation error must be considered. (Related: mathematica.stackexchange.com/a/145853/1871 ) To minimize the influence of interpolation error, I believe we should compare at integration points in t direction i.e. tlst = tst["Coordinates"][[1]]; Table[(Derivative[1, 0][#1][t, 0] - 2 Derivative[0, 2][#1][t, 0] &)@tst, {t, tlst}] // {Histogram@#, ListLinePlot[#, PlotRange -> All]} & We can see the error in most position is actually small. As to the few steep peaks, I suspect it's a sign of internal method changing. Commented Aug 6 at 11:11
• Do you really believe that NDSolve can accept such a boundary condition derived from L'Hôpital's rule? ;) Commented Aug 6 at 12:01

Now that the bounty has ended, let me post my answer. Still, this isn't a complete answer because I haven't figured out why the default ODE/DAE solver of NDSolve behaves like this, but I believe it's safe to say "it's an ODE/DAE solver issue" now, because I manage to find a Method combination that handles the problem properly.

As mentioned in the question:

…the simple $$f'\left(x_i\right)\simeq \frac{f\left(x_{i}\right)-f\left(x_{i-1}\right)}{h}$$ with a dense enough grid will result in the desired solution, too.

$$f'\left(x_i\right)\simeq \frac{f\left(x_{i}\right)-f\left(x_{i-1}\right)}{h}$$ in $$t$$ direction is backward Euler method, which is a type of implicit Runge–Kutta method (see e.g. this post for more info), this gives me a hint. After some trial and error, I found that combination of

1. "FixedStep" controller method

2. "ImplicitRungeKutta" time integration method

3. A proper "Coefficients" setting for "ImplicitRungeKutta"

resolves the problem!:

sol = NDSolveValue[{u'[t] == rhs@u[t], u[0] == iclst}, u, {t, 0, tend},
Method -> {"FixedStep",
Method -> {"ImplicitRungeKutta",
StartingStepSize -> tend/25];

ListPointPlot3D[sol, PlotRange -> All, DataRange -> {{0, R}, {0, tend}},
AxesLabel -> {"r", "t"}]


Definitions of rhs, etc. are the same as in the question so omitted here.

Remark:

1. According to my test, "ImplicitRungeKuttaLobattoIIICCoefficients", "ImplicitRungeKuttaRadauIACoefficients" and "ImplicitRungeKuttaRadauIIACoefficients" all work well for the problem.

2. "ImplicitRungeKuttaGaussCoefficients" can also be used, but only for a not-that-dense and not-that-coarse grid. For example Method -> {"FixedStep", Method -> {"ImplicitRungeKutta", "Coefficients" -> "ImplicitRungeKuttaGaussCoefficients"}}, StartingStepSize -> tend/1000 causes ndcf warning and fail at about t = 0.193, StartingStepSize -> tend/25 results in a solution whose error is more obvious compared with the above three (just check the solution at t == tend), StartingStepSize -> tend/100 gives a solution that looks good.

3. For a really dense grid in t direction e.g. StartingStepSize -> tend/10^4, one may need to manually set the "DifferenceOrder" sub-option of "ImplicitRungeKutta" to e.g. 2.

• (+1) This is perfect solution. Also, it is wander why LSODE solver used by NDSolve in PDE form is not working in ODEs form? :) Commented Aug 14 at 8:09
• @AlexTrounev Yeah, very mysterious. Perhaps I should try reporting this to WRI. Commented Aug 14 at 13:13

Before we answer all the questions, we note that not every numerical method is stable when solving the heat equation. Especially when we have boundary conditions in the form of ODE. In such cases, as is well known, it will be good to use implicit numerical method. Now we show 2 of them. First is clear FDM with implicit step

grid = Range[0, 1, 1/20]; dt = 1/100; nmax = 101; nn =
Length[grid]; U1 = ConstantArray[0, {nn, nmax}]; M1 =
NDSolveFiniteDifferenceDerivative[1, grid, "DifferenceOrder" -> 4][
"DifferentiationMatrix"]; M2 =
NDSolveFiniteDifferenceDerivative[2, grid, "DifferenceOrder" -> 4][
"DifferentiationMatrix"]; var1 = Table[Unique[u1], {nn}];
U1[[All, 1]] = 1 - grid^2;
Do[var = U1[[All, i]];
eq = (var1 - var)/dt - M2 . var1 - (M1 . var1)/grid // Quiet;
eq[[1]] = ((var1 - var)/dt - 2 M2 . var1)[[1]];
eq[[-1]] = var1[[-1]]; {vec, mat} = CoefficientArrays[eq, var1];
sol = LinearSolve[mat // N, -vec];
U1[[All, i + 1]] = sol;, {i, 1, nmax - 1}] // AbsoluteTiming


It takes about 0.1 s on my Silver Pentium. Second one is same as above, but we use NDSolve as follows

xgrid = Range[0, 1, 1/10]; M1 =
NDSolveFiniteDifferenceDerivative[1, xgrid, "DifferenceOrder" -> 4][
"DifferentiationMatrix"]; M2 =
NDSolveFiniteDifferenceDerivative[2, xgrid, "DifferenceOrder" -> 4][
"DifferentiationMatrix"];
U = Table[
uu[i][t], {i, Length[xgrid]}]; eq = -D[U, t] + M2 . U +
M1 . U/xgrid // Quiet;
eq[[1]] = (-D[U, t] + 2 M2 . U)[[1]];
eq[[-1]] = D[U, t][[-1]]; ini = U - (1 - xgrid^2) /. t -> 0;

sol2 = NDSolveValue[
Flatten@Table[{eq[[i]] == 0, ini[[i]] == 0}, {i, Length[eq]}],
U, {t, 0, 1}, StartingStepSize -> 10^-2,
Method -> {"LinearlyImplicitEuler"}]


Also we have two solutions demonstrated by @xzczd in a form of solfem and tst. We can show all solutions in one plot at r=0

sol1 = Interpolation[
Table[{(i - 1) dt, U1[[1, i]]}, {i, 1, nmax}]];

LogPlot[{sol1[t], sol2[[1]], solfem[t, 0], tst[t, 0]}, {t, 0, 1},
PlotStyle -> {{Blue}, {Orange, Dashed}, {Green}, {Blue, Dotted}},
PlotLegends -> {"FDM", "ImplicitEuler", "FEM", "MOL"}]


As we can see from the figure above, all 4 solutions coincide at small times, and then diverge so that solfem and tst solutions coincide, and sol1, sol2 also coincide. This is probably due to the fact that sol1 and sol2 have same order of accuracy. The same can be said about the other 2 solutions. Finally, we can discuss the "magic" of NDSolve. Let check

tst[[3]]

Out[]= {{0., 7.23317*10^-8, 1.44663*10^-7, 2.16995*10^-7,
2.89327*10^-7, 3.61659*10^-7, 5.06322*10^-7, 6.50985*10^-7,
7.95649*10^-7, 2.24228*10^-6, 3.68892*10^-6, 5.13555*10^-6,
0.0000196019, 0.0000340682, 0.0000485346, 0.000069204, 0.0000898734,
0.000110543, 0.000135035, 0.000159527, 0.000184019, 0.00022034,
0.000583552, 0.000946763, 0.00130997, 0.00207649, 0.00284301,
0.00360953, 0.00516498, 0.00672043, 0.00827588, 0.0109491,
0.0136224, 0.0162956, 0.020434, 0.0245724, 0.0287108, 0.034924,
0.0411371, 0.0473502, 0.0535634, 0.0625446, 0.0715258, 0.0805071,
0.0894883, 0.104579, 0.119669, 0.134759, 0.14985, 0.16985, 0.18985,
0.20985, 0.22985, 0.24985, 0.26985, 0.28985, 0.30985, 0.32985,
0.34985, 0.36985, 0.38985, 0.40985, 0.42985, 0.44985, 0.46985,
0.48985, 0.50985, 0.52985, 0.54985, 0.56985, 0.58985, 0.60985,
0.62985, 0.64985, 0.66985, 0.68985, 0.70985, 0.72985, 0.74985,
0.76985, 0.78985, 0.80985, 0.82985, 0.84985, 0.86985, 0.88985,
0.90985, 0.92985, 0.94985, 0.96985, 0.984925, 1.}, {0., 0.010101,
0.020202, 0.030303, 0.040404, 0.0505051, 0.0606061, 0.0707071,
0.0808081, 0.0909091, 0.10101, 0.111111, 0.121212, 0.131313,
0.141414, 0.151515, 0.161616, 0.171717, 0.181818, 0.191919, 0.20202,
0.212121, 0.222222, 0.232323, 0.242424, 0.252525, 0.262626,
0.272727, 0.282828, 0.292929, 0.30303, 0.313131, 0.323232, 0.333333,
0.343434, 0.353535, 0.363636, 0.373737, 0.383838, 0.393939,
0.40404, 0.414141, 0.424242, 0.434343, 0.444444, 0.454545, 0.464646,
0.474747, 0.484848, 0.494949, 0.505051, 0.515152, 0.525253,
0.535354, 0.545455, 0.555556, 0.565657, 0.575758, 0.585859, 0.59596,
0.606061, 0.616162, 0.626263, 0.636364, 0.646465, 0.656566,
0.666667, 0.676768, 0.686869, 0.69697, 0.707071, 0.717172, 0.727273,
0.737374, 0.747475, 0.757576, 0.767677, 0.777778, 0.787879,
0.79798, 0.808081, 0.818182, 0.828283, 0.838384, 0.848485, 0.858586,
0.868687, 0.878788, 0.888889, 0.89899, 0.909091, 0.919192,
0.929293, 0.939394, 0.949495, 0.959596, 0.969697, 0.979798,
0.989899, 1.}}


We see that first steps in time are about dt=7.23317*10^-8, while last steps are about 0.02. The implicit solver used in this computation is LSODA. Unfortunately, I was unable to use this solver to solve this problem.

• Why -1, maybe -500? :) Commented Aug 12 at 1:34
• (+1) This is a right answer. BTW, if we make the grid of sol1 and sol2 dense enough, say grid = Range[0, 1, 1/100]; dt = 1/2000; nmax = 2001; for sol1 and StartingStepSize -> 10^-3/2 together with xgrid = Range[0, 1, 1/100]; for sol2, the solutions will be visually the same. Commented Aug 12 at 4:00
• @xzczd Thank you for your observation. I have tested all implicit solvers implemented in NDSolve, and only Method -> {"LinearlyImplicitEuler"} works as expected. From NDSolveMethodData we can trace that solver used for this problem in your tst is LSODA. It is why first several steps are very small. But it looks like that we can't use LSODA as method of solution same way as NDSolve` with MOL option. :) Commented Aug 12 at 6:27
• (+51) I've also posted the solution I found as an answer, whose underlying idea is similar to yours. Have a look :) . Commented Aug 14 at 4:24