1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Nasser Jan 21 '17 at 22:08
4
$\begingroup$

Edit from the author of the answer

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]]

enter image description here

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

Edit

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

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]]

enter image description here

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]&

enter image description here

$\endgroup$
  • $\begingroup$ That really works, thanks. So, the option opts corrected the error? $\endgroup$ – Asatur Khurshudyan Jan 22 '17 at 19:29
  • $\begingroup$ Are there any chances to use RevolutionPlot3D? Because I can`t make any conclusion relying on ContourPlot. $\endgroup$ – Asatur Khurshudyan Jan 22 '17 at 19:38
  • $\begingroup$ I have never used RevolutionPlot3D. In all cases, there are solutions with ContourPlot3D + manual coordinates change. (ContourPlot3D is often slow) $\endgroup$ – andre314 Jan 22 '17 at 19:41
  • $\begingroup$ 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.". $\endgroup$ – Asatur Khurshudyan Jan 22 '17 at 20:17
  • $\begingroup$ The final question, what does res[[1, 1, 2]] mean? And why I see a white line in the middle of the circle? $\endgroup$ – Asatur Khurshudyan Jan 22 '17 at 21:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.