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''=-f'/r+f/r^2+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. This 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=[f,rf'] so that

y'=[0,r(m^2f+2\lambda f^3)]+(1/r)*[[0,1],[1,0]]*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

I read the equation from your source as

f''(r)+f'(r)/r-f(r)/r^2 = 2*lambda*eta^2*f(r)*(f(r)^2-1)

with the proposed parameters lambda=0.2, eta=2.

In the long-term limit, the left side reduces to the second derivative and the equation to a conservative system with a center at f=0 and two saddle points at f=+-1. The task is to find a solution curve that converges to a saddle point. In more practical terms, this is similar to the task to push a rigid pendulum in such a way that it ends up in the upright position, or moves ever closer to that position.

Writing f=1-g(r) for a solution approaching the saddle point at f=1, the equation is approximately

g''(r) = a^2*g(r),  a^2=4*lambda*eta^2=3.2

This again characterizes this equilibrium as a saddle point, the solutions converging toward it satisfy the reduced ODE g'(r)=-a*g(r). This can be used as upper boundary condition. Translated into the state vector this gives

def bc(ya,yb):
    return np.array([ya[0], yb[1]+a*max_x*(yb[0]-1)])

(replace the equilibrium constant 1 with asymptote if you want to stay with your version).

I got good results from that, with the parameters in the paper as well as with your parameters.

However, the singular mode of the solver seems broken, it inserts nodes up to the allowed max_nodes close to zero where the solution should be simply linear. I set the initial guess to

x=np.logspace(-5,np.log10(max_x),10)
x[0]=0

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

so that this maximal node number is not violated from the start.

Not using the singular mechanism, transferring the singular terms back into the ODE function and using that the finite solutions are almost linear in the initial segment, one can use f(r)=r*f'(r) as initial condition, ya[0]-ya[1] == 0. The interval start is then some small positive number. This results in reasonable node numbers in the solution, 25 for the default tolerance 1e-3 and 100 for tolerance 1e-6.