7
$\begingroup$

I am trying to solve this eigenvalue problem:
\begin{align} \mu \Psi(r) & = -\frac{1}{2}\left ( \Psi^{\prime \prime}(r) + \frac{2}{r} \Psi' (r)\right ) -4\pi \Psi(r) \int _0^\infty dr' r'^2 \frac{\Psi(r')^2}{r>}, \end{align} where $\mu$ is the eigenvalue, $\Psi(r)$ the eigenfunction I'm trying to solve, and $r_>$ is the greater one of $r$ and $r'$. The requirement of the system is
\begin{align}\Psi'(0)& = 0, \cr \Psi(\infty) & = 0.\end{align}

Since it involves an integral, the way I deal with it is to decouple it to two equations.

\begin{align} \mu \Psi (r) & = -\frac{1}{2}\left ( \Psi^{\prime \prime}(r) + \frac{2}{r} \Psi'(r) \right ) +\Psi(r) \Phi(r), \\ \nabla^2 \Phi(r) & = 4\pi \Psi(r)^2, \end{align}

with boundary

\begin{align} \Psi'(0)& = 0, \\ \Psi(\infty) & = 0,\\ \Phi'(0)& = 0, \\ \Phi(\infty) & = 0. \end{align}

Any idea how to solve it?

Attempt 1:

I manually tweak the boundary condition $\Psi(0)$ at zero trying to find a solution of $\Psi(r)$ that doesn't blow up at infinity. \begin{align} \mu \Psi (r) & = -\frac{1}{2}\left ( \Psi^{\prime \prime}(r) + \frac{2}{r} \Psi'(r) \right ) +\Psi(r) \Phi(r), \\ \nabla^2 \Phi(r) & = 4\pi \Psi(r)^2, \end{align}

with boundary

\begin{align} \Psi'(0)& = 0, \\ \Psi(0) & = A,\\ \Phi'(0)& = 0, \\ \Phi(\infty) & = 0. \end{align}

My current method is to guess a pair of input data $(\mu, A)$, use NDSolve to solve for $\Psi, \Phi$, and check whether the solution gives me a $\Psi$ whose absolute value decreases monotonically $i.e.$ vanishes at infinity, as suggested here. However, even with this method I haven't had much success. So, a) is there a good way to implement this algorithm; b) is there a better/alternative way to attack this whole problem? --AFAIK, (N)DEigensystem cannot handle this problem.

Edited:

Attemp 2:

So I just tried it out, naively if I use NDEigensystem directly as follows, it won't solve at all, which is not surprising.

rStart1=10^-3;
rEnd1=5;
epsilon=0;
sollst1=NDEigensystem[{(-1/2*(D[ψ[r],r,r]+2/r*D[ψ[r],r])-4Pi*ψ[r]*Integrate[rp^2*ψ[rp]^2/If[rp>r,rp,r],{rp,0,Infinity}])+NeumannValue[0,r==rStart1]},ψ[r],{r,rStart1,rEnd1},8];//AbsoluteTiming

Attemp 3: So this time I have fixed $\mu$, kept the boundary condition as \begin{align} \Psi'(0)& = 0, \\ \Psi(\infty) & = 0,\\ \Phi'(0)& = 0, \\ \Phi(\infty) & = 0. \end{align} and used the shooting method to select normalization for $\Psi(0)$. The following code works but only for very specific choice of rEnd1 and $\mu=-0.4$. Changing any of them gives me a trivial solution again. It seems the code is a bit fine-tuned in this sense.

rEnd1 = 3.7;
rStart1 = 10^-3;
stableHunter3[μ_] := Module[{}, epsilon = 0;
  eqn2 = {-μ*ψ[r] - 
      1/2*(D[ψ[r], r, r] + 2/r*D[ψ[r], r]) + ϕ[
        r]*ψ[r](*-1/8*ψ[r]^3*)== 0, 
    D[ϕ[r], r, r] + 2/r*D[ϕ[r], r] == 4 Pi*ψ[r]^2};
  bc2 = {ψ[rEnd1] == 
     epsilon, (D[ψ[r], r] /. r -> rStart1) == 
     epsilon, (D[ϕ[r], r] /. r -> rStart1) == 
     epsilon, (ϕ[rEnd1]) == epsilon};
  sollst1 = 
   Map[NDSolveValue[
      Flatten@{eqn2, bc2}, {ψ[r], ϕ[r]}, {r, rStart1, 
       rEnd1}, Method -> 
       "BoundaryValues" -> {"Shooting", 
         "StartingInitialConditions" -> {ψ[0] == #}}, 
      Method -> "StiffnessSwitching"] &, Range[-3, 3, 0.1]]
  ]
funclst = stableHunter3[-0.4];
Plot[Evaluate[funclst /. r -> r00], {r00, rStart1, rEnd1}]

Attemp 4: So based on the suggestion from @bbgodfrey, an analysis on the asymptotic behavior suggests the following. When $r\rightarrow \infty$, \begin{align} \mu \Psi(r) & = -\frac{1}{2}\left ( \Psi^{\prime \prime}(r) + \frac{2}{r} \Psi' (r)\right ) -4\pi \Psi(r) \int _0^\infty dr' r'^2 \frac{\Psi(r')^2}{r>}, \cr & \approx -\frac{1}{2}\left ( \Psi^{\prime \prime}(r) + \frac{2}{r} \Psi' (r)\right ) -\frac{N}{r}\Psi(r),\end{align} where $N\equiv \int_0^\infty \Psi(r)^2 4\pi r^2 d r$. Here $\Phi(r)$ is not necessary, but if it is defined then $\Phi(r) \approx -\frac{N}{r}$ at large $r$. Then, I could use the following code to solve the eigenvalue $\mu$.

rStart1 = 10^-3;
rEnd1 = 5;
epsilon = 0;
Nphy1 = 1;

sollst1 = 
   NDEigensystem[{(-1/2*(D[ψ[r], r, r] + 2/r*D[ψ[r], r]) - 
        Nphy1/r*ψ[r]) + NeumannValue[0, r == rStart1](*,
     DirichletCondition[ψ[r]\[Equal]epsilon,
     x\[Equal]rStart1]*)}, ψ[r], {r, rStart1, rEnd1}, 
    8]; // AbsoluteTiming
sollst1
sollst1[[2]] /. r -> rEnd1
Plot[Evaluate[%[[2]]], {r, rStart1, rEnd1}, 
 PlotRange -> All(*,PlotLabels\[Rule]{1,2,3}*)]


{{-0.176753, 0.497306, -0.507498, 1.62112, 3.15356, 5.08955, 7.42803, 
10.1706}, {InterpolatingFunction[{{0.001, 5.}}, <>][r], 
  InterpolatingFunction[{{0.001, 5.}}, <>][r], 
  InterpolatingFunction[{{0.001, 5.}}, <>][r], 
  InterpolatingFunction[{{0.001, 5.}}, <>][r], 
  InterpolatingFunction[{{0.001, 5.}}, <>][r], 
  InterpolatingFunction[{{0.001, 5.}}, <>][r], 
  InterpolatingFunction[{{0.001, 5.}}, <>][r], 
  InterpolatingFunction[{{0.001, 5.}}, <>][r]}}
{0.205938, -0.112739, 0.0259637, -0.094238, -0.0840003, -0.0769606}

enter image description here

I believe the $\mu =-0.507498$ is the ground state eigenvalue. However, when I use the boundary conditions at $rEnd1$ to trade for the ones at infinity, I ended up with trivial solution only. 

stableHunter3[μ_] := Module[{},
  eqn2 = {-μ*ψ[r] - 
      1/2*(D[ψ[r], r, r] + 2/r*D[ψ[r], r]) + ϕ[
        r]*ψ[r] == 0, 
    D[ϕ[r], r, r] + 2/r*D[ϕ[r], r] == 4 Pi*ψ[r]^2};
  bc2 = {ψ[rEnd1] == 
     0.025963699315910634`, (D[ψ[r], r] /. r -> rStart1) == 
     epsilon, (D[ϕ[r], r] /. r -> rStart1) == 
     epsilon, (ϕ[rEnd1]) == Nphy1/rEnd1};
  sollst1 = 
   NDSolveValue[
    Flatten@{eqn2, bc2}, {ψ[r], ϕ[r]}, {r, rStart1, 
     rEnd1}](*Map[NDSolveValue[Flatten@{eqn2,bc2},{ψ[r],ϕ[
  r]},{r,rStart1,rEnd1},
  Method\[Rule]"BoundaryValues"\[Rule]{"Shooting",
  "StartingInitialConditions"\[Rule]{ψ[0]\[Equal]#}}]&,Range[-3,
  3,0.05]]*)
  ]
sollst1 = stableHunter3[-0.5074977775084505`]
Plot[Evaluate[sollst1 /. r -> r00], {r00, rStart1, rEnd1}]

enter image description here

Any thoughts how to proceed from there to solve $\Psi(r)$ at region where $r$ is small?

Attempt 5: Something weird happened. I start believing there is something related to the precision of NDSolve. Here is the code in which only the second element of the list gives me a nontrivial solution. 

stableHunter3[μ_] := Module[{},
  eqn2 = {-μ*ψ[r] - 
      1/2*(D[ψ[r], r, r] + 2/r*D[ψ[r], r]) + ϕ[
        r]*ψ[r] == 0, 
    D[ϕ[r], r, r] + 2/r*D[ϕ[r], r] == 4 Pi*ψ[r]^2};
  bc2 = {ψ[rEnd1] == 
     0.025963699315910634`, (D[ψ[r], r] /. r -> rStart1) == 
     epsilon, (D[ϕ[r], r] /. r -> rStart1) == 
     epsilon, (ϕ[rEnd1]) == -Nphy1/rEnd1};
  sollst1 = 
   NDSolveValue[
    Flatten@{eqn2, bc2}, {ψ[r], ϕ[r]}, {r, rStart1, 
     rEnd1}](*Map[NDSolveValue[Flatten@{eqn2,bc2},{ψ[r],ϕ[
  r]},{r,rStart1,rEnd1},
  Method\[Rule]"BoundaryValues"\[Rule]{"Shooting",
  "StartingInitialConditions"\[Rule]{ψ[0]\[Equal]#}}]&,Range[-3,
  3,0.05]]*)
  ]
sollst1 = 
 Map[stableHunter3[#] \
&,(*Range[-0.5074977775084505-0.02,-0.5074977775084505+0.02,
  0.005]*){-0.5074977775084505`, -0.5074977775084505`-0, \
-0.5074977775084505`+0.01}]
Plot[Evaluate[Flatten@sollst1 /. r -> r00], {r00, rStart1, rEnd1}, 
 PlotRange -> All]

enter image description here

Improve this question
$\endgroup$
13
  • $\begingroup$ Good question. Actually I think Psi's normalization is related to Phi through the second equation, although I don't think they have a simple relation at the boundary. $\endgroup$
    – Boson Bear
    Commented Jan 22, 2018 at 1:12
  • $\begingroup$ I meant $\nabla^2 \Phi = 4\pi \Psi^2$. Sorry is it true that we have two functions $\Psi(r), \Phi(r)$, and both are second order differentiated, so we have six `unknowns' in total to be determined. On the other hand, we have two differential equations, together with four boundary conditions so the system should be determined? $\endgroup$
    – Boson Bear
    Commented Jan 22, 2018 at 3:06
  • $\begingroup$ ah I think I understand what you mean now. So what I did was to trade $\Psi(\infty)=0$ with a boundary condition at zero, $\Psi(0) = ...$. For $\Phi(\infty)=0$ I haven't changed it to an initial condition for $\Phi(0)$ yet. Do you think that would help if I do that? $\endgroup$
    – Boson Bear
    Commented Jan 22, 2018 at 3:10
  • $\begingroup$ Yes m was assumed to be one. I just edited it out so we don't have to worry about m. Do you think I should try a larger rEnd1 say 10~100? It seems Mathematica starts struggling even at rEnd1=5. I'll try and get back. $\endgroup$
    – Boson Bear
    Commented Jan 22, 2018 at 12:50
  • 2
    $\begingroup$ I suggest that you compute symbolically the asymptotic solutions of the two equations, which decouple at large r. Then, use the asymptotic solutions as outer boundary conditions, because they tend to be much more accurate than just setting the functions to 0, especially at moderate rEnd1 == 5. Also, look at stackexchange answers with tag boundary-condition-at-infinity. But, remember that this problem is very hard, because you are, in effect, using the "Shooting" Method with three unknowns. $\endgroup$
    – bbgodfrey
    Commented Jan 22, 2018 at 13:04

1 Answer 1

Reset to default
6
$\begingroup$

It turns out that the `trial-and-error' method is the simplest, as suggested by a few people that I talked to. Below is a piece of sample code. 

epsilon = 0;
rEnd1 = 12;
rStart1 = 10^-5;
stableHunter3[μ_, A_, B_] := 
 Module[{}, 
  eqn2 = {-μ*ψ[r] - 
      1/2*(D[ψ[r], r, r] + 2/r*D[ψ[r], r]) + ϕ[
        r]*ψ[r] == 0, 
    D[ϕ[r], r, r] + 2/r*D[ϕ[r], r] == 4 Pi*ψ[r]^2};
  bc2 = {ψ[rStart1] == A, (D[ψ[r], r] /. r -> rStart1) == 
     epsilon, (D[ϕ[r], r] /. r -> rStart1) == 
     epsilon, (ϕ[rStart1]) == B};
  sollst1 = 
   NDSolveValue[
    Flatten@{eqn2, bc2}, {ψ, ϕ}, {r, rStart1, rEnd1}, 
    Method -> "StiffnessSwitching"]]

sollst = Table[
    sol = stableHunter3[muloop, psiloop, -4][[1]];
    solder[rr_] := D[sol[r], r] /. r -> rr;
    pt = rStart1;
    While[solder[pt] <= 0 && sol[pt] >= 0 && pt < rEnd1,
     pt += (rEnd1 - rStart1)/10;
     ];
    If[pt == rEnd1, {muloop, psiloop, sol}]
    , {muloop, -4, -3, 0.002}, {psiloop, 0.1, 0.4, 
     0.002}]; // AbsoluteTiming

Out[5]= {69.0395, Null}

(*I used twelve remote kernels to run it. 
Running locally it'll take about 15 min. 
You can try a coarse-grained lattice for 
the muloop and psiloop, and decrease rEnd1 
at the same time.*)

In[6]:= Flatten[DeleteCases[sollst /. Null -> Sequence[], {}], 1]

(-3.77  0.1 InterpolatingFunction[Domain: (0.00001  12.)
Output: scalar]
-3.756  0.106   InterpolatingFunction[Domain: (0.00001  12.)
Output: scalar]
-3.742  0.112   InterpolatingFunction[Domain: (0.00001  12.)
Output: scalar]
-3.65   0.152   InterpolatingFunction[Domain: (0.00001  12.)
Output: scalar])

sol = stableHunter3[-3.65, 0.152, -4][[1]];

solder[rr_] := D[sol[r], r] /. r -> rr;
Plot[Evaluate[{sol[r], solder[r]} /. r -> rrr], {rrr, rStart1, rEnd1},
  PlotRange -> All]

enter image description here

Helpful communications with @bbgodfrey, M. Hertzberg, E. Schiappacasse, S. Schwarz, L. Titus, @xzczd are appreciated.

Improve this answer
$\endgroup$
0

Your Answer

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