大规模离散分数阶Cable方程的数值模拟方法

埃及信息学杂志（2015）16，37开罗大学埃及信息学杂志www.elsevier.com/locate/eijwww.sciencedirect.com原创文章大规模离散分数阶Cable方程新罕布什尔Sweilam*，Hatem Moharram，N.K.阿卜杜勒·莫尼姆埃及开罗大学理学院数学系接收日期：2014年7月17日;修订日期：2014年12月17日;接受日期：2014年12月20日2015年3月29日在线发布摘要本文提出了一种分数阶Cable方程的数值模拟方法在大尺度域中。特别关注的是分数加权平均有限差分法（FWA-FDM）的分布式系统上的并行执行与显式消息传递，其中分数导数定义在Riemann-Liouville意义下。采用预条件共轭梯度法（PCG）对得到的庞大方程组进行了研究，并在其上实现了集群计算，该方法通过最小化进程间通信，满足了在Linux PC集群上实现的适用性。为了验证该方法的有效性和准确性，采用不同数目的Linux PC集群节点进行了数值测试实验.性能指标清楚地显示了在Linux PC集群上使用所提出的方法在执行时间减少和加速方面相对于在单个PC中顺序运行的好处©2015制作和主办由Elsevier B.V.代表计算机与信息学院开罗大学。 这是一篇基于CC BY-NC-ND许可证的开放获取文章（http：//creativecommons。org/licenses/by-nc-nd/4.0/）。1. 介绍从医学上已知，神经细胞是通过电化学信号处理和传递信息的电可兴奋细胞。神经细胞之间的这些信号通过突触发生，突触是与其他细胞的专门连接。大脑是一个作为神经*通讯作者。电 子 邮 件 地 址 ： nsweilam@sci.cu.edu.eg （ 新 罕 布 什 尔 州 ）Sweilam）。开罗大学计算机和信息系负责同行审查。系统大脑中神经细胞的数量因物种而异。据估计，人类大脑大约有1000亿个神经细胞和100万亿个突触。大脑的功能是对身体的其他器官进行集中控制。扩散起着关键作用大脑功能分数电缆方程起着重要的作用，模拟神经细胞中的离子的电扩散与异常亚扩散沿和跨神经细胞。PC机群系统是一种低成本的通用并行计算系统。集群计算是目前处理这些具有挑战性的问题，通常需要高计算能力的最成功的替代方案之一。著名的消息传递接口（MPI）是集群并行编程范式的少数代表。然而，在集群环境中实现它会产生一些问题，特别是http://dx.doi.org/10.1016/j.eij.2014.12.0011110-8665 © 2015由Elsevier B. V.代表开罗大学计算机与信息学院制作和主办。这是一篇基于CC BY-NC-ND许可证的开放获取文章（http://creativecommons.org/licenses/by-nc-nd/4.0）。制作和主办：Elsevier关键词加权平均有限差分法;分数阶Cable方程;预条件共轭梯度法;并行计算; LinuxPC集群工作站38新罕布什尔Sweilam等人~不¼2.x;tmj半]2JHKJJGruünwald权重，定义为wa-1k在第二步中，J不xj;t mtjMMm提出了一种基于MPI的、适用于集群基础设施的高性能计算算法及其有效的映射、调度和跨节点执行。近年来，各种生物系统都表现出反常扩散，其速率不能用扩散常数这一单一参数来表征[1]。这些生物系统中的反常扩散偏离了布朗运动的标准Fichean描述，其主要特征是其均方位移是关于时间的非线性增长，例如x2（t）t1。 作为例如，单个 颗粒 跟踪 实验揭示了蛋白质和脂质在各种细胞膜中的亚扩散（0 a1）[1<<在神经细胞粘附分子中也观察到异常亚扩散[5]。实际上，异常扩散发生在许多其他物理情况下，例如多孔介质中的流体输运[6]，以及粘弹性介质中的机械扩散波的传播[7]。由于其与布朗运动动力学的显著偏离，上述生物系统中的异常扩散不能用传统连续且y0∈x∈i在区间al[0，x]中可积，对于每一阶0a 1，n都是<01网格点为xj<$ajDx和tm1/4mDt，以及其中n是超过a的最小整数，C（.）是伽马函数。如果am N，则（1）与经典的m阶导数y（m）（x）一致。定义2. Gruünwald-Letnikov定义的 近似解 ux;tonthese 网格点uxj;tmuUm：有关分数微积分见[16，17]。在第一步中，常微分算子离散化如下[16]：函数y的a>0阶分数阶导数为@u。m1um1-um定义如下：DayxlimX1 XH wayx-hk;x≥0;2@tDt2和2¼dtujþOðDtÞ÷DtODt;6h！0 HaHKk¼0@u。u-2U-2U2xxaK.1/4天K ：umODx2j-1jj1ODx2：7哪里表示的整数部分，w是归一化的. 一个@x2xj;tmxxjDxGruünwald-Letnikov定义的一般化具体如下：的普通离散化公式的整数阶导数-提维斯在相对较弱的条件下，方法一致;如果y=x，D1-cux;tj哪里1/4d1-cumO分数阶Cable方程的集群计算39tm117-KJ不0K2¼J22.！--X1/4 Dt，N一Rr1不JKJk¼0焕光u-m-k;DtDt拉瓜权重w k的生成函数（一）（一）（二）（三）（一）（四）（五）（六）（七）（八）（八）（627我的天啊; bm; bm]，M12N-2N-1J不xj;tm不J不J不J不JJJj-1jj 1Dxw-1-b-1374X1N2N-2N-1半]d1-cummw1-c u; t-kD tMXmhi1XDtk¼0RUmNkw1-b1-kw1-br¼0M1XmDt1- c-1-c-1-c-1- cJhUm-r-2Um-rUm-ri-lNaXw1-aU m-r：18r¼0其中tm表示tm的整数部分 为了简单起见，等式（17）可以转换成矩阵形式，如下所示：选择h¼Dt。权重wk有多种选择[16，18]所以上面的公式不是唯一的。 让我们表示拉瓜2 阿萨姆邦6ddmam.0......0 0 07..X1拉克拉克A¼6. .. .. .0wð z;aÞ¼k¼0 wkz：10006丁姆河N-2阿萨姆邦N-2ðmÞ丁姆河ðmÞ7如果wz;a1-za;11h然后（9）给出第一个的向后差分公式其中bM1dN-1aN-15顺序，这被称为Gruünwald-Letnikov公式。的系数w_a_b可以通过递归公式计算u<$hum;um;. . ; um; um]; ak1 2/; j我也是。1-一个100万美元的wa;wa1：121/4; 2;.. . ; m-1; m <$1; 2;.. . ; M-1; d-1/d-1/dkkk-1 0对于c1，算子D1-c成为恒等算子，因此，等式1的一致性是成立的。（8）和（9）要求w=0 1，且对k P 1要求w= 01，这又意味着w=z; 0 1。现在，为了获得电缆方程的有限差分格式。（3），我们在方程的中间点处对该方程进行gridxj;t m：分数加权平均差分格式是Eq.（十七）、幸运的是，Eq。方程（17）是三对角方程组，可以用共轭梯度法求解。在k1的情况下，和k1/4，我们有向后欧拉分数阶求积方法和Crank-Nicholson分数阶求积方法，分别已经研究过，例如，在[19]中，但在k0处，该方案被称为全隐式。稳定性分析见[20]。ut然后，我们用前向差分公式（6）代替一阶时间导数，用三点中心公式（7）在时间tm和tm 1dtum-nkd1-bdxx um1-kd1-bdxx um1old1-a um¼ Tm;现在的目标是引入Eq。（17），以求解线性方程的三角系统，其中右手侧利用到该时间为止的计算解的所有历史。我们使用PCG方法来解决这个系统。共轭梯度法简介j t jt j tjJð14Þ1952年的算法可以在[21]中找到。该方法用于求解以下形式的线性系统其中k是权重因子，Tm是所得到的截断误差。标准差公式由下式给出：dtUm- fkd1-bdxxUm1-kd1-bdxxUm1gld1-aUm 0：150现在，通过代入由（6）、（7）和（9）给出的差分算子，我们得到Ax¼b19对于（对称，正定）和bRn：在这种情况下，线性系统的解等价于最小函数[22]：爱因斯坦：R n！ R; x！ 1 x T Ax-x T bUm1-UmDt1XmDtr¼0.Um-r-2Um-rUm-r！R2这意味着x解Ax½b当且仅当E有全局极小x处的mum：1m-1 -k×w-1-b-1Um1-r2Um1-rUm1-rj-1jj 12这种等价性是共轭梯度的基本思想算法不是以典型的方式求解线性系统nDtRr¼0Dx我们搜索函数E的最小值。设x¼x02Rbe1m带1-aUm-r1/40：160的带任意的起始向量。我们搜索E的最小值，线DtR Jr¼0g：R n！ R; a ¼ x ap：把NbBDxa<$Dt，u<$D 1-k<$Nb，在某些情况下，为了近似地获得最小值，我们使用迭代通过简化，我们可以得到以下形式：-uUm112uUm1-uUm1R;17用不同的搜索方向搜索。共轭梯度法在矩阵ðmÞ2JBDt¼j-1Jj1RJ222.2. PCG方法-k40新罕布什尔Sweilam等人j-1哪里jj1条件良好（即条件数不太大）;然而，在实际应用中，大多数矩阵是病态的分数阶Cable方程的集群计算41KKKKKK~e0K1=2¼ BCIBBRK1不不12T不不K不不k1KKrTBrKKK2KKKK KKT¼KK1=2(i.e.条件数很大），从而降低了算法的效率。为了提高算法的效率，我们采用了预处理技术来提高条件数。ak1¼pTCyk-bp CpkpTrkpTCpkrT rkTCpk一个矩阵的BER。通过使用前提条件[23]，我们可以迭代求解方程。（19）比原来的问题更快。预处理背后的想法是在一个[24]见[25]，见[26]。可以使用以下公式迭代地计算残差rk 1/2b-Cyk1等效系统因此，代替求解Ax¼b，我们求解一个相关的问题Aex~ 1/4b，其中Ae被选择为使得rk1 ←rk -ak1Cpk残差它的条件数更接近于1;换句话说，A接近于恒等式。2.3. PCG方法设Ax~<$b~是Ax~<$b~的变换系。两因为rk-ak1Cpkb-Cykapkb-Cyk1<$rk1：方程的初始残差r~0之间的关系（20）0¼b;可以通过以下方式找到：e e初始残差r这些系统通过以下关系相互关联[24]：B1=2 b-Cy的Ax1= 21= 2~1= 2AB1= 21/4~r0乘以B-1=2，从左起，我们得到b-AB1=2y1 /4Ae¼BA;x¼By;b¼b y1= 2B-1=2r~0或b-Ax0/3B-1=2r~0。 因此r01/4B-1=2r~0。推广我们得到rk/B-1=2r~k：y/B-1=2x：B称为预条件子。换一种符号，设C 1/4B1= 2AB1=2。然后，我们没有解决类似地，我们可以得到p^B-1=2p~k：Ax=b，我们必须解决以下相关问题：Cy¼B1=2b;x¼B y：ð20Þ在下面的所有计算中，我们将B-1=2~rk代入在RK。通过将r~k替换为B1=2rk，残差方程变为：B1=2rk1/4B1=2rk-~ak1Cpk或B1=2rk1最简单的预条件子是对角矩阵，对角线条目与本文中的条目相同，我们应用¼B rk-~ak1B1=2AB1=2p~k这个预条件，被称为对角预条件或雅可比预条件[25]。PCG算法具有以下两个部分，重复执行这两个部分，直到通过检查误差准则来执行PCG方法的收敛，即，残差向量的欧几里德范数应当从左边开始，将两边乘以B-1=2，我们得到：rk1<$$>rk-~ak1AB1=2p~k;或rk1<$rk-~ak1Ap2.3.2.部分（b）：计算新的搜索方向小于规定的公差。由于C是正定的，因此rk1 Crk2.3.1. 部分（a）：计算新的平均值参数ak;和残差rk1k1 ;搜索可以迭代地创建搜索方向将至pk1<$rk 1<$bk1pk和p0<$r0：搜索方向为了计算y，我们设置任意的起始向量x，计算最小值的更精确近似值，从左起两边乘以B1=2，我们得到：每次迭代B1=2p~¼B1=2~r哇~B1=2p~i：e：p¼B1=2~rb~pyk1←yk akpk：或pk1k11= 21= 2k1Þ þb~kk1p 或p1/4溴k1哇~k1kp如果我们替换yk 其中B-1=2xk：当计算k1k1k1kk1k1k1kyk≠ 1，则方程变为：B-1=2x←B-1=2x~ap~我们还找到了一个新的改进公式b~k≠1k1k k k-1=2bB1= 2r不 B1 = 2rÞrTB1= 2B 1= 2rk1两边都乘以B从左边开始，yk的值为1~k≤ 1¼k1Tk1k进一步转化为B1=2rrTB1= 2B1= 2rkk k kRT Brk1x← xB1=2p~;或或者b~¼k 1：xk1←xk~akp~k我们可以找到一个k的解析公式。对于固定的y k 和p k，2.4. PCG的迭代公式如下所菲乌伊阿罗帕¼ ðyþapÞ Cðyap1/4apCpap Cy-ap b.. .2kkkk kk当导数为零时，f相对于a的最小值pT Cy所以普纳pk-pkb¼0rk<$1<$ rk-ak<$1Cpk：xk1←xka~kp~k：pk1¼Brk1b~k1pk：B~k1¼RT 第一章：k 1//第k+1//返回//搜索方向：//在步骤krTBrkKak1¼rT rkkpTCpk：//steplength：K00K¼r~0;将B1的两侧相乘=2b-B1=2我们选择B对称正定矩阵Kk1KK：K42新罕布什尔Sweilam等人J“u2222半个小时后H1-BJ8xxjTm¼ O hpJTJH1-BRJJ21JJ22.5. [26]第26话我的世界并行算法与串行算法相同，但并行实现的一些例外是但dxxum1-r<$uxxDx12 uxxxxDt2xxtDx12 u xxxxt.. .#þðDtÞu我... ;1. 从input_le中读取数据并将其划分为cessors（使用MPI_Bcast和MPI_Scatter）。2. 在每个处理器在本地计算内积之后，需要跨所有处理器进行求和归约（使用和8xxttDx12Dt“DxMPI_Allreduce）。3.向量矩阵积需要收集所有局部将向量的部分转换为单个向量，然后每个处理器执行乘法（使用MPI_Allgather）。dxxum1-r<$uxxþðDtÞuuxxxx2我... ;uxxt12 u xxxxt.. .哪里 的 部分 衍生物 是 评价 在 的点下面是PCG的并行代码片段，求解形式为Ax 1/4b的线性方程组。xj;tM·K 你好。 将这些表达式插入Eq. （21）和. 1ΣJDxSSTm¼Ohp- -kDtD1-buxxt-2DtD1-buxxxx12-D1-bu -1-ku0.. . ;8sxxttH1-Bm1xxj与S¼tm，即. 1Σþ1 w1-bdu0H 1-Bm1xxj其中阶O的项 b p> 2的DtaDxb hp未包括在内。H3. 稳定性分析和截断误差定理3.1.方程的截断误差（3）由4. 性能分析和PC集群描述并行计算的主要目的是使计算运行得更快。更快显然意味着为了了解什么是可能的，以及我们可以期望实现什么，我们使用几个指标来衡量并行性能，每个指标都有自己的优点和缺点。第一个是执行时间。执行时间TP是指一个并行程序在P个处理机上的净执行时间，不包括初始操作系统、I/O等费用。第二个是加速。加速比SP，是顺序程序的执行时间除以计算相同结果的并行程序的执行时间具体地，S P= T S/T P，其中T S是顺序表达式。Tm¼Ohp1日本语2-kODtODt2ODx21w-1-b-1-dm1您的位置：TP是P个处理器上的并行执行时间。第三是效率。效率E P是加速的标准化测量：E P=S P/P。理想情况下，加速比应按比例缩放-证据从方程给出的截断误差的定义。（14），得到：早期的P，这意味着效率应该有一个恒定的值1。当然，由于性能损失的各种来源，效率通常低于1，并且它会降低Tm¼d um- fkd1-bd1-bdμm≤ 1g;随着处理器数量的增加， 效率jt j即ttxx jMttxx j大于1表示超线性加速[26]。并行PCG已经在开罗大学科学学院的计算机集群上实现[27]，Tm¼dum-1Xw1-bh1-kdum1-rkdum-ri17个工作站，一个主站和八个从站，类似于图1.一、集群中的每个从机都有一个英特尔®微处理器-1-kw1-bdu0：21Core（TM）2 i7- 2600@3，40 GHz和8 Gb DDR3 RAM。的H1-Bm1xxj主处理器具有英特尔（R）酷睿（TM）2四核CPU微处理器Q6700 2.66 GHz和4 Gb DDR2 RAM。 奴隶们r_local =b_local//初始残差rho= Allreduce（r_local以后用对于k=1：itermax如果k=1p_local=r_localelseβ= ρ/oldrho//初始方向。//在步骤k处的改进，即Bp_local=r_local+ beta*p_local//搜索方向结束p=Gather（p_local）//所有搜索方向的总和处理器v_local=C_local *p//要使用的中间项alpha= rho/Allreduce（p_local所有处理器x_local=x_local+ alpha *p_localr_local=r_local- alpha *v_localoldrho= rho//中间项需要用作term//从上一步开始。//vector// new residualrho= Allreduce（r_local在当前步骤端2222xxtt考虑到Eqs。（3）和（9），我们可以得到：-kOD tODt2 OD x2r¼0xxxx分数阶Cable方程的集群计算43ð Þ¼ ð Þ400040不不不tt本例的数值解见42N.H. Sweilam等人图3Dx¼1的近似解和Dt1，在a 1/40： 2和b=0.7。4000 40图1 PCG算法的数据依赖性.p2 tb1ta1 Σfx;t2t2b2a中国人通过星型局域网（1 Gbps）连接到主机。网络由交换机、每个工作站中的网卡和相应的Cat.6 UTP网络导线(see图 2）。集群的操作系统是Linux，其中包括用于控制并行应用程序执行的工具。选择C++作为编程语言是因为它易于管理 典 型 格 模 型 数 据 的 大 数 组 ， 并 结 合 消 息 传 递 接 口（MPI）库。5. 实验结果在下面的部分中，我们选择两个例子，我们有他们的精确解，所以我们可以检查我们的结果的正确性。在图3中，Dx¼1;Dt¼1 a：1; b： 2;b： 3显示在三维图中其中边界条件u= 0;t=1;t=0，初始条件u=x;0 =0。Eq的精确解（19）是u x;tt2 sin px.用FWA-FDM方法给出的分数Cable方程（19）的精确解的性质如图所示。 四、实施例2.考虑以下分数阶Cable方程的初边值问题：utx;t1-buxxx;t-0：5D1-aux;t;0x 10;0t6T;<<<其中u= 0;t=u= 10;t= u =0且u=x;0 = u =10d=x-5，其中d=x= d是狄拉克δ函数。显示分数电缆方程的模拟在3-D。而在图3中，我们将精确解与数值解进行了比较。例1.考虑分数阶Cable方程图 八比十二在图8中，大规模问题的近似解以3-D形式示出，作为该示例的Cable方程在3-D中的行为的说明图。在图9中，不同的数值解在不同的值，测试和计算，以说明电缆等式的行为，ux;tD1-bux;t对于这些价值观。总之，图 5-7和10-12显示了平行的扩展。在有限域0x 1上，其中06t6T，0a;b 1和以下源项：<<<<对于不同的问题规模（N=4000，6400，8000）和增加的过程数（p=1，2，图2 PC集群描述。xx分数阶Cable方程的集群计算435030图7缩放 效率 为 不同 问题 尺寸 在a1/4：2;b1/4： 7;以及不同数量的处理。图4 19的精确解和数值解在k<$0时的行为，a<$0：2，b<$0： 7，Dx<$1，Dt<$1，T1/2。100 40图8Dt¼1、Dx¼1的近似解 得双曲余切值.图5不同问题大小的缩放执行时间a 1/40： 2和b 1/40： 7以及不同数量的处理。a 1/4 0：2和b =0.7。40 4000图6 不同问题大小的按比例加速，a1/4：2和b1/4：7以及不同数量的处理。.. . ，8，16）。图图5和图10表示相对于进程数的并行执行时间。可以观察到，对于大的N（N=6400，8000），并 行 执 行 时 间 Tp 随 着 p 减 小 ， 而 对 于 小 的 问 题 大 小（N=4000），对于p=2，4，8，在实施例1中，p=1，2，8，16，在实施例2中。两个示例之间处理器数量的差异是由于所使用的集群的内存限制。图图5和图10示出了使用增加数量的过程来解决不同尺寸的问题的实验加速曲线。正如预期的那样，对于给定数量的进程，加速比随着问题大小的增加而此外，对于A图9精确解和近似解，其中k为1/4，a：1/40：5，b：1/4 0：5，D：1/4 x 1/4，D：1/4 t 1/4。图10不同问题大小的缩放执行时间a1/4：2和b1/4：7以及不同数量的处理。在给定问题大小的情况下，加速不会随着进程的数量而继续增加，而是趋于饱和。44N.H. Sweilam等人图 11缩放 加速比 为 不同 问题 尺寸.a1/4：2和b1/4：7以及不同数量的处理。图 12缩放 效率 为 不同 问题 尺寸.a1/4：2和b1/4：7以及不同数量的处理。图图6和图11示出了使用增加数量的过程来解决不同尺寸的问题的效率曲线。同样清楚的是，效率往往随着过程的数量而下降。6. 结论本文研究了基于消息传递接口（MPI）的机群上的分数阶加权平均有限差分算法FWA-FDM。由于神经细胞数量庞大，需要通过Cable方程进行模拟，因此集群计算在这一领域是必不可少的，并可以在许多生物学应用中提供帮助引用[1] [10]黄文辉，李文辉，李文辉.多光子荧光光漂白恢复法测量溶液中分子扩散。J Biophys1999;77：2837-49.[2] [10]张晓刚，张晓刚. 细胞表面的限制扩散或不动部分：一种新的解释。J Biophys1996;70：2767-73.[3] 戈什河人皮肤成纤维细胞表面单个低密度脂蛋白受体分子的迁移和聚集。博士论文Ithaca，NY：Cornell University; 1991.[4] 作者声明：R.自动检测和跟踪单个和成簇细胞表面低密度脂蛋白受体分子。J Biophys1994;66：1301-18.[5] Simson R，Yang B，Moore S，Doherty P，Walsh F，JacobsonK. 质膜亚微米尺度上的结构镶嵌现象。J Biophys1998;74：297-308.[6] [10]黄文辉，李文辉，李文辉.多光子荧光光漂白恢复法测量溶液中分子扩散。J Biophys1999;77：2837-49.[7] 迈纳尔迪湾粘弹性固体中的分数扩散波。非线性波固体1995：93[8] Henry BI，Langlands TAM，Wearne SL.多刺神经元树突的分数电缆模型。物理学评论快报2008;100（12）：128103。[9] 波德鲁布尼岛分数阶微分方程北京：人民出版社，1999.[10] 刘芳，安武，特纳，庄平。时间分数阶对流扩散方程。ApplMath Comput 2003，p. 233-246。[11] 刘F，杨Q，特纳I.分数阶Cable方程两种新隐式数值方法的稳定性和收敛性。JComput Nonlinear Dyn 2011;6（1）. 文章ID01109.p. 第七章[12] Quintana-Murillo J，Yuste SB.分数阶电缆方程的显式数值方法。Int J Juberren Equat 2011，p. 57比69[13] 作者：Langlands TAM，Henry B，Wearne S.解分数电缆方程 ： 有 限 情 况 ; 2005 。 。[14] Langlands TAM，Henryand BI，Wearne SL.神经细胞异常电扩散的分数电缆方程模型：有限域解。J Math Biol 2009;59（6），p. 761- 808[15] Rall W.核心导体理论和神经元的电缆特性。在：Poeter R，编辑.生理学手册：神经系统，第1卷。马里兰州贝塞斯达：美国生理学会，1977年。第39-97页[第3章]。[16] 卢比奇角离散分数阶微积分。SIAM J Math Anal1986;17，p.704- 719[17] Morton KW，Mayer DF.偏微分方程的数值解。剑桥：剑桥大学出版社，1994.[18] Sweilam NH，Khader MM，Nagy AM.双边空间分数阶波动方程的有限差分数值解J Comput Appl Math 2011;235，p. 2832-2841。[19] Palestine，Cuesta EA. Banach空间中一个带记忆积分微分方程的数值方法：定性性质。SIAM J Numer Anal 2003;41，p.1232-1241。[20] Sweilam NH，Khader MM，Adel MM.关于神经元动力学建模的基本方程。2013年。[21] Hestenes MR，Stiefel E.解线性方程组的共轭梯度法。J ResNat Bur Stand1952;49：409-36.[22] 小休丘克 介绍了共轭梯度法，没有痛苦的痛苦。匹兹堡：卡内基梅隆大学计算机科学学院; 1994年。[23] 本 齐 ·M 大 型 线 性 系 统 的 预 处 理 技 术 综 述 。 J ComputPhys2002;182：418-77.[24] O’Leary Dianne P. Notes on some methods for solving linearsystems;[25] 萨阿德·Y稀疏线性方程组的迭代法。伦敦：国际汤姆森出版社; 一九九六年。[26] Thaoma R，Gudula R.多核和集群系统的并行编程。BerlinHeidelberg：Springer-Verlag;2010.[27] Sweilam NH，Moharram HM，Sameh Ahmed.基于MPI集群的二维Poisson方程的并行迭代有限差分算法在：第八届国际会议信息流与系统（INFOS 2012），IEEE Explorer; 2012年。