椭圆方程五点格式的迭代法与快速算法的比较.pdf

丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌丌保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报保山学院学报椭圆方程五点格式的迭代法与快速算法的比较焦岑孙唯唯聂家升（苏州大学应用技术学院 通识教育学院，江苏 苏州 215325）摘要 主要讨论在椭圆方程五点格式的问题中，分别使用Gauss-Seidel迭代法与快速Poisson算法对其求解，并对二者求解该线性方程组的速度进行比较。在系数矩阵是稀疏的大型线性方程组中，迭代法是一个很好的求解该类型的算法，主要是因为给定一个初始向量，通过一定的迭代公式，可以求得之后任意一次迭代的结果，且运算简便，但是，对于迭代法所求得的近似解是否收敛于精确解，并且，在线性方程组有快速算法的情况下，迭代法是否还能在求解方程组中占优势，还需进一步比较。通过比较不同的系数、不同的步长以及不同的误差要求，来判断Gauss-Seidel迭代法与快速Poisson算法的优劣。关键词 椭圆方程五点格式；Gauss-Seidel迭代法；快速Poisson算法中图分类号 O13文献标识码 Adoi:10.3969/j.issn.1674-9340.2024.02.009文章编号 1674-9340（2024）02-0056-10收稿日期：2023-09-19基金项目：2022年江苏省哲社一般项目“基于OBE-CDIO教育理念的新工科大学数学课程体系重构与教学内容改革研究”（项目编号：2022SJYB1538）。第一作者简介：焦岑（1996），女，汉族，江苏扬州人，硕士，讲师，研究方向为应用统计分析研究。1 背景介绍本文从椭圆方程五点格式问题出发，对其分别使用Gauss-Seidel迭代法和快速Poisson算法进行求解，并将迭代法与快速算法的速度进行比较，深入探究迭代法是否更方便，从而为以后的解题拓展新的思路。由于此问题一方面涉及偏微分方程，而另一方面又涉及数值计算，大多数文献都只是考虑了一半的内容，因此，探究该问题就显得特别地重要，主要是因为当所需求解的线性方程组的系数矩阵是大型稀疏矩阵时，迭代法是求其解的一个重要方法1-2，6-7。此外，比较Gauss-Seidel迭代法和快速Poisson算法对于求解大型线性方程组的快慢，了解迭代法的优势与不足，可以为今后考虑问题提供一个指路灯。综上所述，探讨迭代法对于解决问题有着切实可行的意义。而用迭代法与快速算法进行比较，了解迭代法的特点，也为改进迭代法作了铺垫。大多数文献在关于椭圆方程五点格式的问题上，主要都是从概念的意义上去理解，很难求得精确解，而且没有直观的数据显示。因而讨论椭圆方程五点格式的Gauss-Seidel迭代法和快速Poisson算法的比较，就显得尤为重要。一方面，数值解比概念上的描述显得更为直观有效，更有说服力；另一方面，椭圆问题的实现，为今后探讨其他偏微分方程问题，提供了一个借鉴，即使是大型线性矩阵，也能从数值解上判断优异，更有利于解决一些实际的问题。从椭圆方程五点格式入手，主要通过使用Gauss-Seidel迭代法以及使用快速Poisson算法计算不同条件下的误差。通过修改不同的步长M2、M1，不同的系数来探究Gauss-Seidel迭代法与快速Poisson算法的优劣。2 预备知识2.1 椭圆方程的介绍本文中，考虑的是二维Poisson方程在Dirichlet边值条件下的问题焦岑，孙唯唯，聂家升：椭圆方程五点格式的迭代法与快速算法的比较u-u=f()x,y,()x,y （2.1.1）u=()x,y,()x,y （2.1.2）其中，u=2ux2+2uy2。为简单起见，只考虑为正方形区域=()x,y|0 x 2,0 y 2。2.2 差分格式的建立首先，把横轴上的区间0,2进行M2等分，记h2=2/M2为x方向的步长，且有xi=0+ih2，0i M2，然后，把纵轴上的区间0,2进行M1等分，记h1=2/M1为y方向的步长，且有yj=0+jh1，0 j M1。接着，用两簇等距的平行线3x=xi，0i M2，y=yj，0 j M1，将区域划分成M2 M1个小矩形，以两簇直线的交点为结点（xi,yj）。记h=(xi,yj)|0 i M2,0 j M1为属于的结点，其中，h=(xi,yj)|1 i M2-1,1 j M1-1，被称为h的内结点，而称位于上的结点为边界结点，且有h=hh。显然，有h=h h。为方便来看，记=(i,j)|(xi,yj)h，=(i,j)|(xi,yj)h记Sh=v|v=vij|0 i M2,0 j M1为h上的网格函数，设v=vij|0 i M2,0 j M1 Sh，引进如下记号：Dxvij=1M2(vi+1,j-vij)，Dx vij=1M2(vi,j-vi-1,j)，Dyvij=1M1(vi,j+1-vij)，Dy vij=1M1(vi,j-vi,j-1)，2xvij=1M2(Dxvij-Dx vij)，2y=1M1(Dyvij-Dy vij)，v=max|vij|，称 v为v的无穷范数。在结点处考虑边值问题u-u=f()x,y,(x,y)，u=()x,y,(x,y)则有u()xi,yj-2ux2()xi,yj+2uy2()xi,yj=f()xi,yj,(i,j)（2.2.1）u()xi,yj=()xi,yj,()i,j （2.2.2）定义h上的网格函数U=Uij|0 i M2,0 j M1其中Uij=u()xi,yj,0 i M2,0 j M1引理1：如果g(x)C4c-h,c+h，则有g()c=1h2g()c+h-2g()c+g(c-h)-h212g()4()4,c-h 4 c+h由引理1，有2ux2()xi,yj=1h22u()xi-1,yj-2u()xi,yj+u()xi+1,yj-h22124u()ij,yjx4=2xUij-h22124u()ij,yjx4,xi-1 ij xi+12uy2()xi,yj=1h21u()xi,yj-1-2u()xi,yj+u()xi,yj+1-h21124u()xi,ijy4=2yUij-h21124u()xi,ijy4,yi-1 ij 0，则由（2.3.2）知，存在(i0,j0)使得|ui0,j0=M，且|ui0-1,j0，|ui0+1,j0，|ui0,j0-1，|ui0,j0+1中至少有一个小于M。考虑（1.3.1）中()i,j=()i0,j0的等式，有()+2h22+2h21ui0,j0=1h22()ui0-1,j0+ui0+1,j0+1h21()ui0,j0-1+ui0,j0+1将上式两边取绝对值，可得()+2h22+2h21M 1h22()|ui0-1,j0+|ui0+1,j0+1h21()|ui0,j0-1+|ui0,j0+1 0，所以上式与假设M 0矛盾。故M=0。因而差分格式（2.2.6），（2.2.7）是唯一可解的。2.4 差分格式的求解差分格式（2.2.6），（2.2.7）是以uij|1 i M2-1,1 j M1-1为未知量的线性方程组。（2.2.6）可以改写为-58焦岑，孙唯唯，聂家升：椭圆方程五点格式的迭代法与快速算法的比较1h22ui-1,j-1h21ui,j-1+()+2(1h22+1h21)ui,j-1h22ui+1,j-1h21ui,j+1=f()xi,yj1 i M2-1,1 j M1-1（2.4.1）记uj=u1ju2juM2-1,j,0 j M1利用（2.2.7）可将（2.4.1）写为Duj-1+Cuj+Duj+1=fj,1 j M1-1（1.4.2）其中C=+2()1h22+1h12-1h22-1h22+2()1h22+1h12-1h22-1h22+2()1h22+1h12-1h22-1h22+2()1h22+1h12D=-1h12-1h12-1h12-1h12,fj=f()x1,yj+1h22()x0,yjf()x2,yjf()xm-2,yjf()xm-1,yj+1h22()xm-2,yj可进一步可以写为CDDCDDCDDCu1u2uM1-2uM1-1=f1-Du0f2fM1-2fM1-1-DuM1由上式可以看出，它是一个大型的线性方程组，系数矩阵是一个三对角矩阵，且矩阵的每一行至多有5个非零元素。在数学上，称这种系数矩阵为大型的稀疏矩阵，因为该类矩阵中大部分的元素都是0。一般情况下，会使用Jacobi迭代法、Gauss-Seidel迭代法或者超松弛迭代法对其求解4-5。显然，可以运用已学的知识证明上式的系数矩阵是对称正定的。3 Gauss-Seidel求解及快速算法的求解3.1 Gauss-Seidel迭代法3.1.1 Gauss-Seidel迭代法的求解首先，考虑一下下面的方程组：Ax=b（3.1.1）-59第 43 卷第 2 期保山学院学报2024 年 4 月其中，A=a11a22ann，A=a11a12a1na21a22a2nan1an2ann,x=x1x2xn,b=b1b2bn然后，把A拆分成三个部分：A=a11a22ann-0-a210-an1-an20-0-a12-a1n0-a2n0 D-L-U通过选取所需的M为系数矩阵A的下三角部分，即选取M=D-L（下三角矩阵），A=M-N，就得到了使用Gauss-Seidel迭代法求解Ax=b的途径：x()0,初始向量,x(k+1)=Bx(k)+f,k=0,1,（3.1.2）其中B=I-()D-L-1A=()D-L-1()D-L-A=()D-L-1U G,f=()D-L-1b称G=()D-L-1U为求解Ax=b的Gauss-Seidel迭代法的迭代矩阵6-7。下面，给出的是使用Gauss-Seidel迭代法求解方程组所得的未知向量的具体某一个分量的计算公式：x(k)=(x()k1,x()ki,x()kn)T从式（2.2.2）可以得到：(D-L)x()k+1=Ux()k+b或Dx()k+1=Lx()k+1+Ux()k+b即：aiix(k+1)i=bi-j=1i-1aijx()k+1j-j=i+1naijx()kj,i=1,2,n于是，使用Gauss-Seidel迭代法来求解线性方程组Ax=b，具体的迭代法的计算公式如下：x(0)=(x()01,x()0n)T,初始向量,x(k+1)i=()bi-j=1i-1aijx()k+1j-j=i+1naijx()kjaiii=1,2,n;k=0,1,.(3.1.3)通过式（3.1.3）可以了解到，当计算所求向量x(k+1)的第i个分量x(k+1)i时，Gauss-Seidel迭代法会利用之前已经计算出的前i-1最新分量x()k+1j(j=1,2,i-1)，并将其带入迭代公式中进行计算，且Gauss-Seidel迭代法每迭代一次只需要计算一次矩阵与向量的乘法。对于要求解的椭圆方程1h22ui-1,j-1h21ui,j-1+()+2(1h22+1h21)ui,j-1h22ui+1,j-1h21ui,j+1=f()xi,yj1 i M2-1,1 j M2-1(3.1.4)使用Gauss-Seidel迭代公式（3.1.3），可以得到：ui,j()k+1=()f()xi,yj+1h22ui-1,j()k+1+1h21ui,j-1()k+1+1h22ui+1,j()k+1h21ui,j+1()k()+2()1h22+1h211 i M2-1,1 j M1-13.1.2 Gauss-Seidel迭代法的收敛性定理2：解线性方程组Ax=b的Gauss-Seidel迭代法收敛的充分条件是A正定。证明：因为所求的椭圆方程中的系数矩阵A是三对角矩阵，并且A严格的对角占优，故Gauss-Seidel迭代法收敛4，8。-60焦岑，孙唯唯，聂家升：椭圆方程五点格式的迭代法与快速算法的比较3.2 快速Poisson算法引理2：若A为三对角矩阵，即：A=bcabacbacb则A的特征值j=b+2ac cos()jn+1，特征向量 xj=()ac12sin()jn+1()ac22sin()2jn+1()ac32sin()3jn+1()acn2sin()njn+1考虑Dirichlet边值问题的标准中心差分格式u+u=f()x,y,(x,y)（3.2.1）记=0,1 0,1，h=1(M+1)是等距网格步长，记=h2，则数值格式等价于()ui-1,j+ui+1,j+()12-2uij+()ui,j-1+ui,j+1+()12-2uij=f()x,y为书写简单，记=2 h2，则格式等价于2()ui-1,j+ui+1,j+()12-uij+2()ui,j-1+ui,j+1+()12-uij=f()x,y即：()ui-1,j+ui+1,j+()1-2uij+()ui,j-1+ui,j+1+()1-2uij=2f()x,y（3.2.2）记T=1-2-1-2-1-2-1-2-，u=u1ju2junj,f=2f1j2f2j2fnj,j=1,2,n则式（3.2.2）的矩阵形式可写为：TU+UT=F(3.2.3)由引理2可知，Tsj=Ejsj,j=1,2,n其中，Ej=()1-2+2cos()jn+1，sj=sin()jn+1,sin()2jn+1,sin()njn+1T记S=s1,s2,sn,E=E1E2EnS=ST,TS=SE,S*S=12hI=n+12I注：S=ST，但是SE ES，SE=(ES)T。由T=TT，则ST=STTT=(TS)T=(SE)T=ETST=ES，因而由（2.3.3）可得-61第 43 卷第 2 期保山学院学报2024 年 4 月等式左右两边同时乘上SSTUS+SUTS=SFS因为ST=ES,TS=SE，故有ESUS+SUSE=SFS令V=SUS，从而有EV+VE=SFS =F因为E是对角阵，因而V可解Vij=Fij(Ei+Ej)进而，对V=SUS，等式两边同时乘上S，则有S2US2=SVS，因为S*S=12hI，故而有14h2U=SVS，即U=4h2*SVS。此时，U即为椭圆方程标准中心差分格式的所求解。具体的算法（一个简单的快速Poisson算法）如下：1、h=1m+1,F=()f()jh,khmj,k=1,S=()sin()jkhmj,k=1,=()sin2()jh2mj=12、G=()gj,k=SFS3、U=()uj,kmj,k=1，whereuj,k=h4gj,k()j+k4、V=SUS输出是在复杂度为O()n3 2运算中计算的平方上的离散泊松方程的精确解，而且存储只需要几个m m的矩阵。4 数值试验现在，研究椭圆方程问题u-u=f()x,y，运用五点差分格式对其求解，即：1h22ui-1,j-1h21ui,j-1+()+2(1h22+1h21)ui,j-1h22ui+1,j-1h21ui,j+1=f()xi,yju()0,0=0,1 i M2-1,1 j M1-1且该方程的精确解为sin()x sin()y。即探究椭圆方程u-()2ux2+2uy2=(+22)sin()x sin()y,0 x 1,0 y 1,u()0,0=0对上述椭圆方程五点格式使用上面介绍的Gauss-Seidel迭代法，可以得到如下的误差图（见图1、图2）（其中图2中的圈表明的是等高线）：000.511.520.511.5201234567810310376543211.81.61.41.210.80.60.40.20.20.40.60.811.21.41.61.8图1 Gauss-Seidel迭代法的误差分布图图2 Gauss-Seidel迭代法的误差等高线图-62焦岑，孙唯唯，聂家升：椭圆方程五点格式的迭代法与快速算法的比较4.1 数值算例4.1.1 数值算例1当=1，M2,M1为变量时（1）对于Gauss-Seidel迭代法，计算时间以及精确解与数值解的误差如表1所示：表1 Gauss-Seidel迭代法的相关数据结果1M22050100200300M12050100200300迭代步数K2929922 9279 38824 483计算时间0.018 3110.205 4472.067 19025.429 317174.608 371精确解与数值解之差0.007 10.001 23.1918e-046.4363e-053.3215e-05（2）对于快速Poisson算法，计算时间与误差如表2所示：表2 快速Poisson算法的相关数据结果1M220501002005001 0001 500M120501002005001 0001 500计算时间0.000 3980.002 5410.002 7190.006 8350.055 6610.258 0030.679 240精确解与数值解之差0.007 10.001 20.000 306 937.7503e-051.2475e-053.1250e-061.3898e-06从表1和2可以看出，当不变时，随着M2、M1的增加，Gauss-Seidel迭代法的计算时间和迭代步数都在增加，计算时间增加得尤为明显，当M2、M1增加到102数量级时，使用Gauss-Seidel迭代法计算时，就会发现其计算速度就明显地变得非常慢，比快速Poisson算法慢了103个数量级，到了300时，计算时间已达到了6分钟，此时，迭代法便不再适用了。而，快速Poisson算法的计算时间虽然也在增加，但是它的计算时间远远没有Gauss-Seidel迭代法增加得那么快，即使到了103数量级，计算时间仍然很小，变化不是很明显。同时，也可绘制出Gauss-Seidel迭代法和快速Poisson算法的计算时间随M2、M1的变化的图像。可以从图3了解到，当M2、M1很小时（小于102），Gauss-Seidel迭代法的计算时间很短，肉眼几乎看不出来差距，这时，利用 Gauss-Seidel 迭代法是相当有效的，而当M2、M1超过 100 时，Gauss-Seidel迭代法的迭代时间迅速增长，这时，Gauss-Seidel就不再适合求解此问题了。而从图4可以看到，当M2、M1增大到500时，快速Poisson算法的计算时间仍然小于0.1 s，虽然同M2=M1=图3 Gauss-Seidel迭代法的计算时间变化图4 快速Poisson算法的计算时间变化-63第 43 卷第 2 期保山学院学报2024 年 4 月20时相比，计算时间增加了很多，但是，比起Gauss-Seidel迭代法，快速Poisson算法所使用的计算时间就显得非常地短，因而，快速Poisson算法更适合求解椭圆方程。而且，可以看到，即便M2=M1=1 500时（这时，Gauss-Seidel迭代法已经运行得非常慢了，大约需要几十分钟），快速Poisson算法的计算时间仍然不超过1s。因而，可以得出：当M2，M1 100时，Gauss-Seidel迭代法和快速Poisson算法对于求解椭圆方程均有效，而当M2，M1 100时，快速Poisson算法更适合求解椭圆问题。4.1.2 数值算例2当M2=M1=40，为变量时，（1）对于Gauss-Seidel迭代法，计算时间以及精确解与数值解之间的误差如表3所示：表3 Gauss-Seidel迭代法的相关数据结果125501002005001 0005 00010 000迭代步数K51048241414891552919108计算时间0.092 2890.078 9710.051 5060.024 5540.024 4480.014 6710.010 4450.004 8430.003 8650.001 688精确解与数值解之差0.001 90.001 80.001 65.5258e-043.2034e-041.7412e-047.3570e-053.7525e-057.6629e-063.8351e-06（2）对于快速Poisson算法，计算时间和误差如表4所示：表4 快速Poisson算法的相关数据结果125501002005001000500010000计算时间0.0003410.0004140.0003270.0003220.0003220.0003330.0003270.0004020.0004200.000323精确解与数值解之差0.00190.00180.00165.5300e-043.2201e-041.7544e-047.4167e-053.7800e-057.6787e-063.8469e-06从表3和表4可以看出，随着的增加，Gauss-Seidel迭代法的迭代步数、计算时间和精确解与数值解的误差都在逐渐减少，而在快速Poisson算法中，计算时间几乎不随的变化而变化，稳定在0.0 003-0.0 005之间，但，精确解与数值解的误差随着的增大而减少。因而，在步长不变的情况下，快速Poisson算法要明显地优于Gauss-Seidel迭代法。5 总结本文研究的是椭圆方程问题，通过对矩形区域上的带有Dirichlet边界条件下的Poisson方程建立五点差分格式，然后，分别利用Gauss-Seidel迭代法和快速Poisson算法求解，并对二者从M2、M1，不同的系数以及不同的误差精度，来探究Gauss-Seidel迭代法与快速Poisson算法的优劣。通过数值算例，可以明显地看出，相比快速Poisson算法，Gauss-Seidel迭代法的计算时间-64焦岑，孙唯唯，聂家升：椭圆方程五点格式的迭代法与快速算法的比较复杂度更高，而快速Poisson算法不失为解决相应问题的一种良策。参考文献：1 李庆扬，王能超，易大义.数值分析（第5版）M.北京：清华大学出版社，20082 张传林.数值方法M.北京：中国科学文化出版社，2001:80-150.3 V.L.Makarov,S.V.Makarov,M.M.Moskal kov.A fast algorithm for solving the dirichlet problem on hexagonal templet for the poisson equation in a rectangleJ.Journal of Mathematical Sciences,1995,77(05).4 Manideepa Saha,Jahnabi Chakravarty.Convergence of Generalized SOR,Jacobi and GaussSeidel Methods for Linear SystemsJ.International Journal of Applied and Computational Mathematics,2020,6(1-2).5 郝艳花.Jacobi迭代法与Gauss-Seidel迭代法J.山西大同大学学报（自然科学版），2017，33（05）：3-5.6 李焕荣.浅谈线性方程组的迭代求解J.重庆工商大学学报（自然科学版），2012，29（07）：28-32.7 杜衡吉，徐昆良.Jacobi和Gauss-Seidel迭代法求解线性方程组的分析及应用J.曲靖师范学院学报，2011，30（03）：46-50.8 白红梅.Jacobi迭代法与Gauss-Seidel迭代法的收敛性比较分析J.呼伦贝尔学院学报，2009，17（06）：55-58.Comparison of Gauss-Seidel Iterative Method and Fast Poisson Algorithm in Five-Point Format for Elliptic EquationsJIAO Cen,SUN Weiwei,NIE JiaSheng(Applied Technology College,Suzhou University,Suzhou Zhejiang 215325,P.R.China)Abstract:This paper mainly discusses the solving of the five-point scheme of the elliptic equation byusing the Gauss-Seidel iteration method and the fast Poisson algorithm,and compares the speeds of thetwo equations for solving the linear equations.We all know that in the large linear system of equationswhere the coefficient matrix is sparse,the iterative method is a good solution to this type of algorithm,mainly because given an initial vector,we can obtain a result of an arbitrary iteration,and the operationis simple,but whether the approximate solution obtained by the iterative method converges to the exactsolution,and whether the iterative method can also solve the system of equations if the linear equationshave a fast algorithm to take advantage,we must further compare.This paper mainly judges the advantages and disadvantages of Gauss-Seidel iterative method and fast Poisson algorithm by comparing different,different step sizes and different error requirements.Keywords:Five-point scheme for elliptic equation;Gauss-Seidel iterative method;Fast Poisson algorithm-65