Solving a system of coupled bvp ode equations using scipy's solve_bvp for stellar modeling

Question

I'm trying to solve the basic equations for stellar modeling using python. I have the system of ODEs and the variables needed:

from scipy.integrate import odeint
from scipy.integrate import solve_bvp
import numpy as np
import matplotlib.pyplot as plt

#Constants (for now)

G = 6.67e-8
X = 0.71
Y = 0.28
Z = 0.01
#Mr is the mass contained within a certain radius
#kappa is opacity
#nabla is the energy transport assuming the whole star is radiative
#rho is density
#r, P, L, T is max radius, pressure, luminosity, and temperature respectively 
def coupled_differential_equations(P, T, r, L, Mr):

  #rho = (P*(2.1945*1e-3)/(0.0821*T))

  #kappa = (4e25*(1+X)*(Z+1e-3)*(P*(2.1945*1e-3)/(0.0821*T)/(T**3.5)))

  #nabla = ((3 * 4e25*(1+X)*(Z+1e-3)*(P*(2.1945*1e-3)/(0.0821*T)/(T**3.5)) * L * P)/(16 * np.pi * 1.38e-23 * T**4)*(1 / G * Mr))

  dr_dMr = (1 / (4 * np.pi * r**2 * P*(2.1945*1e-3)/(0.0821*T)))
  dP_dMr = -((G * Mr) / (4 * np.pi * r**4))
  dL_dMr = 9.5 * 10e-37 * X**2 * (P*(2.1945*1e-3)/(0.0821*T))**2 * T**4
  dT_dMr = -(((3 * 4e25*(1+X)*(Z+1e-3)*(P*(2.1945*1e-3)/(0.0821*T)/(T**3.5)) * L * P)/(16 * np.pi * 1.38e-23 * T**4)*(1 / G * Mr))*(T/P))*((G * M)/(4 * np.pi * r**4))
  return np.vstack([dr_dMr, dP_dMr, dL_dMr, dT_dMr])

 #***************************************************************************

def boundary_conditions(ya, yb, P, T, L, Mr):
  # Boundary conditions at the center
  center_conditions = np.array([ya[0], ya[2], ya[3]])
    # Boundary conditions at the surface
  surface_conditions = np.array([yb[1], yb[3] - 1e-6, yb[2] - Mr])
  return np.concatenate([center_conditions, surface_conditions])
#Current issue is that P, T, r, L, Mr, kappa is not defined outside of the function??? seems to require initial guess but it's a bvp so it shouldn't need one???

# Initial values and range for r
initial_r = 1e-3
final_r = 3.9e5

x = np.linspace(initial_r, final_r, 100)
y_a = np.zeros((4, x.size))
# Solve BVP using solve_bvp 
sol = solve_bvp(coupled_differential_equations, boundary_conditions, x, y_a)


# Extract the solution
r_solution = sol.x
y_solution = sol.y

# Plot the results
plt.plot(r_solution, y_solution[0], label='r')
plt.plot(r_solution, y_solution[1], label='P')
plt.plot(r_solution, y_solution[2], label='L')
plt.plot(r_solution, y_solution[3], label='T')
plt.xlabel('r')
plt.legend()
plt.show()

The boundary conditions are the center of the star and the surface. At the center (Mr=0), r=0, and l=0. At the surface (Mr=M, all the mass), rho=0, T = (L/(8pi(r^2)sigma)^(1/4), where sigma is the Boltzman constant.

When I run the code, it is currently saying that the coupled_differential_equations parameter in the sol = solve_bvp(coupled_differential_equations, boundary_conditions, x, y_a) is missing arguments, but when I added in (Mr, P, T, L, r, kappa), it says that the varibles are not defined. The scipy handbook (https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html) seemed to not need arguments when plugging the system of odes into the solve_bvp.

Do I need to have an initial guess (something like the shooting method)? I'm a bit lost on what to do right now...

I'm still rather new to coding, but if anyone can give some pointers on how to get the solve_bvp to work, or honestly anything else I'm doing wrong, I'd appreciate it a lot!

Answer 1

There is still a lot missing, so this answer will be incomplete, but I will address some of your issues (mostly syntax) and tell you what you're missing.

First of all, for scipy.integrate.solve_bvp the ODE function signature should be f(x, y). You currently do not have that because you're passing all the y values independently and have x at the end. Here is what it should look like:

def coupled_differential_equations(x, y):
    Mr = x
    P, T, r, L = y    
    dr_dMr = (1 / (4 * np.pi * r**2 * P*(2.1945*1e-3)/(0.0821*T)))
    dP_dMr = -((G * Mr) / (4 * np.pi * r**4))
    dL_dMr = 9.5 * 10e-37 * X**2 * (P*(2.1945*1e-3)/(0.0821*T))**2 * T**4
    dT_dMr = -(((3 * 4e25*(1+X)*(Z+1e-3)*(P*(2.1945*1e-3)/(0.0821*T)/(T**3.5)) * L * P)/(16 * np.pi * 1.38e-23 * T**4)*(1 / G * Mr))*(T/P))*((G * Mr)/(4 * np.pi * r**4))
    return np.vstack([dr_dMr, dP_dMr, dL_dMr, dT_dMr])

(I also replaced the M with Mr since M is undefined, but maybe M should be defined somewhere.

Second of all, the boundary condition function should have the function signature bc(ya, yb), but again you're passing more. Here, ya and yb will contain the values of P, T, r, and L (in that order) at the initial and final mesh locations, respectively. You have a system of 4 first-order ODEs, so you can only specify 4 boundary conditions. That's a mathematical requirement (otherwise you overdefine the system). Here is what your boundary condition should look like:

def boundary_conditions(ya, yb):
    # At the center (Mr=0), r=0, and l=0
    # At the surface (Mr=M, all the mass), rho=0, T = (L/(8pi(r^2)sigma)^(1/4),
    # where sigma is the Boltzman constant.
    center_conditions = [ya[2], ya[3]]
    surface_conditions = [yb[1] - (L/(8*np.pi*yb[2]**2*sigma))**0.25, ??]
    return np.array(center_conditions + surface_conditions)

The issue here is (1) you're missing one of the boundary conditions (I'm not sure how you get rho from the other values) and (2) you need to define sigma somewhere (that's simple).

Once you have these, you need to fix how you're calling the function. Since Mr is your independent variable and not r (I'm not sure why you've done it this way, but you have), your mesh x must be the range of values for Mr instead. I don't know what those are, so you will have to fill them in.

Mr_initial = ??
Mr_final = ??
Mr = np.linspace(Mr_initial, Mr_final, 100)

The remaining lines will be:

y0 = np.zeros((4, Mr.size))
sol = solve_bvp(coupled_differential_equations, boundary_conditions, Mr, y0)

y_solutions = sol.y

fig, axes = plt.subplots(1, 4)
labels = ["P", "T", "r", "L"]
for y_sol, label, ax in zip(y_solutions, plot_labels, axes):
    ax.plot(Mr, y_sol)
    ax.set_xlabel("Mr")
    ax.set_ylabel(label)
fig.show()

Note: I changed the order of the plot labels to match the ordering in coupled_differential_equations (the order doesn't matter as long as you're consistent), and I put the plots in their own subplots since they will probably have different scales.