Avoiding divergent solutions with odeint? shooting method

I am trying to solve an equation in Python. Basically what I want to do is to solve the equation:

$$\frac{1}{x^2}\frac{d}{dx}\left(Gam \frac{dL}{dx}\right)+L\left(\frac{a^2x^2}{Gam}-m^2\right)=0$$

This is the Klein-Gordon equation for a massive scalar field in a Schwarzschild spacetime. We know $m$ and $Gam=x^2-2x$.

The initial/boundary condition that I know are $L\rvert_{2+\epsilon}=1$ and $L\rvert_{\infty}=0$.

Notice that the asymptotic behavior of the equation is

$$L(x\to\infty)\to \frac{e^{\pm(m^2-a^2)x}}{x}$$

Then, if $a^2>m^2$ we will have oscillatory solutions, while if $a^2 < m^2$ we will have a divergent and a decay solution. What I am interested in is the decay solution; however when I am trying to solve the above equation transforming it as a system of first-order differential equations and using the shooting method in order to find the $a$ that can give me the behavior that I am interested about, I am always having a divergent solution (for $0<a^2<m^2$). I suppose that it is happening because odeint is always finding the divergent asymptotic solution.

Is there a way to avoid or tell to odeint that I am interested in the decay solution? If not, do you know a way that I could solve this problem? Maybe using another method for solving my system of differential equations? If yes, which method?

Basically what I am doing is to add a new system of equations for $a$

$$\frac{d^2a}{dx^2}=0,\\ \frac{da}{dx}(2+\epsilon)=0,\\ a(2+\epsilon)=a_0)$$

in order to have $a$ as a constant. Then I am considering different values for $a_0$ and asking if my boundary conditions are fulfilled.

EDIT 1.

Hello, let me try to explain a little more my issue

I am incorporating the value at infinity considering the assimptotic behavior, it means that I will have a relation between the field and its derivative. I will post the code for you if it is helpful:

from IPython import get_ipython
get_ipython().magic('reset -sf')
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from math import *
from scipy.integrate import ode


These are initial conditions for Schwarzschild. The field is invariant under reescaling, then I can use $L(2+\epsilon)=1$

def init_sch(u_sch):
om = u_sch[0]
return np.array([1,0,om,0]) #conditions near the horizon, [L_c,dL/dx,a,da/dx]


These are our system of equations

def F_sch(IC,r,rho_c,m,lam,l,j=0,mu=0):
L = IC[0]
ph = IC[1]
om = IC[2]
b = IC[3]

Gam_sch=r**2.-2.*r

dR_dr = ph
dph_dr = (1./Gam_sch)*(2.*(1.-r)*ph+L*(l*(l+1.))-om**2.*r**4.*L/Gam_sch+(m**2.+lam*L**2.)*r**2.*L)
dom_dr = b
db_dr = 0.
return [dR_dr,dph_dr,dom_dr,db_dr]


Then I try for different values of "om" and ask if my boundary conditions are fulfilled. p_sch are the parameters of my model. In general what I want to do is a little more complicated and in general I will need more parameters that in the just massive case. Howeve I need to start with the easiest which is what I am asking here

p_sch = (1,1,0,0) #[rho_c,m,lam,l], lam and l are for a more complicated case
ep = 0.2
ep_r = 0.01
r_end = 500
n_r = 500000
n_omega = 1000
omega = np.linspace(p_sch[1]-ep,p_sch[1],n_omega)
r = np.linspace(2+ep_r,r_end,n_r)
tol = 0.01
a = 0

for j in range(len(omega)):
print('trying with $omega =$',omega[j])
omeg = [omega[j]]
ini = init_sch(omeg)
Y = odeint(F_sch,ini,r,p_sch,mxstep=50000000)
print Y[-1,0]
#This is my boundary condition that I am implementing. Basically this should be my conditions at "infinity
if abs(Y[-1,0]*((p_sch[1]**2.-Y[-1,2]**2.)**(1/2.)+1./(r[-1]))+Y[-1,1]) < tol:
print(j,'times iterations in omega')
print("R'(inf)) = ", Y[-1,0])
print("\omega",omega[j])
omega_1 = [omega[j]]
a = 10
break
if a > 1:
break


Basically what I want to do here is to solve the system of equations giving different initial conditions and find a value for "a=" (or "om" in the code) that should be near to my boundary conditions. I need this because after this I can give such initial guest to a secant method and try to fiend a best value for "a". However, always that I am running this code I am having divergent solutions that it is, of course, a behavior that I am not interested. I am trying the same but considering the scipy.integrate.solve_vbp, but when I run the following code:

from IPython import get_ipython
get_ipython().magic('reset -sf')
import numpy as np
import matplotlib.pyplot as plt
from math import *
from scipy.integrate import solve_bvp

def bc(ya,yb,p_sch):
m = p_sch[1]
om = p_sch[4]
tol_s = p_sch[5]
r_end = p_sch[6]

return np.array([ya[0]-1,yb[0]-tol_s,ya[1],yb[1]+((m**2-yb[2]**2)**(1/2)+1/r_end)*yb[0],ya[2]-om,yb[2]-om,ya[3],yb[3]])

def fun(r,y,p_sch):
rho_c = p_sch[0]
m = p_sch[1]
lam = p_sch[2]
l = p_sch[3]

L = y[0]
ph = y[1]
om = y[2]
b = y[3]

Gam_sch=r**2.-2.*r

dR_dr = ph
dph_dr = (1./Gam_sch)*(2.*(1.-r)*ph+L*(l*(l+1.))-om**2.*r**4.*L/Gam_sch+(m**2.+lam*L**2.)*r**2.*L)
dom_dr = b
db_dr = 0.*y[3]
return np.vstack((dR_dr,dph_dr,dom_dr,db_dr))

eps_r=0.01
r_end = 500
n_r = 50000
r = np.linspace(2+eps_r,r_end,n_r)
y = np.zeros((4,r.size))
y[0]=1
tol_s = 0.0001
p_sch= (1,1,0,0,0.8,tol_s,r_end)

sol = solve_bvp(fun,bc, r, y, p_sch)


I am obtaining this error: ValueError: bc return is expected to have shape (11,), but actually has (8,). ValueError: bc return is expected to have shape (11,), but actually has (8,).

• can you verify the equations as I made them better readable using MathJax? – Anton Menshov Jun 16 '18 at 9:43
• This is a general issue with using a shooting method for a BVP with one growing and one decaying: any tiny change (such as due to roundoff error) to $a$, $a\to a+\epsilon$ produces a change to the solution of the ODE of the form $+\delta e^{+|m^2-a^2|x}/x$, which causes $L|_\infty=0$ to never be satisfied even approximately. In this case even though the BVP might be stable, the corresponding IVP is unstable. There are multiple standard ways of working around this, have you tried just using a collocation method like solve_bvp? – Kirill Jun 16 '18 at 11:49
• Hello, I didn't try with solve_bvp because I am not sure how can I tell to solve_bvp that I am only interested in $0<a^2<m^2$. Because it results that when $a^2>m^2$ I will have oscillatory solutions of the form $sin(kx)/x$ which for $x$ big enough I will have the behavior $L|_{\infty}=0$, however I am not interested in that solution, I am interested only when $a^2<m^2$. Is there a way to tell to solve_bvp this? – Luis Enrique Padilla Albores Jun 16 '18 at 19:35
• Yes, the equations are ok. Thanks for the changes. – Luis Enrique Padilla Albores Jun 16 '18 at 19:40
• I added more explanation of what I am doing, could you check it, please? – Luis Enrique Padilla Albores Jun 17 '18 at 4:43