我正在尝试使用 python 求解恒星建模的基本方程。我有 ODE 系统和所需的变量：
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()
边界条件是恒星的中心和表面。在中心（Mr=0），r=0，l=0。在表面（Mr=M，所有质量），rho=0，T = (L/(8pi(r^2)sigma)^(1/4)，其中 sigma 是玻尔兹曼常数。
当我运行代码时，目前显示
sol = solve_bvp(coupled_differential_equations, boundary_conditions, x, y_a)我需要有一个初步的猜测（比如拍摄方法）吗？我现在有点不知道该怎么办...
我对编码还很陌生，但如果有人可以提供一些关于如何让solve_bvp工作的指示，或者说实话我做错了什么，我将非常感激！
仍然有很多遗漏，所以这个答案将是不完整的，但我会解决你的一些问题（主要是语法）并告诉你你遗漏了什么。
首先，对于
scipy.integrate.solve_bvpf(x, y)yxdef 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])
（我还用
MMrMM其次，边界条件函数应该具有函数签名
bc(ya, yb)yaybPTrLdef 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)
这里的问题是（1）你缺少一个边界条件（我不确定你如何从其他值中得到
rho一旦有了这些，您就需要修复调用函数的方式。由于
MrrxMrMr_initial = ??
Mr_final = ??
Mr = np.linspace(Mr_initial, Mr_final, 100)
剩下的几行是：
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()
注意：我更改了绘图标签的顺序以匹配
coupled_differential_equations