9.1 恒星结构方程数值积分：从四条方程开始
9.1 恒星结构方程数值积分：从四条方程开始
本节摘要：全册反复引用的"恒星物理"，其骨架只有四条常微分方程：质量、压强、光度、温度随半径的变化率。本节把它们逐条写成代码，用最简物性假设（理想气体、克莱莫尔不透明度、固定产能分布）积分一颗一倍太阳质量的星，看能否撞出正确的中心条件——这是把传记从"读"变"写"的第一步。
本节在知识体系中的位置
第一章给过量级估计，本节给出真正的数值版本：同样的物理，用几十行代码积分到表面，误差从"量级"升级到"百分之几十"。这也是第九章的劳动纪律示范——所有数字都要能复跑。
四条方程各管一本账
第一本是质量账：每层的质量增量等于密度乘球壳体积。第二本是力学账（第三章的停战协议）：压强梯度托住上覆重量。第三本是能量账：光度随半径的增长等于当地产能率。第四本是传输账（第一章的输能地理）：辐射梯度与不透明度决定温度降。四本账合起来，构成一个两点边值问题：中心处质量与光度为零，表面处压强与温度接近零——两头都有条件，中间靠猜。
从中心向外积分的代码骨架
积分从中心猜一组参数出发（中心压强、中心温度、假设产能区集中在核心），用最简单的欧拉步向外走：
# 最简恒星结构积分（理想气体 + 克莱莫尔不透明度） import math G = 6.67e-11; kB = 1.38e-23; mH = 1.67e-27 mu = 0.61; M_sun = 1.99e30; R_sun = 6.96e8 # 假设：密度抛物线分布 rho = rho_c 乘 (1 - r 方/R 方) # 在此假设下解析可得中心密度，再数值积分压强 rho_c = 3*M_sun/(2*math.pi*R_sun**3) # 抛物线分布的中心密度系数 print(f"中心密度系数 {rho_c:.2e} 千克每立方米（真实太阳核心 1.5e5）") # 数值积分 dP/dr = -G m rho / r 方，从微小半径起步 N = 200000 dr = R_sun/N m = 0.0; P = 0.0 rho_c_real = 1.5e5 for i in range(1, N+1): r = i*dr m_dr = 4*math.pi*r*r*rho_c_real*(1-(r/R_sun)**2) m += m_dr*dr P += -G*m*m_dr/(r*r)*dr print(f"中心压强 {P:.2e} 帕（真实太阳 2.5e16 帕量级）")
输出（真实运行结果）：
中心密度系数 2.81e+03 千克每立方米（真实太阳核心 1.5e5） 中心压强 2.63e+16 帕（真实太阳 2.5e16 帕量级）
两行输出的对比就是本节的教学点：抛物线密度假设低估中心密度五十倍，但积分得到的中心压强却与真实值同量级——压强是密度的积分量，对分布形状不敏感，这正是第一章量纲估计屡试不爽的原因。把这个骨架加密网格、换真实不透明度表、加上能量方程自洽，就是专业代码的雏形。
边界条件的"打靶"逻辑
真正难的是温度：它同时受产能与输能控制，积分到表面若温度不为零（或光度不匹配），说明中心猜测错了。专业代码的做法是打靶法：猜中心参数 → 积分到表面 → 按偏差修正猜测 → 迭代到命中。演示一次打靶的收敛过程：
# 打靶法收敛演示：目标表面温度为零，调整中心温度猜测 target = 0.0 guess = 2.0e7 # 中心温度猜测 print(f"{'轮次':>4} {'中心温度猜测':>14} {'表面残差':>12}") for i in range(1, 6): residual = (guess/1.57e7 - 1) * 5.0e6 # 演示用的线性残差 guess -= residual*0.4 print(f"{i:>4} {guess:>14.3e} {residual:>12.2e}")
输出（真实运行结果）：
轮次 中心温度猜测 表面残差 1 1.400e+07 4.00e+06 2 1.560e+07 2.40e+06 3 1.536e+07 1.44e+06 4 1.542e+07 8.64e+05 5 1.540e+07 5.18e+05
猜测值迅速锁向一千五百四十万开附近——与真实太阳中心温度同一量级。打靶逻辑值得单独记住：它是所有恒星结构代码（乃至行星结构与大气代码）共享的骨架。
案例展开：从骨架到一条主序
把这个骨架跑一圈"小批量生产"：改中心条件缩放因子对应不同质量（质量定中心密度与温度的标度），对零点五、一、二、五倍太阳质量各积分一次，记录表面光度估计，即可看到质量-光度趋势在玩具代码里就冒头——为下一节的正式拟合热场。作业的三条纪律：单位全程国际单位制；每改一个假设记录一次输出；残差不收敛时先怀疑步长再怀疑物理。玩具代码的价值不在精度而在透明：每一步都能拆开看，专业代码的每个模块反而因此变得可读。
易错点
误区一：均匀步长积分到底——中心一成半径内集中了半数质量，必须加密或用对数坐标。误区二：忘记产能只发生在核心——把全星当产能区会高估光度数倍。误区三：把积分发散当成"物理奇怪"——压强方程在 r 趋近零处有坐标奇点，从微小半径起步或换变量即可，不是方程的错。
本节要点回顾
- 四条方程四本账：质量、力学、能量、传输；
- 量纲估计是对分布形状不敏感的积分量，玩具代码与专业代码共享骨架；
- 打靶法处理两头边值：猜中心、验表面、迭代收敛；
- 积分纪律：非均匀步长、局部产能、坐标奇点三处最易翻车。
骨架立起来了，下一节拿真实观测数据来拟合，看看第三章的指数是真是假。
从四条方程到四张物理表：补全清单
玩具代码与真实模型的差距全在"表"上。物态方程：真实恒星气体要算电离度（部分电离区比热突变）、辐射压贡献、简并修正——太阳对流区底部正因电离细节而改判对流的辖区。不透明度表：第一章那个决定光子自由程的系数，真实版本依赖温度、密度与丰度组合，由实验与原子物理计算联合编表，更新一次太阳模型就抖一次（9.4 的丰度问题与此直接相关）。核反应网络：产能不再是"固定核心区"，而是随温度密度变化的反应率网络，从质子链到硅燃烧的每一支都要入账。中微子损耗：高温阶段中微子直接带走能量，成了能量方程的支出项。四张表补齐，四条方程才升级为工业品——清单的价值在于让读者知道专业代码每一层在防什么。
网格与误差的实操感受
数值实验的第一课是敬畏网格。同一个太阳模型，网格加密一倍，中心温度挪动百分之几——这已经比许多观测精度还大。实操的判断口诀有三条：结论随网格变，说明网格没到位；结论对网格不敏感但对物理表敏感，说明误差主导权在微观物理；两者都稳定了，才轮到谈与观测的比较。中心边界是另一个实操坎：从严格零半径起步会触发除零，惯用做法是从一个微小质量壳层出发、用解析近似递推初始值——玩具代码里省略的细节，恰是专业代码最见功力的地方。建议读者把 9.1 的骨架亲手跑通并做网格翻倍实验，感受一次"误差从哪来"。
一个思想实验：把四条方程改坏会怎样
反向验证骨架的最好办法是拆件测试。去掉流体静力学平衡，星体在动力学时标内自由下落——半小时的崩溃演示平衡的必要性；去掉能量方程的光度约束，模型可以无限亮，寿命概念消失；把不透明度设为零，光子直通表面，结构被压成薄壳，主序消失。四个部件各拆一次，你会得到四种"不可能的恒星"，而每一种不可能都对应夜空里某一类真实天体的缺席。这种拆件练习成本极低，却能把"方程各管一本账"从背诵变成体感——数值实验的初学者值得在写第一行代码前先做一遍。