I am trying to solve a boundary value problem with Python. I have been using scipy.integrate.solve_bvp but the result that it is giving me is completely wrong. Basically my code is as follows:

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

r_in = 0.0
r_end = 500.
l = 0.
rho_c = 1.
m = 1.

def system_DE(r,y,p_self,l=l,rho_c=rho_c):    
     R = y[0]
     ph = y[1]
     a = y[2]
     al = y[3]

     om = p_self[0]
     u = p_self[1]

    dR_dr = ph
    da_dr = a*(((2.*l+1.)*r/2.)*(-om**2.*a**2.*R**2./al**2.+ph**2.+l*(l+1.)*a**2.*R**2./r**2.+m**2.*a**2.*R**2.)-(a**2.-1.)/(2.*r))
    dal_dr = al*(da_dr/a-l*(l+1.)*(2.*l+1.)*a**2.*R**2./r-(2.*l+1.)*m**2.*a**2.*r*R**2.+(a**2.-1.)/r)
    dph_dr = -2.*ph/r+(l*(l+1.)*(2.*l+1.)*a**2.*R**2./r+(2.*l+1.)*m**2*a**2*r*R**2-(a**2-1)/r)*ph-om**2.*a**2.*R/al**2.+l*(l+1.)*a**2.*R/r**2.+m**2.*a**2.*R

     return [dR_dr,dph_dr,da_dr,dal_dr]

def bc(ya,yb,p_sch,l=l,rho_c=rho_c,r_in=r_in,r_end = r_end,m=m):

    om = p_self[0]
    u = p_self[1]

    if l==0.:
        return np.array([ya[0]-rho_c,ya[1]-rho_c*r_in,yb[0]*((m**2.-om**2.)**(1./2.)+1./(r_end))+yb[1],ya[2]-1.,ya[3]-u,yb[3]*yb[2]-1.])
        return np.array([rho_c*r**l,rho_c*l*r**(l-1.),1.,u])

n_r = 1000
r = np.linspace(r_in,r_end,n_r)
y = np.zeros((4,r.size))
y[0,0] = rho_c
y[0,1] = rho_c

y[2,:] = 1.

p_self=[0.8,0.85] #[om,u]

sol = solve_bvp(system_DE,bc, r, y, p=p_self, tol=0.01,max_nodes=100000)

x_plot = np.linspace(r_in, 200  , 5000)
L = sol.sol(x_plot)[0]
Ph = sol.sol(x_plot)[1]
a = sol.sol(x_plot)[2]
alpha = sol.sol(x_plot)[3]

plt.plot(x_plot, L)
plt.plot(x_plot, a)
plt.plot(x_plot, alpha)

The idea is that with the parameters that I am giving to the code I should obtain something similar to the figure [1] (up) on this article: https://arxiv.org/pdf/0908.2435.pdf, however if I run the code I am obtaining something completely wrong. Is there another solver that I could use? Is there something wrong with my code?

Thanks for your time. Regards,

Luis P.

  • $\begingroup$ Check your code, I am using Spyder and I am seeing some nan values for Ph, a and alpha. Are you sure your implementation is correct? $\endgroup$ – David Oct 14 '18 at 2:34
  • $\begingroup$ Why do you use the global p_self in bc when the passed parameter vector is p_sch? You should consult the fields "sol.success" and "sol.message", there might be something about a "singular Jacobian". $\endgroup$ – LutzL Nov 9 '18 at 9:45

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.