# Solving heat equation in spherical coordinates

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?

• what happens to your pde when $r=0$ ? Commented Mar 15 at 7:59

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"}]]
`