Solving heat equation in spherical coordinates

Question

I am attempting to solve a heat equation in spherical coordinates. However, I encountered an error message. My Mathematica code is

eqn = D[u[r, θ, ϕ, t], t] == 1/r^2 D[r^2 D[u[r, θ, ϕ, t], r], r] +    
      1/(r^2 Sin[θ]) D[Sin[θ] D[u[r, θ, ϕ, t], θ], θ] + 
      1/(r^2 Sin[θ]^2) D[u[r, θ, ϕ, t], {ϕ, 2}];

sol = NDSolveValue[{eqn, u[1, θ, ϕ, t] == 0,
  (D[u[r, θ, ϕ, t], r]/. r -> 0) == 0,
  (D[u[r, θ, ϕ, t], θ] /. θ -> 0) == 0,
  (D[u[r, θ, ϕ, t], θ] /. θ -> π) == 0,
  u[r, θ, ϕ, 0] == Exp[-r^2]},
  u, {r, 0, 1}, {θ, 0, π}, {ϕ, 0, 2 π}, {t, 0, 1}]

The error message is

"Power::infy: Infinite expression 1/0.^2 encountered."

However, the integral ranges should be [r, 0, 1], [θ, 0, π], and [φ, 0, 2 π] to cover the equivalent volume. Is there a way to solve this issue?

Answer 1

Try correct Laplacian(?) and "MethodOfLines" without the boundary conditions D[u[...],r],u[...],\[Theta]] (which are probably NeumannValue's).

I changed the DirichletCondition at r==1 to DirichletCondition[u[t, r, \[Theta], \[Phi] ] == Exp[-t/.001] to make the condition consistent to the initial condition!

Boundary condition \[Phi] == 0and \[Phi] ==2Pi are PeriodicBoundaryConditions!

eqn = D[u[r, \[Theta], \[Phi], t], t] == 
Laplacian[u[r, \[Theta], \[Phi], t], {r, \[Theta], \[Phi]},"Spherical"] ;


    U = NDSolveValue[{eqn,
   DirichletCondition[u[ r, \[Theta], \[Phi], t ] == Exp[-t/.001], 
    r == 1 && 0 < \[Phi] < 2 Pi], 
   u[ r, \[Theta], \[Phi] , 0] == Exp[-r^2],
   PeriodicBoundaryCondition[u[ r, \[Theta], \[Phi], t ], \[Phi] == 0,
     TranslationTransform[{0, 0, 2 Pi}]], 
   PeriodicBoundaryCondition[
    u[ r, \[Theta], \[Phi], t ], \[Phi] == 2 Pi, 
    TranslationTransform[{0, 0, -2 Pi}]]}, 
  u, {r, 0, 1}, {\[Theta], 0, \[Pi]}, {\[Phi], 0, 2 \[Pi]}, {t, 0, 1},
   Method -> {"MethodOfLines", "TemporalVariable" -> t, 
    "SpatialDiscretization" -> "FiniteElement"}]

Unfortunately `NDSolve doesn't evaluate and gives the error message "LinearSolve::sing: Matrix SparseArray[<2026094>, {41343, 41343}] is singular." and doesn't finish evaluation.

No idea what's wrong here!

Workaround:

In cartesian coordinates the simulation runs without problems:

region = Ball[]
U = NDSolveValue[{Derivative[1, 0, 0, 0][u][t, x, y, z] == 
    Laplacian[u[t, x, y, z], {x, y, z}],
   DirichletCondition[u[t, x, y, z ] == Exp[-t/.001], 
    x^2 + y^2 + z^2 == 1], 
   u[0, x, y, z ] == Exp[-(x^2 + y^2 + z^2)^2]
    }, u, Element[{x, y, z}, region], {t, 0, 1}, 
  Method -> {"MethodOfLines", "TemporalVariable" -> t, 
    "SpatialDiscretization" -> "FiniteElement"}]

    Block[{}, 
 Manipulate[
  SliceDensityPlot3D[U[t, x, y, z], "CenterPlanes", 
   Element[{x, y, z}, region], ColorFunction -> "TemperatureMap", 
   ColorFunctionScaling -> ! True, 
   RegionFunction -> Function[{x, y, z}, x^2 + y^2 + z^2 < 1]] , {{t, 
    0}, 0, 1, Appearance -> "Labeled"}]]