Semi-classical approximation of modified Schrodinger's equation on a sphere leads to ndsz warning

Question

I know this error is very common. I checked out a good number of proposed solutions to this problem, but unfortunately, none of them helped in my case.

I would appreciate any help or hints on how to overcome this issue.

My program is quite simple. It is supposed to solve two coupled time-dependent partial differential equations on a sphere. I don't know if it's relevant, but from the physical point of view, it's actually the semi-classical approximation of a modified Schrodinger's equation on a sphere, with the Newtonian gravitational potential made scale-invariant.

The program is supposed to solve the PDE up to t=100, but it breaks down even at a value as low as 0.01. Here's the code with some comments for further information:

  (* defining the potential that will be used in the equations. *)
V[ϕ, θ] := 
  0.045360921162651446/Sqrt[1 + Cos[ϕ] Sin[θ]] + 
   0.045360921162651446/Sqrt[
   1 - Sin[θ] Sin[π/6 - ϕ]] + 0.045360921162651446/
   Sqrt[1 - Sin[θ] Sin[π/6 + ϕ]];
(* defining the region in which the differential equation is to be solved. θ = π/2 is avoided on purpose because of the singularities of V *)
Π1 = 
  ImplicitRegion[
   0 <= ϕ <= 2 π && 
    0 <= θ <= π/2 - 0.001, {ϕ, θ}];
Π2 = 
  ImplicitRegion[
   0 <= ϕ <= 
     2 π && π/2 + 
      0.001 <= θ <= π, {ϕ, θ}];
Π = RegionUnion[Π1, Π2];
(* defining the two coupled PDEs. *)
firsteq = {-D[Subscript[S, 0][ϕ, θ, t], t] == -V[ϕ, θ] + (Derivative[0, 1, 0][Subscript[S, 0]][ϕ, θ, t]^2 + 
       Sec[θ]^2Derivative[1, 0, 0][Subscript[S, 0]][ϕ, θ, t]^2)}; 
secondeq = {-D[Subscript[S, 1][ϕ, θ, t], t] == -((-Tan[θ])Derivative[0, 1, 0][Subscript[S, 0]][ϕ, θ, t] + 
        Derivative[0, 2, 0][Subscript[S, 0]][ϕ, θ, t] + Sec[θ]^2Derivative[2, 0, 0][Subscript[S, 0]][ϕ, θ, t]) + 
      2(Derivative[0, 1, 0][Subscript[S, 0]][ϕ, θ, t]Derivative[0, 1, 0][Subscript[S, 1]][ϕ, θ, t] + 
        Sec[θ]^2Derivative[1, 0, 0][Subscript[S, 0]][ϕ, θ, t]*Derivative[1, 0, 0][Subscript[S, 1]][ϕ, θ, t])};
(* defining the initial conditions *)
initial = {Subscript[S, 0][ϕ, θ, 0] == Sin[θ], 
   Subscript[S, 1][ϕ, θ, 0] == Sin[θ]^2};
Equations = Join[firsteq, secondeq, initial];
(* Solving the equations ... *)
Solution = 
  NDSolve[Equations, {Subscript[S, 0], Subscript[S, 1]}, {t, 0, 
    100}, {ϕ, θ} ∈ Π, 
   Method -> "Automatic"];

I get the error:

NDSolve::ndsz: At t == 0.011684614824200409`, step size is effectively zero; singularity or stiff system suspected.

The system is indeed singular. I plotted the incomplete solution and it appears that very soon the solution diverges over some regions near the poles and equator. Is there some way to overcome this problem?

Thank you very much.

Edit:

I modified the potential to remove its singularities:

V[ϕ_, θ_] := 
      0.045360921162651446/Sqrt[1.00001 + Cos[ϕ] Sin[θ]] + 
       0.045360921162651446/Sqrt[
       1.00001 - Sin[θ] Sin[π/6 - ϕ]] + 0.045360921162651446/
       Sqrt[1.00001 - Sin[θ] Sin[π/6 + ϕ]];

Following the suggestions in the comments, I tried to correct another issue in the code and define periodic boundary conditions for ϕ. Here's the code I added. I'm grateful to @xzczd for this.

  δ = 10^-5;
periodicity = {PeriodicBoundaryCondition[
    Subscript[S, 0][ϕ, θ, t], ϕ == 0 && 
     0 <= θ <= π , Function[x, x + {2 π, 0}]], 
   PeriodicBoundaryCondition[
    Subscript[S, 0][ϕ, θ, 
     t], ϕ == 2 π + δ && 0 <= θ <= π , 
    Function[x, x + {-2 π, 0}]], 
   PeriodicBoundaryCondition[
    Subscript[S, 1][ϕ, θ, t], ϕ == 0 && 
     0 <= θ <= π , Function[x, x + {2 π, 0}]], 
   PeriodicBoundaryCondition[
    Subscript[S, 1][ϕ, θ, 
     t], ϕ == 2 π + δ && 0 <= θ <= π , 
    Function[x, x + {-2 π, 0}]]};

But this time, I got this error:

NDSolve::ndcf: Repeated convergence test failure at t == 0.`; unable to continue.

Following @xzczd's helpful suggestion, I decided to forget about FiniteElements and use the TensorProduct method instead. Here's the final code. I also changed the initial conditions.

    firsteq = {-D[Subscript[S, 0][ϕ, θ, t], t] == -V[ϕ, θ] + (Derivative[0, 1, 0][Subscript[S, 0]][ϕ, θ, t]^2 + Csc[θ]^2*Derivative[1, 0, 0][Subscript[S, 0]][ϕ, θ, t]^2)}; 
    secondeq = {-D[Subscript[S, 1][ϕ, θ, t], t] == -((-Cot[θ])*Derivative[0, 1, 0][Subscript[S, 0]][ϕ, θ, t] + Derivative[0, 2, 0][Subscript[S, 0]][ϕ, θ, t] + Csc[θ]^2*Derivative[2, 0, 0][Subscript[S, 0]][ϕ, θ, t]) + 
          2*(Derivative[0, 1, 0][Subscript[S, 0]][ϕ, θ, t]*Derivative[0, 1, 0][Subscript[S, 1]][ϕ, θ, t] + Csc[θ]^2*Derivative[1, 0, 0][Subscript[S, 0]][ϕ, θ, t]*Derivative[1, 0, 0][Subscript[S, 1]][ϕ, θ, t])};
initial = {Subscript[S, 0][ϕ, θ, 0] == Sin[θ]^2, 
           Subscript[S, 1][ϕ, θ, 0] == Sin[θ]^2};



Periodicity = {Subscript[S, 0][0, θ, t] == 
    Subscript[S, 0][2 π, θ, t], 
   Subscript[S, 1][0, θ, t] == 
    Subscript[S, 1][2 π, θ, t]};
        Equations = Join[firsteq, secondeq, initial, periodicity];

        Solution = 
              NDSolve[Equations, {Subscript[S, 0], Subscript[S, 1]}, {t, 0, 
                10}, {ϕ, 0, 2 π}, {θ, 0, π}, 
               Method -> {"PDEDiscretization" -> {"MethodOfLines", 
                   "SpatialDiscretization" -> {"TensorProductGrid", 
                     "MaxPoints" -> 500, "MinPoints" -> 500}}}];

However, after about half an hour, I got this error:

NDSolve::ndsz: At t == 0.00003716693557750961`, step size is effectively zero; singularity or stiff system suspected.
NDSolve::eerr: Warning: scaled local spatial error estimate of 1.0182820651808292*^12 at t = 0.00003716693557750961 in the direction of independent variable ϕ is much greater than the prescribed error tolerance. Grid spacing with 500 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

I tried to make the discretization finer, but the kernel crashed for MaxPoints, MinPoints = 1000.

Answer 1

First of all, I'd like to emphasize again that it's important to make sure the equation system itself is correct. Anyway, after the correction the system is less demanding, so let me post an answer.

The idea is simple: the system is nonlinear, derivative of $S_0$ in $\theta$ direction of firsteq is of 1st order, derivative of $S_1$ in $\theta$ direction of secondeq is of 1st order, the solution of the system seems to be stiff, so the experience obtained in e.g. this post may help. Let's try artificial viscosity.

After adding artificial viscosity, something related to the issue discussed in this post seems to show up, so I use pdetoode to discretize the system.

V[ϕ_, θ_] := 
  0.045360921162651446/Sqrt[1.00001 + Cos[ϕ] Sin[θ]] + 
   0.045360921162651446/Sqrt[1.00001 - Sin[θ] Sin[π/6 - ϕ]] + 
   0.045360921162651446/Sqrt[1.00001 - Sin[θ] Sin[π/6 + ϕ]];
Clear[points, grid, domain]
points@ϕ = points@θ = 200;
eps = 0(Pi/10^2);
domain@ϕ = {0, 2 Pi}; domain@θ = {eps, Pi - eps};
(grid@# = Array[# &, points@#, domain@#]) & /@ {ϕ, θ};
difforder = 2;
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = 
  pdetoode[{Subscript[S, 0], Subscript[S, 1]}[ϕ, θ, t], t, 
   grid /@ {ϕ, θ}, difforder, {True, False}];
μ0 = 10^-2; μ1 = 10^-2;
Unevaluated[
   eq = {μ0 (D[Subscript[S, 0], {θ, 2}] + 
          D[Subscript[S, 0], {ϕ, 2}]) - 
       D[Subscript[S, 0], 
        t] == -V[ϕ, θ] + (D[Subscript[S, 0], θ]^2 + 
         Csc[θ]^2 D[Subscript[S, 0], ϕ]^2), μ1 (D[Subscript[S, 
           1], {θ, 2}] + D[Subscript[S, 1], {ϕ, 2}]) - 
       D[Subscript[S, 1], 
        t] == -(-Cot[θ] D[Subscript[S, 0], θ] + 
          D[Subscript[S, 0], {θ, 2}] + 
          Csc[θ]^2 D[Subscript[S, 0], {ϕ, 2}]) + 
       2 (D[Subscript[S, 0], θ] D[Subscript[S, 1], θ] + 
          Csc[θ]^2 D[Subscript[S, 0], ϕ] D[Subscript[S, 1], ϕ])};
   ic = {Subscript[S, 0] == Sin[θ]^2, Subscript[S, 1] == Sin[θ]^2} /. 
     t -> 0;
   bc = {D[Subscript[S, 0], θ] == 0, 
      D[Subscript[S, 1], θ] == 
       0} /. {{θ -> domain[θ][[1]]}, {θ -> 
        domain[θ][[-1]]}}] /. {Subscript[S, 0] -> 
    Subscript[S, 0][ϕ, θ, t], 
   Subscript[S, 1] -> Subscript[S, 1][ϕ, θ, t]};
del = #[[2 ;; -2]] &;
ode = Map[del, ptoofunc@eq, {2}] // Quiet;
odebc = ptoofunc@diffbc[t][bc];
odeic = ptoofunc@ic;
varlst = Outer[#@##2 &, {Subscript[S, 0], Subscript[S, 1]}, grid@ϕ, 
   grid@θ];
tend = 10;
showStatus[status_]:=LinkWrite[$ParentLink,
  SetNotebookStatusLine[FrontEnd`EvaluationNotebook[],ToString[status]]];
clearStatus[]:=showStatus[""];
clearStatus[]
jianshi[t_]:=EvaluationMonitor:>showStatus["t = "<>ToString[CForm[t]]]
 sollst = NDSolveValue[{ode, odeic, odebc}, varlst, {t, 0, tend}, jianshi[t], 
   Method -> {"EquationSimplification" -> "MassMatrix"}(,SolveDelayed->
   True)]; // AbsoluteTiming
(* {1984.24, Null} *)
{sols0, sols1} = rebuild[#, grid /@ {ϕ, θ}, -1] & /@ sollst;

Remark

Adding PrecisionGoal -> 3, AccuracyGoal -> 3 to NDSolveValue will reduce the timing to about 277.04 seconds! But the beginning part of $S_1$ (for about $t<0.2$) changes. (Not widely different, but quite observable. ) Anyway, if the accuracy of transient behavior isn't important for you, consider adding these options.

The showStatus is from this post.

I've solved $S_0$ and $S_1$ all in one, but since firsteq only involves $S_0$ (and further trial suggests it's easier to solve than secondeq), you may consider solve firsteq first.

Let's check the solution. $S_0$:

piclst0 = Table[
   Plot3D[sols0[f, th, t], {th, #3, #4}, {f, #1, #2}, PlotPoints -> 50] & @@ 
    Flatten[domain /@ {ϕ, θ}], {t, 0, tend, tend/25}];
ListAnimate[piclst0]

$S_1$:

piclst1 = Table[
   SphericalPlot3D[sols1[f, th, t], 
     {th, #3, #4}, {f, #1, #2}, PlotPoints -> 50] & @@
       Flatten[domain /@ {ϕ, θ}], {t, 0, tend, tend/25}];
ListAnimate[
 Show[#, PlotRange -> {{-1.2, 1.2}, {-1.2, 1.2}, {-3.2, 3.2}}] & /@ piclst1]

The solution around $\theta=\pi/2$ looks a bit unusual, I'm not sure if it's the nature of the solution or numeric error. You may adjust points@ϕ, points@θ, μ0, μ1 further to see how the solution changes.