本发明涉及核工程领域，特别涉及一种核反应堆物理燃耗方程的计算方法。
背景技术：
燃耗计算是为了计算核燃料的组成成分，这对反应堆的运行、冷却以及核电站的放射性防护有着重大的意义。该计算对计算速度和精度有着比较高的要求。现阶段，计算燃耗主要有两种方法：第一种是通过解燃耗链的方法，比如tta方法，这种方法的计算精度比较高，但是效率比较低，适合单个核素的求解。另一种方法是通过矩阵的形式求解燃耗方程，用矩阵指数法来做快速求解，这种方法的计算效率比较高，但是数值结果受到数学方法上的影响比较大。
近几年来，矩阵指数法的研究逐步深入[1]，效率和精度上都有了比较大的提高。比较有代表意义的是泰勒展开法、帕德有理展开法和切比雪夫有理展开法。泰勒展开和帕德展开是基于指数函数在原点附近的展开计算的，运用在矩阵上，主要使用与范数比较小的矩阵。当矩阵的范数增大的时候，计算精度会下降，甚至出现计算错误。采用缩放指数的方法能在一定程度上解决问题。切比雪夫有理展开方法是现阶段一个比较合适的算法，计算精度满足工程实际的需求，但是其计算速度不够快。
克雷洛夫子空间法是通过将一个大型矩阵投影到一个小的子空间内，然后在这个小空间上进行求解以得到加速算法的目的。由于燃耗计算中，燃耗矩阵的刚性比较大，是一个接近奇异的矩阵，一般的克雷洛夫子空间法在精度上并不能满足要求。
技术实现要素：
本发明的目的在于克服现有技术的缺点与不足，提供一种克雷洛夫子空间加速求解燃耗方程的数值计算方法，在采用有理展开的矩阵指数法的基础上，通过耦合新型的克雷洛夫子空间法—广义最小残差法，在保证一定的精度要求的前提上，提高了燃耗求解速度。
本发明的目的通过以下的技术方案实现：一种克雷洛夫子空间加速求解燃耗方程的数值计算方法，包括：
s1、矩阵指数有理展开：
指数函数在一定的区间内可以用有理展开的方法进行近似计算，其形式为：
这里v是有理展开式的分子分母的阶数；
采用部分分式写法可以改写成：
将其运用在燃耗计算上，可化成：
式中：n0为初始核素组成向量；v为有理展开式的阶数；ξj为有理展开式中的极点；τj为有理展开式中各极点对应的留数；τ0为一个固定常数；a为燃耗矩阵；i为单位矩阵；
s2、利用克雷洛夫子空间法—广义残差法加速：
用广义残差法对公式(1)进行加速求解，其流程如下所示：
s2-1输入一个n阶的燃耗矩阵a，一个n阶的单位矩阵i，一个初始核素组成向量n0，一个时间步长t；任意选取一个初始迭代值x0并且计算b＝ta-ξji和r0＝n0-bx0；
s2-2设置v1＝r0/||r0||2，β＝||n0||2；
s2-3循环：forj＝1,2,…,m,…,
hi,j＝(bvj,vi),i＝1,2,…,j，
其中(bvj,vi)代表求其内积；
在上面的每一次循环中都做残差计算：
建立汉斯博格矩阵：
解最小二乘法问题：
建立相应的：yj和||rj||2＝||vj(βe1-hjyj)||2，
当残差||rj||2小于设定阈值时，退出循环；
s2-4计算获得单个解：
xm＝x0+vmym
s2-5把v个线性系统的解组合计算，可以得到最终核素组成向量n(t)：
优选的，在建立克雷洛夫子空间的过程中，对矩阵的对角线做偏移处理不会改变其子空间中的向量，即：km(a,v)＝km(a-εji,v)，在这种情况下，有理展开式子的v个线性系统中，只需要做一次子空间的建立，并在这个基础上做各个线性方程的残差迭代计算。
优选的，在进行有理展开时，进行ilu分解预处理：
l-1(a-ξji)r-1u＝l-1n0；(2)
yj＝r-1u；(3)
将要求解的矩阵进行ilu分解，然后将得到的矩阵l和u的逆分别乘到矩阵两边，方程的另一端左乘l的逆，然后通过广义残差法迭代得到的结果再左乘r的逆，得到最后的结果。
优选的，在求解燃耗的时候引入重启动技术，在进行了若干步迭代之后若未收敛，则将得到的解和残差作为初始值，重新开始建立克雷洛夫子空间并进行迭代，具体如下：
(1)任意选取初始值x0和内迭代数k，并且计算b＝ta-ξji和r0＝n0-bx0；
(2)令v1＝r0/||r0||2,β＝||n0||2；
(3)循环：forj＝1,2,…,k
按照步骤s2广义残差法的流程计算：令||rj||2达到最小的yj值；
(4)构建近似解：
xk＝x0+vkyk
(5)重启动：
计算rk＝n0-bxk；如果该值满足计算需求则结束迭代；否则令x0＝xk,v1＝rk/||rk||2，然后重复步骤(3),(4),(5)。
本发明与现有技术相比，具有如下优点和有益效果：
在满足相应的计算精度需求的前提下，通过建立一个子空间并在该空间内求解，与整个n维空间上相比，大大减少了计算量，极大地提高了计算效率。
附图说明
图1是gs/gmres迭代需要的时间对比图。
图2是gs/gmres的迭代步数对比图。
具体实施方式
下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。
实施例1
在中子通量保持不变的一个燃耗步长内，采用新型的克雷洛夫子空间来加速求解燃耗方程。为了能在求解燃耗的过程中加速矩阵指数有理展开算法，将新型克雷洛夫子空间法—广义残差法[2]耦合到原有的矩阵指数有理展开法，在保证一定的精度要求的前提上，提高了燃耗求解速度。主要是将一个n维的燃耗矩阵和初始核素浓度投影到一个m维子空间上，然后在这个子空间中寻找最优解。由于子空间的维度m比燃耗矩阵自身的维度n小很多，所以在求解的时候计算量大大减少，计算速度因此得到较大的提升。对该算法做了进一步优化，包括矩阵偏移技术，预处理技术和重启动技术的应用，使得该算法更加适合燃耗方程的快速求解运算。
具体采用的技术方案是：
1、矩阵指数有理展开方法
指数函数在一定的区间内可以用有理展开的方法进行近似计算，其形式为：
这里v是有理展开式的分子分母的阶数；
采用部分分式写法可以改写成：
将其运用在燃耗计算上，主要可化成：
式中：
n0为初始核素组成向量；
v为有理展开式的阶数；
ξj为有理展开式中的极点；
τj为有理展开式中各极点对应的留数；
τ0为一个固定常数；
a为燃耗矩阵；
i为单位矩阵；
2、新克雷洛夫子空间法—广义残差法加速
用广义残差法对公式(1)进行加速求解。其算法流程如下所示：
输入：一个n阶的燃耗矩阵a，
一个n阶的单位矩阵i，
一个初始核素组成向量n0，
一个时间步长t。
输出：克雷洛夫子空间维度m，
一个线性系统迭代求解的结果xm。
具体过程如下：
(1)任意选取一个初始迭代值x0并且计算b＝ta-ξji和r0＝n0-bx0；
(2)设置v1＝r0/||r0||2，β＝||n0||2；
(3)循环：forj＝1,2,…,m,…,
hi,j＝(bvj,vi),i＝1,2,…,j，
其中(bvj,vi)代表求其内积；
在上面的每一次循环中都做残差计算；
建立汉斯博格矩阵：
解最小二乘法问题：
建立相应的:yj和||rj||2＝||vj(βe1-hjyj)||2,
yj是一个n*1向量，使||rj||2达到最小值时有对应向量值yj，当残差||rj||2足够小时，退出循环。
(4)计算获得单个解：
xm＝x0+vmym
(5)把v个线性系统的解组合计算，可以得到最终核素组成向量n(t)：
3、新克雷洛夫子空间法的进一步优化处理
3.1做矩阵的偏移处理
由于矩阵的偏移特性[3]在建立克雷洛夫子空间的过程中，对矩阵的对角线做偏移处理不会改变其子空间中的向量，既是：km(a,v)＝km(a-εji,v)。在这种情况下，有理展开式子的v个线性系统中，只需要做一次子空间的建立，并在这个基础上做各个线性方程的残差迭代计算。这样充分节省了计算时间，提高了计算效率。特别适用于燃耗求解。
3.2做矩阵的预处理
由于燃耗计算中，不同核素的半衰期差异很大，导致燃耗矩阵的刚性比较大，自身接近一个奇异矩阵，在进行有理展开时，不利于线性方程的求解。所以需要进行预处理使得方程容易求解。做一个ilu分解可以很好的解决这个问题。
l-1(a-ξji)r-1u＝l-1n0；(2)
yj＝r-1u；(3)
将要求解的矩阵进行ilu分解，然后并将得到的矩阵l和u的逆分别乘到矩阵两边，方程的另一端左乘l的逆。然后通过广义残差法迭代得到的结果在左乘r的逆，得到最后的结果。这样做的收益是，所求矩阵在两边乘上l和u的逆的时候，矩阵刚性正好得到了改善，不再是一个奇异矩阵，因此可以很快的进行迭代求解。
3.3重启动处理
当线性方程迭代次数比较多时，克雷洛夫子空间的维度增加会造成计算量的不断增大。为了解决这个问题，本发明在求解燃耗的时候引入了重启动技术。其主要思想就是在进行了若干步迭代之后若未收敛，则将得到的解和残差作为初始值，重新开始建立克雷洛夫子空间并进行迭代。其算法如下：
(1)任意选取初始值x0和内迭代数k，并且计算b＝ta-ξji和r0＝n0-bx0；
(2)令v1＝r0/||r0||2,β＝||n0||2；
(3)循环：forj＝1,2,…,k
按照前文广义残差法的算法流程计算：令||rj||2达到最小的yj值；
(4)构建近似解：
xk＝x0+vkyk
(5)重启动:
计算rk＝n0-bxk；如果该值满足计算需求则结束迭代；
否则令x0＝xk,v1＝rk/||rk||2，然后重复步骤(3),(4),(5)。
对一个1696阶的燃耗矩阵进行计算分析，主要以切比雪夫有理展开方法为基础，并与用高斯赛德尔迭代法进行全空间迭代求解的算法进行比较分析。图1、图2计算结果中容易看出通过用新型克雷洛夫子空间法加速求解的算法，计算效率提高了接近5倍。由于在用迭代法解线性方程时，规定了相同的收敛判据，因此认为全空间和子空间的算法在计算精度基本一致。
上述的实例是以一个范数为1015的1696阶燃耗矩阵为代表进行说明的，但本发明提出的加速求解方法理论上对于燃耗矩阵均适用。
[1]moler,c.andc.vanloan,nineteendubiouswaystocomputetheexponentialofamatrix,twenty-fiveyearslater.siamreview,2003.45(piis00361445024180101):p.3-49.
[2]saad,y.andm.h.schultz,gmres:ageneralizedminimalresidualalgorithmforsolvingnonsymmetriclinearsystems.1986.
[3]botchev,m.a.,krylovsubspaceexponentialtimedomainsolutionofmaxwell'sequationsinphotoniccrystalmodeling.journalofcomputationalandappliedmathematics,2016.293:p.20-34.
上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。