恒星讲义6

第六章 恒星的结构和演化模型
第一节 恒星结构演化基本方程组
1. 基本假设条件
恒星是一个巨大的气体球，其自身的引力和气体的压强相平衡而维持在一个相对稳定状态。

在其中心区域，温度非常高，由此引发的热核反应提供了恒星不断向外辐射的能量。

在恒星内部，辐射传热过程将热核反应在恒星中心释放的热量转移到恒星的表面，同时也确定了温度在恒星内部的分布。

但是在某些区域，热核反应释放的能量非常大，或者是物质的不透明度非常大，使得温度梯度超过绝热温度梯度而引发对流运动。

由于运动的尺度非常大，流动将以发展非常充分的湍流的形式出现。

在对流区内，湍流非常有效地传递着热量，温度的分布也由对流传能效率决定。

恒星结构演化问题的简化近似条件：
(i)恒星是球对称的结构，不考虑自转、非球位型磁场的效应；
(ii)恒星处于准稳态，其结构满足流体静力学平衡条件；
(iii)除了最外部的大气层，在恒星内部物质处于局部热动平衡状态。

在上述近似假设下，恒星结构演化的基本方程组是：
2. 质量连续性方程
取一个半径为r 的球壳，其中包含的物质的质量满足方程：
dr
r dM
r
ρπ2
4=
其中M r 是半径为r 的球面内包含的质量，球壳的厚度为dr 。

一般取M r 为自变量比较方便。

于是，质量连续性方程可以写成为：
ρ
π2
41r M
r r
=
∂∂
3. 流体静力学平衡方程
上述球壳受到的重力应该与作用在球壳上的压力平衡，以保持球壳为静止状态：
dP
r gdM
r
2
4π=-
于是，流体静力学平衡方程写成为：
4
2
44r
GM r
g M
P r r
ππ-
=-
=∂∂
其中，重力加速度g ＝GM r /r 2。

4. 能量守恒方程
设单位质量物质的热核产能率为εN ，中微子能量损失率为εν，而恒星的光度为L r ，那么能量守恒定律要求：
r
r N M
L t
s T ∂∂-
-=∂∂ νεε
利用热力学关系，可以将上式的左边改写成为：
()g ad P ad P P
t P P T c t P t
T
T c t P
t T c t s T ερδ-=∂∂∇-∇=
⎪⎭⎫ ⎝⎛∂∂∇-∂∂=∂∂-
∂∂=∂∂ln ln
这一项常常被称为引力能。

于是，能量守恒方程变成：
g
N P
N r
r t
P t
T c M
L εεερδεενν+-=∂∂+
∂∂--=∂∂
5. 能量传递方程
在恒星内部，能量向外传递有三种可能的方式： (a)辐射传热
根据以前的讨论，在光学厚的区域内，辐射热通量为：
r
T acT F r
L R
R R ∂∂-
==ρκ
π344 3
2
(b)电子气体热传导
当自由电子气体处于强简并状态时，电子气体的费米能量有可能非常高，使得电子气体有很大的平均速度。

这时电子气体的热传导将称为向外传递能量的主要方式。

引入一个等效不透明度κD ，可以将热传导通量写成为和辐射热通量类似
的热扩散形式：
r
T acT F r
L D
D D ∂∂-
==ρκ
π344 3
2
(c)对流传热
根据前面的讨论，湍流对流的热输运通量为：
''4 2
T u c F r
L l P C C ρπ==
于是，总光度可以写成为：
r
T acT
T u c r T acT T u c F F F r
L l P D R l P C D R r ∂∂-=⎪⎪⎭
⎫ ⎝⎛+∂∂-
=++=ρκρκκρ
ρπ34''1
134''4 3
3
2
其中，κ是包括辐射和传导作用的总不透明度。

利用温度梯度的定义，得到：
''434
T u c acT
H l P P R ρρκ-
∇=∇
在辐射平衡区，上式右边只有第一项。

在对流区，计算出湍流对流通量后，温度梯度也就知道了。

利用温度梯度，可以得到：
∇-
=∂∂P
r T GM M
T r r
2
4π
6. 元素丰度演化方程
元素数丰度演化方程的一般形式可以写成为：
∑=
∂∂l
k l
k kl
i Y Y R
t
Y ,
在恒星结构演化计算中，通常使用质量丰度X i ＝Y i A i ，其中A i 是i 元素的原子量。

于是，元素丰度演化方程写成为：
∑
=
∂∂l
k l
k kl l
k i i X X R A A A t
X .
如果热核反应发生在对流区，一般假定对流运动的输运效率非常高，于是对流区内的物质是完全混合的。

此时元素丰度的变化为：
⎰⎰∑
=
∂∂c
r
c r
l
k l k kl l
k i i dM
dM
X X R A A A t
X .
第二节 边界条件和初始条件
基本方程组中，有四个方程含有对位置的偏导数，于是需要4个边界条件。

1. 中心边界条件
在恒星的中心，M r ＝0。

显然必须满足条件：
0=r
0=r L
但是，恒星中心的温度和压强则是自由的。

2. 表面边界条件
恒星表面的情况较复杂。

严格地说，随着高度的增加，物质的密度连续地下降。

从这个意义上来讲，恒星没有确定的表面。

我们所看到的恒星圆面，只是在此波段光深为1附近的地方。

不同波段光深不同，从而看到的也是太阳大气中不同深度的地方。

这使得外边界的位置无法确定，也就无法规定此处必须满足的条件了。

精确的作法是将恒星内部结构模型与恒星大气模型相结合，将内部结构确定的重力加速度和光度作为确定大气结构的参数，反过来将大气底部的温度压强作为内部结构的外边界条件，最终寻找到自洽的解。

作为一个近似的作法，我们以灰大气解作为大气模型，来寻找与内部结构自洽的外边界条件。

灰大气温度分布为：
⎪⎭⎫ ⎝
⎛
+=
3243
44
τeff T T
其中，有效温度T eff 与恒星光度L 和半径R 的关系是：
4
24eff T R L σπ=
在外边界处τ＝0,于是得到： 4
421eff
T T =
4
463eff
T a T
a
P ==
于是，恒星结构方程组构成了一个两点边值问题。

3. 初始条件
恒星热状态的变化和元素丰度的变化涉及到对时间的偏导数，因此需要给出恒星内部结构和元素丰度的初始分布。

初始丰度的选取可以根据我们感兴趣的对象的性质直接给定。

一般选择一组丰度比，如星族I 型恒星常取太阳的元素丰度比，而星族II 型恒星一般选择一个较小的金属丰度Z 。

初始结构模型的确定就没有那么简单了。

一种有效的方法是利用恒星在主序演化时标非常长，可以近似认为其处于定常状态而忽略恒星热状态的变化，求解静态结构模型得到初始状态。

还有一种方法是利用基本方程组的近似解，如Emden 模型来模拟恒星早期的某些状态，从而得到近似的初始条件。

4. 解的唯一性问题
虽然在恒星内部物理量(元素丰度除外)显然是连续和有界的，但是由于边界条件是分别在两个边界上给出的，因此恒星结构演化问题不满足Lipschitz 唯一性条件，也就是说可以存在多个解。

这一点已经被大量的计算所证实。

但是，这些多重解之间的差别相当大。

因此，在HR 图上一个相对来说很小的范围内，基本方程组的解是唯一的。

第三节 Virial 定理
一个满足流体静力学平衡的自引力系统，其所包含的引力能与内能之间存在着关系。

这种关系就称为Virial 定理。

将流体静力学平衡方程乘以4πr 3,然后在整个恒星内部积分得到：
⎰
⎰⎰
⎰
⎰
-=-
==
-=-
M
r
M
M M
M
r
r
M
r
r dM
P
Pdr
r P
r dP r dM
r
GM dM r r
GM 0
2
3
30
3
43 124444ρ
πππππ
对于单原子理性气体，如果忽略辐射对内能的贡献，其单位质量的内能为：
ρ
μ
2323P T u =
ℜ=
代入上式，可以得到：
2=+T G E E
其中，总势能E G 和总内能E T 分别定义为：
⎰
-=M
r
r
G dM
r
GM E 0
⎰=
M
r
T dM
P
E 0
23ρ
上式表明，恒星具有的引力能是其内部包含的内能的两倍。

另一方面，恒星的总能量为-E T ，这表明恒星处于引力作用下的束缚态。

当恒星内部没有热核能源时，它辐射掉的能量只能来自于其所具有的总能量，即：
()
dt
dE dt
dE dt
E E d L G T T
G 21-
==
+-
=
由于引力能是负的，于是恒星必须收缩，其释放的引力能一半加热了恒星物质，一半从恒星表面辐射掉。

对于一般的理想气体，利用热力学定律可以得出：
ρ
γP
u 11
-=
其中γ＝c P /c V 。

如果γ是一个常数，那么Virial 定理可以写成为：
()013=-+T G E E γ
恒星的总能量为：
()T
T G E E E 43--=+γ
由此可以得到恒星结构稳定性的一个重要判据，即：只有当γ>4/3时，恒星才能够处于流体静力学平衡。

第四节 几个重要的时标
恒星内部发生的各种物理过程是以不同的速度进行的。

下面三个时标是与恒星结构演化问题关系最密切的。

1. 动力学时标
当流体静力学平衡遭到破坏，恒星将进行动力学调整。

这时，运动方程将取代流体静力学方程来描述所发生的运动过程。

而运动所需要的典型时间，被称为动力学时标。

考虑两种特殊情况来说明动力学调整所需要的时间。

(a)自由坍缩过程
如果压强消失，物质将做自由落体运动。

显然，这时运动方程可以写为：
g dt
r d -=22
在做数量级分析时，取各个量的典型值，于是得到：
GM
R
g
R ff
3
=
=
τ
其中，M 是恒星总质量，R 是恒星半径。

τff 称为自由流体时标。

(b)绝热声传播过程
如果流体静力学平衡结构受到小的扰动，流体元将围绕其平衡位置振动。

设扰动沿z 方向，并且是绝热的，则扰动方程组可以写为：
z
P z
P t s ∂∂⎪⎪⎭⎫ ⎝⎛∂∂-
=∂∂-
=∂∂'
1'
1v 00ρρρρ 0
v '
10=∂∂+∂∂z
t
ρρ
其中下标“0”代表平衡结构的物理量，带撇的量是对平衡量的扰动。

将上述方程组合并，得到绝热声波的波动方程：
22
2222
022
v v '1v z c z
P
t z P t
s s s ∂∂=∂∂⎪
⎪⎭⎫ ⎝⎛∂∂=∂∂∂⎪⎪⎭⎫ ⎝⎛∂∂-
=∂∂ρ
ρρρ
其中绝热声速c s 为：
ρρ
P
P c s
s 12
Γ=⎪
⎪⎭⎫
⎝⎛∂∂=
而Γ1是绝热指数。

考虑声波穿越恒星所需要的时间τs ：
ff
R
R
R
s
s GM
R
dr u
dr P
c
dr
τ
ρτ≈Γ≈
Γ=
Γ=
=
⎰
⎰
⎰13
10
10
2323
可以看到，这两种动力学过程对应的时标是相似的。

例如： 太阳： τff ≈27分钟 红巨星(M ＝M ⊙，R ＝100R ⊙)： τff ≈18天 白矮星(M ＝M ⊙，R ＝0.02R ⊙)： τff ≈4.5秒
2. 热时标（Kelvin-Helmholtz 时标）
恒星处于局部热动平衡状态。

当流体偏离局部热动平衡时，将发生驰豫过程使恒星回到局部热动平衡状态。

这一过程的典型时标就称为热时标，有时也称为Kelvin-Helmholtz 时标。

这个时标也是恒星将其全部热能损失掉所需要的时标。

于是，热时标τKH 定义为：
RL
GM L E L E G T KH 2
=≈=
τ
对于太阳，τKH ＝3×107年。

3. 核时标
热核反应维持着恒星稳定地发光。

在其中心，物质从轻的核素转变为重的核素。

一旦核燃料耗尽，恒星将失去光芒而死亡。

于是，将恒星耗尽其内部核燃料所需要的时间叫做核时标τN 。

显然，
L
E N N =
τ
其中E N 是恒星内部核燃料可以释放的能量的总和。

氢是恒星内部最丰富的核燃料，其燃烧维持恒星绝大部分发光时间。

对于氢燃烧，其热值为1.5×1018erg/g 。

因此，核时标大约为：
L
L M
M N ⊕⊕
=11
10
τ 年
第五节 多方模型
1. 多方关系
如果压强和密度之间存在关系：
n
K P 11+=ρ
那么只需要2个方程就可以确定恒星的内部结构。

这种模型被称为多方模型。

满足多方过程的情况：
(1)完全简并气体：在非相对论情况下，
P ∝ ρ5/3；n ＝3/2
在相对论情况下，
P ∝ ρ4/3；n ＝3
(2)绝热对流区：对于完全电离理想气体，
5
2ln ln =
∇==
∇ad P
d T d
于是，利用状态方程可以得到：
P ∝ ρ5/3；n ＝3/2
(3)气体压强与总压强之比β为常数区域：这时辐射压强为
()4
31T
a P P R =
-=β
求出T 后代入状态方程，得到：
()3
/43
/13
/413ρ
βμβ⎥⎦
⎤
⎢⎣⎡-⎪
⎪⎭
⎫ ⎝⎛ℜ=a P ；n ＝3
(4)等温核：这时状态方程为
ρμ
T
P ℜ=
；n ＝∞
2. Emden 方程
利用多方关系，当K 为常数时，引入引力势Φ为：
()
dr
d n K dr
d n K dr dP
g dr d n
n
1
1
11111ρρρρ+-=⎪⎭⎫ ⎝
⎛
+-=-==Φ
-
积分上式得到：
()n
n K 1
01ρ
+-Φ=Φ
在恒星表面处，显然有ρ＝0和Φ＝0，于是Φ0＝0。

于是，上式可以写成为：
()n
n K ⎪⎪⎭
⎫
⎝⎛+Φ-=1ρ
将此结果代入引力势的泊松方程，得到：
()n
n K dr d r dr d r ⎪⎪⎭
⎫ ⎝⎛+Φ-==⎪⎭⎫
⎝⎛Φ144122ππρ
引入变量：
()
r n K G Ar z n
n c
1
14-+=
=ρπ
n
c c
u 1
⎪⎪⎭
⎫ ⎝⎛=ΦΦ=
ρρ
其中下标c 代表中心的值，上述方程可以变换为：
0222
=++
n
u dz
du z dz
u d
上述方程被称为指数为n 的Emden 方程。

边界条件为：
恒星中心z ＝0时：u=1,u '=0 恒星表面z ＝z n 时：u ＝0
3. Emden 方程的解
由于z ＝0是方程的一个奇点，设n 为有限数，将方程的解用级数展开为：
+++++=4
43
32
211z a z a z a z a u
代入方程中，比较同次幂的系数，可以得到：
++
-
=4
2
120
611z n z u
由上式可以看到，当z 很小时，u 有极大值1。

当n 为0、1、5时，Emden 方程有解析解：
n ＝0时：2
611z
u -=
n ＝1时：z z u sin =
n ＝5时：2
3
111z
u +
=
对于任意的正数n ，方程只能通过数值方法求解。

另外，只有当n <5时，方程的解可以满足外边界条件，即模型的半径是有限的。

下表给出了一些z n 的值。

由表中数据可以发现，n 值越大，z n 的值也越大，表明模型的半径越大，而此时模型的中心聚度却越高。

4. 恒星的多方模型
在给定恒星的质量M 和半径R 后，可以由Emden 方程的解构造出恒星的结构模型，即确定温度T 、压强P 和密度ρ的分布。

在根据恒星内部物质的具体情况确定指数n 后，可以得到：
n
c u ρρ=
n
n
c
n
u
K K P ++
+
==11111ρρ
可以看到，只要能够确定中心密度ρc ，就可以确定恒星的内部结构。

恒星的总质量为：
n
n
n
z z c z c z n c R
n c R
dz du z R dz du z d z r dz z u z r dr r u dr r M =⎪
⎭⎫ ⎝⎛
-=⎪⎭⎫
⎝⎛-⎪⎭⎫ ⎝⎛=⎪
⎭⎫ ⎝⎛===⎰⎰⎰⎰14 444430
230
2
302
02πρπρπρπρρπ
于是，中心密度可以写成为：
1
1
3
314-=-=⎪⎭⎫
⎝⎛-=⎪⎭⎫ ⎝⎛-=n
n z z z z c dz du z dz du z R M
ρπρ
其中ρ是恒星的平均密度。

利用这些关系，可以估算恒星中心的物理状态。

例如，太阳的平均密度为1.4g/cm 3。

取n ＝3,查表得到z n ＝6.9和ρc 为54.2×1.4＝76.4g/cm 3。

于是，A ＝z n /R ⊙＝9.9×10-11。

这样可以算出K ＝3.85×1014。

于是，P c ＝1.24×1017dyn/cm 2和T c ＝1.2×107K 。

第六节 恒星结构演化模型的数值求解方法
一般说来，基本方程组是一组互相耦合的非线性微分方程组。

对于相对比较准确的物质状态方程和不透明度，分析方法求解基本方程组是不可能的。

数值方法求解微分方程组是一种非常有效的方法。

数值解虽然不能给出刻画恒星结构演化问题的普遍规律的结论，但是可以帮助我们准确了解每一个特定恒
星的演化图像，进而总结出恒星演化的总体概貌。

1. 差分方案
数值求解微分方程组的基本思想就是用差分代替微分，从而将微分方程变为代数方程，求解代数方程而得到原来微分方程的近似解。

(1)对形如
)(x f dx
dy =
的微分方程，在(x 1,x 2)区间上积分上述方程，得到：
⎰
=
-2
1
)(12x x dx x f y y
构造无限可微函数f(x)的多项式逼近：
()()
+-+
-+=2
11111)(''!
21)(')()(x x x f x x x f x f x f
当积分步长∆x ＝x 2-x 1充分小时，有限项逼近就可以给出任意指定精度的解函数y 。

特别重要的是只采用零阶近似f(x 1)，差分的精度是O(∆x)。

(2)对形如
),(y x f dx
dy =
的微分方程，在(x 1,x 2)区间上可以构造如下差分方程： (i)显式格式
),(111
212y x f x x y y =--
(ii)隐式格式
),(221
212y x f x x y y =--
(iii)中心差分格式
[]),(),(2
122111
212y x f y x f x x y y +=--
在上述差分格式中，显式格式和隐式格式的误差是O(∆x)，而中心差分格式的误差为O(∆x 2)。

显式格式与中心差分格式会传递截断误差，且不稳定；而隐式格式不传递截断误差，且是无条件稳定的。

(3)对形如
),(y x f dx
dF
dx dy dx d -=-
=⎪⎭⎫ ⎝⎛λ
的微分方程，它在物理上代表了扩散过程，而F 代表了扩散流的通量。

于是，一个重要的条件是：在每一格点左右两边计算差分时，F 的计算必须完全相同，这样可以保证不会因为差分而引入虚假的热源。

这种格式被称为守恒格式，因为在没有热源的情况下，它保证了热通量是一个常数。

一般将上述方程分为两个一阶微分方程来处理，即：
),(y x f dx dF = λ
F
dx
dy -
=
并且在(x 1,x 2)区间上，采用中心差分格式来近似原微分方程：
[]),(),(2
122111
212y x f y x f x x F F +=--
⎥⎦
⎤
⎢⎣⎡+-
=--22111
21221λλF F x x y y
由此可以解出F ，并且得到：
2
112
2121
21211
122λλλ+-----=
f x x x x y y F
2
112
1121
21221
1
22λλλ+
-+---=
f x x x x y y F
同理，在(x 2,x 3)区间上，差分给出2格点位置处的热通量为：
3
223
3232
32321
122λλλ+-----=
f x x x x y y F
于是可以看到，在2格点左右两边，通过2格点的热通量是不同的，因此上述差分方法不是守恒格式。

在某些场合，这种非守恒格式会带来计算的误差，甚至是不稳定。

一种合理的差分方法是，将热通量设定在分解面上。

于是，对于n 格点，积分原微分方程得到：
n n x
x
n n n n x f dx y x f F F dx dy dx dy n n ∆≈=-=⎪
⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛⎰
+
-
-++-2
12
1),(2
1212
1
2
1λλ
而热通量F 的差分取为：
n
n n
n n n x x y y F
---=+++
+112
12
1λ
这种差分方法称为二阶导数的中心差分格式。

它不但有O(∆x 2)的精度，而且还保证了分解面处的热通量是由分解面两边的物理量计算得到的，因此满足守恒格式的条件。

2. 恒星结构演化的离散化方程组
将恒星由外边界到中心分为M 层，共M+1个格点，其中中心记为0格点,而外边界记为M 格点。

取变量组为：
⎪⎭
⎫
⎝⎛
-
=M M r 1ln ξ （自变量）
P y ln 1=
T y ln 2=
⎪⎭⎫ ⎝
⎛
+=L L y r α1ln 3
（α是一个常数）
r y ln 4=
采用差分格式代替微分，得到离散化方程组：
()1,,11
1,,,,,----=--j k j k j j j j j i j i y y x x f x x y y ， i,k ＝1,2,3,4,j ＝1,2,…，M
注意，能量守恒方程中对时间的微商一般采用隐式格式离散化：
t
P P t
T T c t
P t T c P
P
g ∆-+
∆--≈∂∂+
∂∂-=00
ρ
δρδε
其中下标0代表上一时刻的值，∆t 是演化的时间步长。

定义函数G i,j 为：
()()0
,,,1,,111,,,=---=----j k j k j j j j j i j i j i y y x x f x x y y G
原方程的解可以由上述代数方程组的解近似。

上述代数方程共有4M 个，而未知数有4(M+1)个。

其余4个方程由边界条件补充。

首先是在内边界上存在边界条件r ＝0和L r ＝0。

显然，r ＝0的条件是无法直接使用的，因为中心是基本方程组的一个奇点。

考虑0和1格点范围内中心球的性质。

在恒星中心附近一个很小的区域内，密度ρc 和产能率εc 可以被看成是常数。

于是，根据基本方程组可以得到：
1/31/3
43r
c
M
r ⎪⎪⎭
⎫ ⎝
⎛=πρ
3/23
/14
342r
c
c M
G P P ⎪⎪⎭
⎫
⎝⎛-=πρ
r
c r M
L ε=
3/23
/14342r
c
c c a
d c M
P T G T T ⎪⎪⎭
⎫
⎝⎛∇-
=πρ
如果恒星中心是对流区
3/23
/223438r
c c
c c M
acT
T T ⎪⎪⎭
⎫ ⎝
⎛-
=πρκε
如果恒星中心是辐射平衡区
这些关系将中心附近的物理量与恒星中心的边界条件联系起来了，如果给定了中心的压强P c 和温度T c ，则中心附近的解就已经确定了。

这样，0格点处的未知数只剩下P c 和T c 2个。

外边界处边界条件的处理更加复杂。

在光学深度很小的区域内，温度的分布服从灰大气模型；而在光学深度很大的地方，扩散近似被用来确定温度的分布。

因此，需要寻找一个自洽的方法来确定温度在不同区域的光滑的过渡。

一般来说，设定某一个光学深度处为内部解和大气解的拟合点，例如τ＝2/3，并且要求在拟合点处内外解光滑的衔接。

设给定恒星的L 和R ，于是在τ＝0处有： σπ2
4
4821R L T T eff =
= σ
π24243
R aL
T a P =
=
根据灰大气温度分布和流体静力学平衡方程：
⎪⎭⎫ ⎝⎛+=
321632
4
τσπR L
T
κ
τ
2
R GM d dP =
可以得到恒星大气的温度和压强分布。

显然，拟合点处的温度T M 和压强P M 是恒星光度L 和半径R 的函数。

利用这组（T M ，P M ）作为边界条件，那么格点M 处只剩下L 和R 是未知数。

于是就可以封闭内部区域的差分方程组。

在拟合点处每次求解大气模型很费时间，而且容易造成计算的不稳定性。

如果内部迭代给出的L M 和R M 修正不大，可以认为P M 和T M 是L M 和R M 的线性函数：
P M
eff P M P M c T b L a P ++=,ln ln ln T
M
eff
T M T M c T b L a T ++=,ln ln ln
于是，可以在HR 图上取3个相距很近的点（或者3组L,T eff ），通过求解灰大气模型得到在大气底部的P M 和T M ，代入上式就可求出拟合系数a 、b 和c 。

上述关系式可以用来作为外边界条件。

3. 解非线性代数方程组的牛顿迭代法—Henyey 方法
对于一个非线性方程：
0)(=x f
在任意x 0处进行展开：
()
+-+=000)(')()(x x x f x f x f
就可以构造迭代格式：
'
1f f x x n n -
=+
来寻找满足原方程的根。

这种方法被称为牛顿迭代法。

它是收敛速度最快的迭代算法。

对于恒星结构演化问题，离散化代数方程组为：
()0,1,,,=-j k j k j i y y G ，
i ＝1,2,3,4；j ＝1,2,……，M
进行展开得到：
04
1
1,1
,,,4
1
,,,=∆∂∂+
∆∂∂+
∑∑
=--=k j k j k j
i j k k j
k j i j i y y
G y y G G
给定一组试探解y ，通过求解上述方程组就可以得到改正量∆y 。

将改正量与试探解相加，得到新的试探解。

反复进行上述迭代过程，直到改正量小于预先指定的精度要求，并且试探解对方程组的满足程度达到预先指定的精度要求，就得到了收敛的解。

4. 解线性代数方程组的追赶法
上述代数方程组的系数矩阵是一个分组对角矩阵，可以采用下述技巧进行求解。

当j ＝1时，方程的数目是4个，而未知数为：P 0,T 0,P 1,T 1,L 1,r 1。

于是，可以将其中4个的改正量用其余2个的改正量表示出来，例如：
1,11,21,11,11,10,1C y B y A y +∆+∆=∆ 1,21,21,21,11,20,2C y B y A y +∆+∆=∆
1,31,21,31,11,31,3C y B y A y +∆+∆=∆ 1,41,21,41,11,41,4C y B y A y +∆+∆=∆
当j ＝2时，方程数是4个，而未知数是8个。

但是，已经有2个（y 3,1和y 4,1）用另外2个（y 1,1和y 2,1）表达出来了。

于是独立的未知数也是6个。

可以用类似的方法将其中4个表达为其余2个的函数：
1,2
2,21,22,11,21,1C y B y A y +∆+∆=∆ 2
,22,22,22,12,21,2C y B y A y +∆+∆=∆
2,32,22,32,12,32,3C y B y A y +∆+∆=∆ 2,42,22,42,12,42,4C y B y A y +∆+∆=∆
对于2＜j ＜M ，依此类推。

当j ＝M 时，方程数是4个，而未知数是6个（P M-1,T M-1,L M-1,r M-1 ,L M ,R M ），其中有2个（L M-1,r M-1）又可以用另外2个（P M-1,T M-1）表达出来。

于是，未知数只剩下4个，可以完全求出。

利用得到的结果，可以将解的形式依旧写成为：
M
1,M ,2M 1,M ,1M 1,1,1C y B y A y M +∆+∆=∆- M 2,M ,2M 2,M ,1M 2,1,2C y B y A y M +∆+∆=∆-
M 3,M ,2M 3,M ,1M 3,,3C y B y A y M +∆+∆=∆ M
4,M ,2M 4,M ,1M 4,,4C y B y A y M +∆+∆=∆
于是，代入前面得到的递推关系，就可以得到完全格点上的改正量。

5. 元素丰度的演化
为简单起见，元素丰度随时间的演化一般采用显式格式进行离散化：
∑=
∆-≈∂∂l
k l k kl i i i Y Y R
t
Y Y t
Y ,0,0,0
,0
,
如果对元素丰度与热核反应的自洽性要求很高时，也可以采用隐式格式进行离散化。

这时，元素丰度将与其它物理量一起迭代求解。