Python中子星状态方程代码无输出问题排查求助

问题根源分析 #

你的代码出现空文件或写入异常，主要由以下几个问题导致：

• 文件覆盖与写入冲突：同时使用open()手动写入和np.savetxt()操作同一文件，且np.savetxt默认以写入模式打开文件，每次调用都会清空原有内容；手动打开的文件句柄与np.savetxt的文件指针不共享，导致写入混乱。
• 循环逻辑错误：for循环嵌套while循环的逻辑存在问题，当i接近数组长度N时，i+1会超出数组索引范围，触发IndexError导致程序提前终止，无法执行f.close()完成文件写入；且while循环内未更新i，会陷入无限循环（或触发索引错误）。
• 文件句柄未及时刷新：手动写入文件后未调用flush()，缓冲区内容可能未写入磁盘。

修复后的代码 #

"""
Created on Wed Oct 12 20:44:10 2022
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

### Constants
h_bar = 1.055e-27      # reduced planck constant (kg m2 s-1)
c = 2.99e+10           # m/s
G = 6.67430e-8        # Gravitational constant m^3/(kg*s)
pi = 3.1415926535897
gamma = 5/3
m_sun = 2.998e33
K = 1e10
rho_c = 2e15           # Density central kg/m^3
P_c = K * rho_c ** gamma
km = 1e5
dr = 10e3
M0 = (4 * pi * (dr**3) * rho_c) / 3

#### Initial values for r = 0
N = 1600
m = np.zeros(N)
P = np.zeros(N)
rho = np.zeros(N)
r = np.zeros(N)
GR = np.zeros(N)

r[0] = 0
m[0] = 0
P[0] = P_c
rho[0] = rho_c
GR[0] = 0

####### 使用单一方式写入文件，避免冲突
# 初始化文件，写入表头或初始数据
with open("radial.out", "w") as f:
    # 先写入列名（可选，方便后续读取）
    f.write("r,rho,P,m,GR\n")
    # 写入初始行数据
    initial_line = f"{r[0]},{rho[0]},{P[0]},{m[0]},{GR[0]}\n"
    f.write(initial_line)

# 处理r[1]的情况
r[1] = dr
rho[1] = rho_c
P[1] = P_c
m[1] = 4 * pi * r[1]**3 * rho[1] / 3 
GR[1] = 0

# 写入r[1]的数据
with open("radial.out", "a") as f:
    line = f"{r[1]},{rho[1]},{P[1]},{m[1]},{GR[1]}\n"
    f.write(line)

####### 修正循环逻辑，避免索引越界和无限循环
for i in range(1, N-1):  # 只循环到N-2，避免i+1超出数组长度
    if P[i] <= 0:
        break  # 压力小于等于0时停止计算
    
    # 计算下一个步长的参数
    m[i+1] = m[i] + 4 * pi * r[i]**2 * rho[i] * dr
    rho[i+1] = (P[i]/K) ** (3/5)
    # 避免分母为0：当m[i]或r[i]为0时跳过GR计算（但这里i从1开始，r[i]不为0）
    denominator = 1 - (2 * G * m[i]) / (r[i] * c**2)
    if denominator <= 0:
        # 避免GR计算出现除以0或负数开方，直接停止
        break
    GR[i+1] = (1 + (4 * pi * r[i]**3 * P[i])/(m[i] * c**2)) * (1 + P[i]/(rho[i] * c**2)) / denominator
    P[i+1] = P[i] - (G * m[i] * rho[i] * GR[i]) / r[i]**2
    r[i+1] = r[i] + dr
    
    # 写入当前步长的所有数据
    with open("radial.out", "a") as f:
        line = f"{r[i+1]},{rho[i+1]},{P[i+1]},{m[i+1]},{GR[i+1]}\n"
        f.write(line)

关键修改说明 #

1. 统一文件写入方式：放弃混合使用np.savetxt和手动文件写入，改用with open()以追加模式（"a"）逐行写入，避免文件覆盖和指针冲突。
2. 修正循环范围：将for i in range(1, len(P))改为range(1, N-1)，确保i+1不会超出数组最大索引N-1，避免IndexError。
3. 替换while循环为条件判断：将嵌套的while P[i] > 0改为if P[i] <=0: break，避免无限循环，同时在压力小于等于0时终止计算，符合物理意义。
4. 增加分母校验：在计算GR[i+1]时检查分母是否合法，避免除以0或负数导致的计算错误，防止程序崩溃。
5. 使用with语句管理文件：with语句会自动处理文件的打开和关闭，确保缓冲区内容写入磁盘，无需手动调用close()和flush()。

内容的提问来源于stack exchange，提问作者Trissa Shamp