# Boundary value problem with singularity and boundary condition at infinity

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?

• Are you sure this is a simple BVP rather than an eigenvalue problem? Does this come from the separation of a wave equation in cylindrical coordinates? Mar 24, 2022 at 14:20
• @Endulum this is definitely a boundary value problem. It does, although there is no time-dependence. See e.g equation 11 and related plots here arxiv.org/pdf/2203.00894.pdf Mar 24, 2022 at 19:06
• @Endulum sorry Eq. 23 and plot here is clearer: researchgate.net/publication/… Mar 24, 2022 at 19:18
• Cross-posted in stackoverflow.com/questions/71600263/…, the problem is similar to push a ball up a hill with the aim of it stopping at the top. Just having it at the top at some time will not have it with zero velocity there. Forcing it to be close to a heteroclinic curve at the upper boundary will give correctly looking solutions. Mar 25, 2022 at 11:43

## 1 Answer

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.

• I didn't realize the solver singularity treatment is broken! Thanks for pointing that out. Additionally there was the issue with the boundary condition at infinity, setting the function to its asymptote isnt sufficient for stability, I guess. Mar 25, 2022 at 21:30
• Not "broken" broken, it arrives at some result that looks correct. But it terminates with an error condition, which is not the desired behavior. And the refinement cycles slow the solver down without providing substantial improvement. Mar 25, 2022 at 21:45