# Numerical solving diffusion equation in spherical coordinates

Mathematica nicely solves Poisson's equation in spherical coordinates as

eqn=Laplacian[V[r,\[Theta]],{r,\[Theta],\[Phi]},"Spherical"]==-Sin[\[Theta]];
a=1;b=10;
sol=NDSolveValue[{eqn,V[a,\[Theta]]==1,V[b,\[Theta]]==0},V,{r,a,b},{\[Theta],0,\[Pi]}];


but failed to solve a similar time-dependent problem

eqn = D[V[t, r, \[Theta]], t] ==
Laplacian[V[t, r, \[Theta]], {r, \[Theta], \[Phi]}, "Spherical"] +
Sin[\[Theta]];
a = 1/2; b = 1;
sol = NDSolveValue[{eqn, V[t, a, \[Theta]] == 1,
V[t, b, \[Theta]] == 0, V[0, r, \[Theta]] == (r - b)/(a - b)},
V, {t, 0, 1}, {r, a, b}, {\[Theta], 0, \[Pi]}];


Why does Mathematica complains for this case?

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

NDSolveValue::ndnum: Encountered non-numerical value for a derivative at t == 0..

I know that the pole is a special point. However, how did Mathematica get around that in the first case?

I got the answer (see below) which means I have to call explicitly FiniteElement. However, I've faced with even more strange 'by default' behavior. Here they are two independent equations which can be solved without problem

eqn = -D[V[t, r, \[Theta]], t] +
Laplacian[V[t, r, \[Theta]], {r, \[Theta], \[Phi]}, "Spherical"] +
Cos[\[Theta]];
eqn1 = -D[W[t, r, \[Theta]], t] +
Laplacian[W[t, r, \[Theta]], {r, \[Theta], \[Phi]}, "Spherical"] +
Cos[\[Theta]];
a = 1/2; b = 1;
sol = NDSolveValue[{eqn == 0,
eqn1 == NeumannValue[-W[t, r, \[Theta]]/r, r == b],
DirichletCondition[W[t, r, \[Theta]] == 0, r == a],
W[0, r, \[Theta]] == 0,
DirichletCondition[V[t, r, \[Theta]] == 0, r == a],
DirichletCondition[V[t, r, \[Theta]] == 0, r == b],
V[0, r, \[Theta]] == 0}, {V, W}, {t, 0, 1}, {r, a, b}, {\[Theta],
0, \[Pi]}];


But if you slightly change an order of equations (eqn == 0 moved to 4th position)

sol = NDSolveValue[{eqn1 ==
NeumannValue[-W[t, r, \[Theta]]/r, r == b],
DirichletCondition[W[t, r, \[Theta]] == 0, r == a],
W[0, r, \[Theta]] == 0, eqn == 0,
DirichletCondition[V[t, r, \[Theta]] == 0, r == a],
DirichletCondition[V[t, r, \[Theta]] == 0, r == b],
V[0, r, \[Theta]] == 0}, {V, W}, {t, 0, 1}, {r, a, b}, {\[Theta],
0, \[Pi]}];


then you will get an error

LinearSolve::parpiv: Zero pivot was detected during the numerical factorization or there was a problem in the iterative refinement process. It is possible that the matrix is ill-conditioned or singular.

I'm puzzled what is a rule or it is a bug?

Change your boundary conditions to DirichletConditions

eqnt = D[V[t, r, \[Theta]], t] ==
Laplacian[V[t, r, \[Theta]], {r, \[Theta], \[Phi]}, "Spherical"] +Sin[\[Theta]] ;
a = 1/2; b = 1;

sol = NDSolveValue[{ eqnt,
DirichletCondition[V[t, r, \[Theta]] == 1, r == a],
DirichletCondition[V[t, r, \[Theta]] == 0, r == b] ,
V[0, r, \[Theta]] == (r - b)/(a - b)},
V, {t, 0, 1}, {r, a, b}, {\[Theta], 0, \[Pi]}]


time dependeent solutions

Plot3D[Evaluate@Table[sol[t, r, \[Theta]], {t, {0, .1, 1  }/5}], {r,a, b}, {\[Theta], 0, Pi}, PlotStyle -> Opacity[1], Mesh -> False]


Hope it helps!

• nice use of DirichletCondition! Nov 13, 2023 at 15:56
• @UlrichNeumann, wow, it is so simple and so confusing. Thank a lot. Could you please explain this feature? Nov 13, 2023 at 16:01
• @Nasser Thanks! The alternative ...Method -> {"MethodOfLines", "TemporalVariable" -> t, "SpatialDiscretization" -> {"FiniteElement"}} was too long for me ;-) Nov 13, 2023 at 16:07
• @RodionStepanov Using DirichletConditionsor NeumannValue forces NDSolve to call Finite Element solver. Instead you could get the same result if you use NDSolve without DirchletCondition and add the option Method -> {"MethodOfLines", "TemporalVariable" -> t, "SpatialDiscretization" -> {"FiniteElement"}} Nov 13, 2023 at 16:10
• Try sol["ElementMesh"]` for your first example , surprisingly you might see a FiniteElement Mesh. Nov 13, 2023 at 18:05