我试图求解修改后的 TOV 方程 dP'(r)/dr = -1.474 ε'(r)M'(r)/r^2 (1+P'(r)/ε'(r))(1 +11.210^-6 r^3 P'(r)/M'(r))[(1-2.948M'(r)/r)]^-1 dM'(r)/dr = bo * ε'(r) 其中 M(r)=M'(r)Mo, ε(r)=ε'(r)εο , εο= 1MeVfm^-3, GMo/c ^2=1.474 公里和 4pi/(Moc^2)=0.710^-40 s^2/(kgKm^2),bo=8.910^-7 km^-3,状态方程ε(P)=15.55P^0.666+76.71P^0.247。这些修改后的方程肯定是正确的,因为它们取自两篇不同的论文。我遇到的问题是半径不会随着 Pc 的不同值而改变,而 M 会改变一点点。我使用的Pc范围与论文中相同:Pc=1-1200 Mev*fm^-3。有人可以帮助我吗?这是我使用的代码:
import numpy as np
from scipy.integrate import solve_ivp
# Constants for the modified TOV equations
epsilon_0 = 1 # MeV/fm^3
b0 = 8.9e-7 # km^-3
# Equation of state: ε(P) = 15.55 * P^0.666 + 76.71 * P^0.247
def energy_density(P):
return 15.55 * P**0.666 + 76.71 * P**0.247
# Modified TOV equations: dP/dr and dM/dr
def tov_equations(r, y):
P, M = y
epsilon_r = energy_density(P) * epsilon_0 # Energy density from the equation of state
if P <= 0:
return [0, 0] # Stop when pressure becomes zero
# Modified TOV equations
dP_dr = -1.474 * (epsilon_r * M * (1 + P / epsilon_r)) / r**2 * \
(1 + 11.2e-6 * r**3 * P / M) * (1 - 2.948 * M / r)**-1
dM_dr = b0 * epsilon_r
return [dP_dr, dM_dr]
# Boundary conditions at r=0: central pressure and mass
def tov_initial_conditions(Pc):
M_c = 1e-2 # Small initial mass at the center
return [Pc, M_c]
# Solve TOV equations
def solve_tov(Pc):
# Initial conditions: central pressure and mass
y0 = tov_initial_conditions(Pc)
# Define the radial range for integration (we'll stop when pressure drops to ~0)
r_max = 50 # Maximum radius to integrate out to [km]
r_eval = np.linspace(1e-4, r_max, 10000) # Evaluation points in radius
# Solve TOV equations using solve_ivp (Runge-Kutta method)
solution = solve_ivp(tov_equations, [1e-4, r_max], y0, method='RK45', t_eval=r_eval, rtol=1e-6, atol=1e-6)
# Extract radius, pressure, and mass from the solution
radius = solution.t
pressure = solution.y[0]
mass = solution.y[1]
# Find the radius where pressure becomes negligible (surface of the star)
surface_indices = np.where(pressure <= 1e-10)[0]
if len(surface_indices) == 0:
# If no valid surface is found, use the last point
surface_index = -1
else:
surface_index = surface_indices[0]
R_star = radius[surface_index] # Radius at the surface
M_star = mass[surface_index] # Mass at the surface
# Handle cases where R_star or M_star may still be unrealistic (very small or zero)
if R_star <= 0 or M_star <= 0:
return None, None
return R_star, M_star
# Generate central pressure values from 1 to 1200 MeV/fm^3
Pc_values = np.linspace(1, 1200, 400)
# Arrays to store the results
R_values = []
M_values = []
# Loop over each central pressure value, solve the TOV equations, and store the results
for Pc in Pc_values:
R_star, M_star = solve_tov(Pc)
if R_star is not None and M_star is not None:
R_values.append(R_star)
M_values.append(M_star)
# Print a sample of the results
print("Sample results (Radius in km and Mass in dimensionless units):")
for i in range(0, len(R_values), 50): # Print every 50th result for sampling
print(f"Central pressure: {Pc_values[i]:.2f} MeV/fm^3 -> Radius: {R_values[i]:.2f} km, Mass: {M_values[i]:.2f}")
Sample results (Radius in km and Mass in dimensionless units):
Central pressure: 1.00 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.31
Central pressure: 151.25 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.26
Central pressure: 301.50 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.37
Central pressure: 451.75 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.30
Central pressure: 602.00 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.27
Central pressure: 752.25 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.25
Central pressure: 902.50 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.23
Central pressure: 1052.75 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.22
结果总是相同的,因为半径越大,压力就会越来越大。
因此压力永远不会低于
1e-10
; surface_indices
始终为空; surface_index
始终是-1
;并且半径始终为 0.08
。
由于压力显然应该在中心最高并在较高半径处减小,因此您的方程是错误/编码错误的。
我会首先查看
dP_dr
的标志。