# 1/0 infinity while solving heat equation in polar coordinates

The ndsolve

    uo = 1; m = 1; n = 1;  a = 1; T = 20;  M = 10;
W0[r_, phi_] := Exp[-(1 - r)^2]
ud[phi_] := UnitStep[phi] - UnitStep[phi - Pi/36]
u[t_] := uo (1 - Cos[(2 Pi m)/T t])^n

NDSolve[{a Laplacian[W[r, phi, t], {r, phi},"Polar"]== D[W[r, phi, t],t], (D[W[r, phi, t], r] /. r -> 1) == ud[phi] u[t], W[r, phi,0] ==W0[r, phi]}, W, {r, 0, 1}, {phi, 0, 2 Pi}, {t, 0, T}]


gives the error

Power::infy: Infinite expression 1/0.^2 encountered.
NDSolveValue::ndnum: Encountered non-numerical value for a derivative at t == 0..


No idea on how to proceed. Any hints or ideas?

• When solving this by hand, we know that in heat equation, on a disk, one of the conditions must be that u is bounded at r=0 (center of disk). This causes one of the solutions causing the blow up to be discarded. This is standard method when solving heat PDE on disk. It looks like ndsolve is complaining about this in your example. – Nasser Jan 21 '17 at 22:08

Caution : There something strange in this solution. Not sure it's OK. The strange thing is that the temperature tends towards 0. This is unlikely, except if the Op's target is precisely to do that.
To be analysed.

Here is the solution :

I works only on Mathematica >= 11

uo = 1; m = 1; n = 1;  a = 1; T = 20;  M = 10;
W0[r_, phi_] := Exp[-(1 - r)^2]
ud[phi_] := UnitStep[phi] - UnitStep[phi - Pi/36]
u[t_] := uo (1 - Cos[(2 Pi m)/T t])^n

opts = Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}};

res=NDSolve[{
a Laplacian[W[r, phi, t], {r, phi},"Polar"]== D[W[r, phi, t],t]+NeumannValue[ud[phi] u[t],r==1],
W[r, phi,0] ==W0[r, phi]
}, W, {r, 0, 1}, {phi, 0, 2 Pi}, {t, 0, T},
opts]

ContourPlot[ res[[1,1,2]][Sqrt[x^2+y^2],Pi+ArcTan[-x,-y], 10],{x,-1,1},
{y,-1,1},RegionFunction->Function[{x,y},x^2+y^2<0.99]] My approach :

• As I was aware that the 1/0 problem was solved on V11 and fem (See the penultimate paragraph here.), I added the option FiniteElement
• Then I had a new problem with a error message that let me think that your condition (D[W[r, phi, t], r] /. r -> 1) == ud[phi] u[t], which is a Neumann boundary conditon, was interpreted as a Dirichlet boundary condition (I don't know why), so i tried the explicit Neumann syntax and it worked

The solution above corresponds to the case where there is a adiabatic wall at phi=0. Here is the solution without the wall :

Edit from the Author, 01/01/2020 The following is not correct : The wall is not completely removed. The problem is that the flux that enters the wall is not the same than the flux that get outside the wall. A problem of this kind is solved here
END EDIT

uo = 1; m = 1; n = 1;  a = 1; T = 20;  M = 10;
W0[r_, phi_] := Exp[-(1 - r)^2]
ud[phi_] := UnitStep[phi] - UnitStep[phi - Pi/36]
u[t_] := uo (1 - Cos[(2 Pi m)/T t])^n

opts = Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}};

res=NDSolveValue[{
a Laplacian[W[r, phi, t], {r, phi},"Polar"]== D[W[r, phi, t],t]+NeumannValue[ud[phi] u[t],r==1],
PeriodicBoundaryCondition[W[r,phi,t],phi==2 Pi,Function[v,v+{0,-2 Pi}]],
W[r, phi,0] ==W0[r, phi]
}, W, {r, 0, 1}, {phi, 0, 2 Pi}, {t, 0, T},
opts]

ContourPlot[ res[Sqrt[x^2+y^2],Pi+ArcTan[-x,-y], 0.3],{x,-1,1},
{y,-1,1},RegionFunction->Function[{x,y},x^2+y^2<0.99]] An animation of the same thing :

Labeled[ContourPlot[ res[Sqrt[x^2+y^2],Pi+ArcTan[-x,-y], #],{x,-1,1},
{y,-1,1},RegionFunction->Function[{x,y},x^2+y^2<0.99],
ContourLabels->True,Contours-> {
.9,.8613,.861,.860,.858,.856,.855,.850,.845,.840,.835,.830,
.825,.8,.775,.750,.725,.7,.68,.65,.625,.6,.55,.525,
.5,.475,.45,.425,.4,.38,.36,.34,.32,.30,.26,.24,.20,
.175,.15,.125,.10,.075,.05,
.020,1. 10^-3,1. 10^-8},
ColorFunction-> "Temperature",
ColorFunctionScaling->False
],Row[{"t=",#}]]& /@ Range[0,16,0.1] //
Export["therm00.gif",#, "DisplayDurations" -> 0.3]& • That really works, thanks. So, the option opts corrected the error? – Asatur Khurshudyan Jan 22 '17 at 19:29
• Are there any chances to use RevolutionPlot3D? Because I cant make any conclusion relying on ContourPlot. – Asatur Khurshudyan Jan 22 '17 at 19:38
• I have never used RevolutionPlot3D. In all cases, there are solutions with ContourPlot3D + manual coordinates change. (ContourPlot3D is often slow) – andre314 Jan 22 '17 at 19:41
• When I write RevolutionPlot3D[W0[r, phi], {r, 0, 1}, {phi, 0, 2 Pi}], it plotes what is expected. However, if I want to plot RevolutionPlot3D[W[r, phi, T] /. res, {r, 0, 1}, {phi, 0, 2 Pi}], it says "Dot::rect: Nonrectangular tensor encountered.". – Asatur Khurshudyan Jan 22 '17 at 20:17
• The final question, what does res[[1, 1, 2]] mean? And why I see a white line in the middle of the circle? – Asatur Khurshudyan Jan 22 '17 at 21:19