# Numerical precision on tricky coupled nonlinear boundary value problem on infinite interval

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%).

code

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)]

#solve
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)

plt.plot(x,f,label="f")
plt.plot(x,h,label="h")
plt.xlabel("r")
plt.legend()


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

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$$:

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.

• @LutzLehmann this is the post I was referring to. Oct 9, 2022 at 4:02
• 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. Oct 9, 2022 at 8:53
• 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. Oct 9, 2022 at 9:55
• @LutzLehmann perhaps there is some alternate solution method which allows one to be more agnostic about the boundary conditions? Oct 9, 2022 at 10:04