I am trying to solve with high precision the following coupled system $(f,h)$ on $[0,\infty]$:

$$-h''-\frac{1}{r}h'+\lambda_c h(f^2-1)+4\lambda_h h^3=0$$ $$-f''-\frac{1}{r}f'+\frac{1}{r^2}f+ f(-\lambda_f+2\lambda_c h^2)+\lambda_f f^3=0$$ where $\lambda_h, \lambda_f, \lambda_c $ are positive constants, and with boundary conditions $f(0)=0, f(\infty)=1, h'(0)=0, h(\infty)=0$. The condition $h'(0)$ is necessary for $h$ to be finite at the origin.

The issue

The problem is tricky because of the singularity at $r=0$, and the boundary conditions at infinity which preclude a naive shooting approach. Further the mixed boundary conditions on $h$ preclude a naive relaxation method as $h(0)$ is not known.

In the following I'll fix $\lambda_h=\lambda_f=1$. I am having trouble getting scipy's solve_bvp to converge for this system for arbitrary $\lambda_c$. Solving $f$ by itself by removing the coupling term with $h$ (no "backreaction") works fine. When solving the system with $h$ when $\lambda_c \sim 1 $, the solution converges to 0.1% accuracy, and the results look reasonable.

However, I am having trouble pushing to higher precisions, and the solution does not converge if $\lambda_c$ is too large, and has incorrect behavior for small $\lambda_c$ (setting $\lambda_c=0$, we should have $h=0$, however I have not been able to get that with precision higher than 0.1%).


The following snippet produces the below plot for $\lambda_c=1$ which converges and looks reasonable.

    def solve_h_f_sys(λ_c=1, λ_h=0, λ_f = 1, back_reaction=True, tol=1e-7,nb_rvals=500000,max_rval=100,max_nodes=1000000):
        x = np.logspace(-7,0,nb_rvals)*max_rval
        min_rval = x[0]
        br = int(back_reaction)
        def fun(r,y): # y is (f, rf', h, rh')
            z_h = λ_c*( (y[0])**2 -1) + 4*λ_h*(y[2]**2)
            z_f = λ_f*(y[0]**2-1) + 2*λ_c*(y[2]**2)*br
            return y[1]/r, y[0]/r + r*y[0]*z_f, y[3]/r, r*y[2]*z_h
        def bc(ya,yb):
            return ya[0]-ya[1], ((2/λ_f)**0.5)*yb[1] + max_rval*(yb[0]**2-1), ya[3]/min_rval, yb[2]+yb[3] 
        y=[np.tanh(x), x/np.cosh(x)**2, 1/(1+x), -1/ (1+x)]
        res=scipy.integrate.solve_bvp(fun, bc, x, y, tol=tol, max_nodes=max_nodes)
        print(res.message,"  ", len(res.x))
        f = res.sol(x)[0]
        h = res.sol(x)[2]
        rms_residuals = res.rms_residuals

    return x, h, f, rms_residuals

x, h, f, rms_residuals = solve_h_f_sys(λ_c=1, λ_h=1, λ_f = 1, back_reaction=True, tol=1e-3,nb_rvals=500000,max_rval=200,max_nodes=1000000)


enter image description here

Setting $\lambda_c=1.2$ the solution does not converge and we get:

enter image description here

For small $\lambda_c=0$ I cannot get $h(0)=0$ to be true better than at the 0.1% level. Zoomed in plot for $\lambda_c=0$:

enter image description here

Asymptotic analysis

I find that at small $r$, $f(r) \sim \alpha r $ and $h(r) \sim a +br^2 +cr^4$ for constants $\alpha, a,b,c$. This seems to be respected by the numerical solution. At large $r$, if there is no coupling between $f$ and $h$ then we have $f(r)=1-C / r^2+\mathcal{O}\left(1 / r^4\right)$. Assuming $f$ has the same form when the coupling is turned on, I find the ansatz $h = K/r$ satisfies the first equation if $K^2 = (1+2C\lambda_c)/(4\lambda_h)$. Substituting in the second equation, this requires $C = (1+2\lambda_cK^2)/(2\lambda_f)$. There is a solution $(C,K)$ as long as $\lambda_c^2 < 2\lambda_h\lambda_f$. This suggests the ansatz $h \propto 1/r$ at large $r$ is maybe a poor approximation.

For now I the boundary conditions at infinity I am using are the ones suggested in this answer for $f$: $$f^{\prime}(r)^2=\frac{\lambda_{f}}{2}\left(1-f(r)^2\right)^2$$

These should be appropriate even when $f$ couples to $h$ because $h$ must go to zero at infinity so the backreaction term can be neglected. For $h$, I am imposing that $h \propto 1/r$. This is likely incorrect but better than imposing $h=0$ at infinity which does not give correct-looking solutions.

  • $\begingroup$ @LutzLehmann this is the post I was referring to. $\endgroup$
    – phyphy
    Oct 9, 2022 at 4:02
  • $\begingroup$ Lutz Lehmann won't see this: I don't think at-notifications work for a user that did not participate in the answer or comment thread before. You can ping him in the chat instead though. $\endgroup$ Oct 9, 2022 at 8:53
  • $\begingroup$ It appears, in review of the cited answer, that the far-field approximation is not precise enough. Meaning that the omitted terms are possibly not that small for medium large $r_\max$. That should also be the reason that the variant using the S matrix mechanism of solve_bvp for the singularity does not converge, the derivative gets a visual hook up at the end that is most likely forced by the boundary condition. Some correction terms need to be added. $\endgroup$ Oct 9, 2022 at 9:55
  • $\begingroup$ @LutzLehmann perhaps there is some alternate solution method which allows one to be more agnostic about the boundary conditions? $\endgroup$
    – phyphy
    Oct 9, 2022 at 10:04


Your Answer

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