1
$\begingroup$

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?

$\endgroup$

1 Answer 1

5
$\begingroup$

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] 

enter image description here

Hope it helps!

$\endgroup$
9
  • $\begingroup$ nice use of DirichletCondition! $\endgroup$
    – Nasser
    Nov 13, 2023 at 15:56
  • $\begingroup$ @UlrichNeumann, wow, it is so simple and so confusing. Thank a lot. Could you please explain this feature? $\endgroup$ Nov 13, 2023 at 16:01
  • $\begingroup$ @Nasser Thanks! The alternative ...Method -> {"MethodOfLines", "TemporalVariable" -> t, "SpatialDiscretization" -> {"FiniteElement"}} was too long for me ;-) $\endgroup$ Nov 13, 2023 at 16:07
  • $\begingroup$ @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"}} $\endgroup$ Nov 13, 2023 at 16:10
  • 1
    $\begingroup$ Try sol["ElementMesh"] for your first example , surprisingly you might see a FiniteElement Mesh. $\endgroup$ Nov 13, 2023 at 18:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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