Boundary value problem with singularity and boundary condition at infinity

Question

I'm trying to solve the following boundary value problem on $[0,\infty]$:

$$f^{\prime \prime}=-\frac{1}{r} f^{\prime}+\frac{1}{r^{2}} f+m^{2} f+2 \lambda f^{3}$$

$$f(0)=0 \ ; f(\infty)=\sqrt{-m^2/(2\lambda)} $$

for some constants $m^2<0, \lambda>0$. There is no closed form but we should have $f$ monotonically increasing from $0$ to $\sqrt{-m^2/(2\lambda)}$. There is a removable singularity at $r=0$. The problem is just Bessel's equation plus a term in $f^3$.

I'm trying to solve this with Scipy's integrate.solve_bvp which can solve multi-boundary problems with a singularity at one boundary, defining $y=\begin{bmatrix} f \\ rf' \\ \end{bmatrix}$ so that

$$y'=\begin{bmatrix} 0 \\ r(m^2f+2\lambda f^3) \\ \end{bmatrix}+\frac{1}{r}\begin{bmatrix} 0 \ \ 1 \\ 1 \ \ 0 \\ \end{bmatrix}y$$

I impose the boundary condition at infinity at some large value max_x. Unfortunately my code, following the structure of the example at https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html, gives the wrong solution:

import scipy.integrate
import numpy as np
import matplotlib.pyplot as plt

m_squared=-1
Lambda=1
asymptote=np.sqrt(-m_squared/(2*Lambda))

#evaluate infinity b.c here
max_x=100

def fun(r,v):
    z=(m_squared*v[0]+(2*Lambda)*(v[0]**3) )*r
    return np.vstack((z-z, z))
  
#boundary condition
def bc(ya,yb):
    return np.array([ya[0], yb[0]-asymptote])

# to treat singularity
S=np.array([[0,1],[1,0]])


x=np.linspace(0,max_x,5000)

# guess for vector y at points x  
y=np.zeros((2, len(x)))
y[0,-1]=asymptote
print(y)

#solve
res=scipy.integrate.solve_bvp(fun, bc, x, y, p=None, S=S)
x_plot=np.linspace(0,max_x,1000)
y_plot=res.sol(x_plot)[0]

plt.plot(x_plot,y_plot,label="numerical")
plt.axhline(asymptote,linestyle="--",label="asymptote")

plt.xlabel("r")
plt.ylabel("f")
plt.legend()

I checked that modifying the above code to solve e.g $f''=f-1$ with $f(0)=0, f(\infty)=1$ works fine. There are no singularity in this case, so it suffices to modify fun and set S=None.

Is there an issue with my code or should I use a different boundary value solver?

Answer 1

The left boundary is really singular, and the desired right boundary behavior is unstable. \begin{align} f''(r)+\frac{f'(r)}{r}-\frac{f(r)}{r^2}&=-f(r)+2\lambda f(r)^3 \\ r^2f''(r)+rf'(r)-f(r)&=r^2[-f(r)+2\lambda f(r)^3] \end{align} The second form highlights the regular nature of the singularity at $0$. The indicial equation has roots $\pm1$. For a finite solution the basis solution based on $r^{-1}$ has to be avoided. This means that the solution is mostly linear, $f(r)=c_1r$, close to $r=0$. This is equivalent to the equation $f(r)-rf'(r)=0$. The next correction term would follow from $$ (2^2-1)c_2=0\\ (3^2-1)c_3=-c_1\\ (4^2-1)c_4=-c_2\\ (5^2-1)c_5=-c_3+2λc_1^3\\ ... $$ to be the cubic term $-\frac18c_1r^3$ etc. Thus set the left interval boundary at $a=\sqrt{tol}$ where $tol$ is the relative error tolerance. Then the shift of the left boundary from $0$ to $a$ with condition $af'(a)-f(a)=0$ will give a sufficient accuracy.

In the far field of large $r$ the first equation reduces to $$ f''(r)=-f(r)+2λf(r)^3. $$ The only asymptotically converging solutions are those that reach the saddle points $f_\infty=\pm\sqrt{1/(2λ)}$ via their heteroclinic orbits or stable manifolds. The approximate equation has a first integral $$ f'(r)^2+f(r)^2-λf(r)^4=C $$ For $(f,f')=(f_\infty,0)$ this gives $C=\frac1{4λ}$ which completes the square in $$ f'(r)^2=\frac1{4λ}(2λf(r)^2-1)^2. $$ We are looking for the solution that grows to $f_\infty$ from zero, thus $f'(r)>0$, $2λf(r)^2-1<0$, $$ f'(r)=\frac1{2\sqrtλ}(1-2λf(r)^2), \qquad\left\{ \small\begin{aligned} \sqrt{λ}f'(r)+λf(r)^2&=\frac12, \\ (e^{\sqrt{λ}F(r)})''&=\frac12 e^{\sqrt{λ}F(r)},\\ \sqrt{λ}F(r)&=\ln(Ae^{t/\sqrt2}+Be^{-t/\sqrt2}),\\ f(r)&=\frac1{\sqrt{2λ}}\frac{Ae^{t/\sqrt2}-Be^{-t/\sqrt2}}{Ae^{t/\sqrt2}+Be^{-t/\sqrt2}} \end{aligned}\right. $$ This left identity is a condition that can be imposed at $r=b$ to get a solution that gets close to the saddle point also in the original equation. The far-field solution on the right specializes here to $f(r)=f_\infty\tanh((t-c)\sqrt2)$ where $c$ is determined by the value of $f(b)$.

λ=0.2

#evaluate infinity b.c here
max_x=20

def fun(r,y):
    z=y[0]*(2*λ*y[0]**2-1)
    return y[1]/r, y[0]/r+r*z
  
#boundary condition
def bc(ya,yb):
    return ya[0]-ya[1], 2*λ**0.5*yb[1]+max_x*(2*λ*yb[0]**2-1)

x=np.logspace(-4,0,10)*max_x

# guess for vector y at points x  
y=[np.tanh(x),x/np.cosh(x)**2]

#solve
res=solve_bvp(fun, bc, x, y, tol=1e-4, max_nodes=2500)
print(res.message,"  ", len(res.x))

The initial guess could be brought closer to the far-field solution,... but this works rapidly enough.

Using the singular term mechanism of solve_bvp would be more compact, but there is some problem with the error estimation that forces infinite grid refinement at the singular point (until max nodes exceeded). This run with tol=1e-4 requires 35 nodes, with tol=1e-6 this increases to 147 nodes, far below the default boundary.