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}
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}]
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]
r. Then, use the asymptotic solutions as outer boundary conditions, because they tend to be much more accurate than just setting the functions to0, especially at moderaterEnd1 == 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"Methodwith three unknowns. $\endgroup$