点堆动力学方程的数值解法

1、 报告人： 学号： 对于这些不同的差分方程是否都同样有效，同样可靠，而且能得到同样的计算结果呢？ 答案是否定的。事实上，性质不同的，数值结答案是否定的。事实上，性质不同的，数值结果不同。果不同。差分格式与原微分方程的差分格式与原微分方程的相容性相容性利用差分格式求解的利用差分格式求解的收敛性收敛性差分格式的差分格式的稳定性稳定性如何判断和分析差分方程有效性和可靠性如何判断和分析差分方程有效性和可靠性就成为非常必要和现实的问题了。就成为非常必要和现实的问题了。1、相容性、相容性10 ( , ; ) (2)nnnnyytytty 如果如果00()( )limlim( , ( ) ; )tty tt

2、y tty ttt 则称差分格式与微分方程是相容的。则称差分格式与微分方程是相容的。差分方程相容性是讨论当差分方程相容性是讨论当步长无穷小步长无穷小时，差分方程逼近于微时，差分方程逼近于微分方程的程度，因此，相容性是讨论分方程的程度，因此，相容性是讨论差分方程和微分方程的差分方程和微分方程的关系关系。相容的条件显然是相容的条件显然是ytfytnn,0;,2、收敛性、收敛性 ( 0)(0)1yyy tye此方程的解析解是此方程的解析解是 （1）差分方程收敛性表示差分方程数值解和微分方程精确解逼近程度，）差分方程收敛性表示差分方程数值解和微分方程精确解逼近程度，只有在只有在 差分方程收敛于微分方程

3、时，差分方程解才可能是微分方程精确差分方程收敛于微分方程时，差分方程解才可能是微分方程精确解。解。（2）差分方程相容性是差分方程首先要满足的，差分方程相容性）差分方程相容性是差分方程首先要满足的，差分方程相容性是收敛性的必要性条件，但并不是充分条件。是收敛性的必要性条件，但并不是充分条件。例：例：3、稳定性、稳定性计算误差如果在运算过程中是逐渐缩小的，计算误差如果在运算过程中是逐渐缩小的，那么这种方法就叫做稳定的；那么这种方法就叫做稳定的；如果这些误差在运算过程中是逐渐放大的，如果这些误差在运算过程中是逐渐放大的，则称为不稳定的。则称为不稳定的。所谓稳定性问题，就是误差的积累是否能得到所谓稳定

4、性问题，就是误差的积累是否能得到控制的问题。控制的问题。4、点堆方程的刚性、点堆方程的刚性刚性（刚性（stiff）微分方程）微分方程1( )( )( )( )( )( )( ) 1,2,.,Iiiiiiiidn ttn tC tdtdC tn tC tiIdt堆内缓发中子先驱核的衰变常数以及中堆内缓发中子先驱核的衰变常数以及中子一代时间之间相差数个数量级子一代时间之间相差数个数量级）。）。 7107210721jtjtttjeAneAeAeAntn正反应性，当 k = 1+0.001 时1=0.018817 2=-0.014160 3=-0.062221 ,4=-0.18664 5=-1.20

5、97 6=-3.7650 7=-55.5341=0.0188177=-55.534利用可视化方法，形象地说明点堆动力学方程刚性问题的形利用可视化方法，形象地说明点堆动力学方程刚性问题的形成及来源。成及来源。，计算时间为，计算时间为100 s图图1：基于可视化方法反应性正阶跃输：基于可视化方法反应性正阶跃输入时入时0 100 s内中子密度变化内中子密度变化图图2 f1 f7,在不同坐标系下的变化曲线在不同坐标系下的变化曲线二、三性与点堆动力学方程的刚性二、三性与点堆动力学方程的刚性 4、刚性（、刚性（stiff）微分方程）微分方程00100100111y( )( )y( )( )( ) ( ),

6、( )( ),F(t)(1) (1)( )/00( )/0tTTIIIIId tF ttdty tyy tn t C tC tyn CCIItF t 其中，及是列向量，表示如下矩阵：二、三性与点堆动力学方程的刚性二、三性与点堆动力学方程的刚性 4、刚性（、刚性（stiff）微分方程）微分方程ii( )0F tE11/()/()IIiiiiiii 如果如果是这个矩阵的特征值，那么是这个矩阵的特征值，那么就是特征方就是特征方程程由此可推出：由此可推出：的根。的根。二、三性与点堆动力学方程的刚性二、三性与点堆动力学方程的刚性 4、刚性（、刚性（stiff）微分方程）微分方程二、三性与点堆动力学方程的

7、刚性二、三性与点堆动力学方程的刚性 4、刚性（、刚性（stiff）微分方程）微分方程1I如果111221IIIl 该方程的根全部是实数，并满足该方程的根全部是实数，并满足61( )( )( )( )(1 1)( )( )( ) 1,2,.,6iiiiiiidn ttn tC tdtdC tn tC tidt三、点堆动力学方程的数值解法三、点堆动力学方程的数值解法 满足的初值条件为：满足的初值条件为：000t=0( ),( )= , 1,2,.,6iitn tn C tCi三、点堆动力学方程的数值解法点堆动力学方程的数值解法 (1)后退后退欧拉方法欧拉方法2、点堆动力学方程的数值解法点堆动力学方

8、程的数值解法61(1)( )(1)(1)(1)(1, )( , )( )(1)( ) (1, ) (3-6) ( )( ) ( , ) 1,2,.,60,1,2,jn in iin if itC ijC i jjn ij C ijtf ij C i jji三、点堆动力学方程的数值解法点堆动力学方程的数值解法 61( (1)(1)( ( )(1)/(1)( ) (1)(1, )( , )/(1( ) ( )( ) ( , ) (37)( ) (0)(0, )( )1,2,.,60,1,2,jtin in itf itj n iC ijC i jtjf ij C i jj nCjjji三、点堆动力

9、学方程的数值解法点堆动力学方程的数值解法 (2) 梯形公式法梯形公式法61(1)( )( )(1)( )( )(1)(1)/2(1, )( , )( ) (1)( )(1, )( , )( ) 22( )( ) ( , ) (3-8) 1,2,.,60,1,2,jn in iiin if in if itC ijC i jj n in iC ijC i jjtf ij C i jji三、点堆动力学方程的数值解法点堆动力学方程的数值解法 61( )(1)(2( (1) ( )(1)2( (1)( )( ( )(1)(2( ) ( , )(1, ) 2( )( )( ) ( , ) (3-9)(

10、) (0)(0, )( )1,2,.,60,1,2,jt f it f itin in ititj n in itj C i jC ijtjf ij C i jj nCjjji 61( )( )( )( )(1)( )( )( ) 1,2,.,6iiiiiiidn ttn tC tdtdC tn tC tidt三、点堆动力学方程的数值解法三、点堆动力学方程的数值解法 满足的初值条件为：满足的初值条件为：000t=0( ),( )= , 1,2,.,6iitn tn C tCi三、点堆动力学方程的数值解法点堆动力学方程的数值解法 3、厄密特多项式方法、厄密特多项式方法将点堆方程简记为将点堆方程简

11、记为( )( )iiiiiiidnnCdtdb ta tCnCdt（2）三、点堆动力学方程的数值解法点堆动力学方程的数值解法 3、厄密特多项式方法、厄密特多项式方法 将时间离散化，并以将时间离散化，并以h ht2-t1t2-t1表示从表示从t1t1到到t2t2的时间步长，在这步长的两端，记的时间步长，在这步长的两端，记1122112211221122( ), ( )( ), ( )( ), ( )( ), ( )iiiiiiiinn tnn tCC tCC taa taa tbb tbb t（3）三、点堆动力学方程的数值解法点堆动力学方程的数值解法 3、厄密特多项式方法、厄密特多项式方法11

12、11222211 112222 , , iiiiiiiiiiiiiipa nCpa nCqb nCqbnC( ) )in tC t和 （进一步，记进一步，记它们分别是它们分别是 在步长在步长h两端的导数值两端的导数值 （4）三、点堆动力学方程的数值解法点堆动力学方程的数值解法 3、厄密特多项式方法、厄密特多项式方法11, tttthh即1( ) ( )( ) 1( ) ( )( ) iiiiiiidnanCh ddCbnCh dt为方便起见，作变量代换为方便起见，作变量代换于是于是t t的时间间隔（的时间间隔（t1,t2t1,t2）就变化为关于）就变化为关于的的时间间隔（时间间隔（0,10,1

13、）, ,故点堆方程可写成故点堆方程可写成 （5）三、点堆动力学方程的数值解法点堆动力学方程的数值解法 3、厄密特多项式方法、厄密特多项式方法（6）( ) ( )0,1 3inC设和在区间（）上可以 用 阶厄密特插值多项式表示：231234231234( )( )iiiiinmmmmCdddd所谓厄密特插值，就是要求在插值区间两端，该所谓厄密特插值，就是要求在插值区间两端，该多项式的函数值与导数值与被插值的函数值和导多项式的函数值与导数值与被插值的函数值和导数值相等。数值相等。 三、点堆动力学方程的数值解法点堆动力学方程的数值解法 3、厄密特多项式方法、厄密特多项式方法 （8）112132121

14、42121(2)3()()2()mnmhpmh ppnnmh ppnn 11213212142121(2)3()()2()iiiiiiiiiiiiiidCdhqdh qqCCdh qqCC 可以推出：可以推出： （7） 三、点堆动力学方程的数值解法点堆动力学方程的数值解法 3、厄密特多项式方法、厄密特多项式方法 （9）把（把（6 6）代入（）代入（5 5），并对），并对在（在（0,10,1）区间上积分，得到）区间上积分，得到112323211234123400112323211234123400( )() ()( )() ()iiiiiiiiiiiiiiinnh ammmmdhdddddCCh

15、 bmmmmdhddddd三、点堆动力学方程的数值解法点堆动力学方程的数值解法 3、厄密特多项式方法、厄密特多项式方法2i2 , Cn22222iiiiiiiikikikiRnCnGCRUCFV将将作为未知量对积分后的式子进行整理，得到：作为未知量对积分后的式子进行整理，得到： （10）三、点堆动力学方程的数值解法点堆动力学方程的数值解法 3、厄密特多项式方法、厄密特多项式方法12345611111112222222333333344444445555555666666 RR UVVVV VRVUVVVVRVVUVVVRVVVUVVRVVVVUVRVVVVV6 U2112222332442552662 nCCCCCC123456 FGGGGGG （11）三、点堆动力学方程的数值解法点堆动力学方程的数值解法 3、厄密特多项式方法、厄密特多项式方法223232230231112311223223232301()(32)121()()21211(32)(2)+()2121 ()(32)12()111()212(iiiiiiiiiiiiiiiiiiiiiiiiiihRa h AAAAbhhAA hhFAAAnhp AAACqhRha B