求解二阶波动方程的三次样条差分方法

第27卷第l期 2011年2月 大 学 数 学 COLLEGEMATHEMATICS V01．27，№．1 Feb．2011 求解二阶波动方程的三次样条差分方法 齐远节1， 刘利斌2 (1．池州市第一中学，安徽池州247000；2．池州学院数学计算机科学系。安徽池州247000) [摘 要]有限差分法在求解二阶波动方程初边值问题过程中通常受到精度和稳定性的限制．本文对二 阶波动方程的时间、空间项分别采用i次样条公式进行离散，推导出精度分别为o(r2+hz)，o(rz+ht)， O(r4+h2)和0(r4+h4)的四种三层隐式差分格式。以及与之相匹配的第一个时间步的同阶离散格式，并采 用Fourier方法分析了格式的稳定性．数值实验表明，本文给出的四种格式中的两种与精细时程积分方法、 经典的C—N格式和高精度紧致隐式差分格式相比，计算量减小，精度更高，对长时间计算问题仍然可以获 得精度很高的数值结果．本文格式包含了高精度紧致隐式差分格式的一种格式作为特例． [关键词]波动方程；三次样条；差分方法；稳定 [中图分类号]0241．8[文献标识码]A [文章编号]1672—1454(2011)01—0059—06 1 引 言 二阶波动方程是用于描写振动现象和波在介质中传播过程的数学物理方程．二阶波动方程出现在 机械振动、声波、弹性波的理论中．关于波动方程初值或初边值问题数值解法的研究一直是微分方程数 值解方向的重点课题．目前，对于这类方程的数值解法已经有一些结果，其中最常见的是有限元方法 和有限差分方法【lj．文I-1]中提出了一个三层的显式差分格式，其局部截断误差阶仅为O(r2+h2)．文 E23提出了一类显式辛格式，它比隐式辛格式的计算量小，但格式是条件稳定的，并且为二阶精度．文 [3]针对二阶波动方程提出了精细时程积分法，该方法虽然在时间方向是精确计算，但是，在空间方向 的局部阶段误差阶仅为O(h2)，且在实际计算中h的值不能太小，随着矩阵H阶数的增大，计算量增 大．文[4]的作者考虑了经典的C-N格式，并在此基础上设计了重叠型区域分解的并行算法．由于G N算法的局部截断误差阶仅为O(r2+h2)，所以文[4]的方法虽然适合并行计算，但是其精度太低．文 [5]中作者提出了两种隐式紧致差分格式，其局部误差阶分别为O(r2+h4)，0(r4+h4)． 本文利用三次样条的基本公式，给出了四种求解二阶波动方程的三点、三层隐式差分格式，其局 部截断误差分别为0(r2+h2)、0(≯+h4)、0(r4+h2)和O(r4+h4)．一般利用三次样条构造的差 分格式每个时间步需要耦合联立求解三个三对角方程组，而本文只需求解一个三对角方程组，工作量 大大降低，且本文包含了文[5]给出的精度为0(一+h4)的格式． 2样条基本公式及精度的研究‘6] 2．1 基本思想和均匀划分下的三次样条基本公式． 三次样条近似的基本思想是：分段用三次曲线来逼近真实解，在不同的段上它们一般是不同的， 但必须满足插值条件和连接条件．从而导出在节点处的函数值、一阶导数值及二阶导数值之间的基本 [收稿日期]2008—04—23 [基金项目]引进研究生项目(2009RC009) 万方数据 60 大 学 数 学 第27卷 关系式．把上述三类值作为未知量，直接代入微分方程中，结果得到一个有限的线性方程组，通过解 线性方程组而求出微分方程的近似解． 设置为区间[o，L]上的节点， 0一。o0， (11) df dZ “(z，o)一fl(Jc)，堕唔盟一^(z)，00． (13) 其中u(x，￡)是关于z与t的二元未知函数．取时间步长r，空间步长h=L／N，在节点(z，，t。) =(ih，nr)处的数值解记为“7．下面利用三次样条公式(1)和(10)式，将(11)式分别在时间和空间上离 散，可以得到如下四种个差分格式： 3．1．1二阶差分格式． 利用三次样条公式(1)，将(11)式在时间和空间上离散，得 (1一r2)“nHq-1+(4+2r2)Mrl+(1一r2)“尊} 一(2+4r2)U”+-l一(8r卫一8)“?+(2+4，)“耸l 一(1一r2)“=}一(4+2r2)“r1一(1一r2)“搿， (14) 其中，．一r肛，i=1，2，⋯，N一1，易知格式(14)的局部截断误差阶为O(r2+h2)． 3．1．2四阶差分格式(I)． 将(11)式在时间上采用(1)式离散，在空间上采用(10)离散，得 (1—2，上)“r，d-lI+(10+4r2)“P1+(1—2r2)“搿 =(2+8r2)“翟1一(16，一20)“?+(2+8r2)M备1一(1—2r2)“身 一(10+4r2)“r1一(1—2r2)扎斜， (15) 其中i=1，2，⋯，N一1，易知格式(15)的局部截断误差阶为0(r2+h4)． 3．1．3四阶差分格式(Ⅱ)． 将(11)式在时间上采用(10)式离散，在空间上采用(1)离散，得 (r2—2)“r卜e'+-1I一(8+2r2)“广1+(产一2)“刿 一一(4+lOr2)甜翟l+(20r2—16)“?一(4+lOr2)“～1 +(2一r2)“爿+(8+2r2)“r1+(2一r2)“T件r-}，(16) 其中i=l，2，⋯，N一1，易知格式(16)的局部截断误差阶为O(r4+h2)． 3．1．4四阶差分格式(111)． 将(11)式在时间和空间上都采用(10)离散，得 (1一r2)U“r-11+(10+2，)“r1+(1一r2)“搿 =(2+lOt2)“譬l一(20r2—20)M?+(2+10r2)“聋l 一(1一r2)U—r--11一(10+2r2)“●1一(1一r2)“料，(17) 其中i一1，2，⋯，N一1，易知格式(17)的局部截断误差阶为O(r4+h4)，该格式与文献[5]提出的精度 也为0(r4+h4)的紧致差分格式一致． 3．2初始条件的离散L5J． 因为以上四个格式是三层的，即在进行每一次时间推进时都需要知道前2步上的值．初始时刻的 值由(12)式的第一式给出，第一个时间步上的值将由(12)式的第二式给出．为此，需要对(12)式的第 二式也进行离散． 对于格式(14)一(17)，为了不影响其整体精度，并避免在第一时间步上引入高频误差，对(12)式 的第二式采用高阶离散．为此，利用泰勒公式将“；在U?处做四阶展开，可得 万方数据 62 大 学数 学 第27卷 一1卅r(au)0r2一(雾)；+吾(a"3u／,、0i+叫，． 利用(11)一(13)式可将(18)式改写为 一1 cmmc肌+萼(券)：+吾未(象)：+O(r4+列，． 将(19)式中的时间导数项用向前差分、空间导数项用中心差分离散后，略去高阶项，整理可得 一r2“Ll+2(3+r2)Uj—rz“～1 =2(3—2，)(，1)i+6r(厂2)f+2r2((，1)斗1+(，1)卜1)． 上式即为与上述四种格式相匹配的初始条件的离散格式． 4稳定性分析 (18) (19) (20) 为了证明上述四种格式的稳定性需要引用如下引理． 引理(Miller准则)[83实系数二次方程 Az2+Bx+C=0，A>0 的两个根按模小于等于1的充要条件是 A—C≥0，A+B+C≥0，A—B+C≥0． 定理格式(14)格式(17)稳定的必要条件是r≤1，而格式(15)稳定的必要条件是r≤√2，格式 (16)稳定的必要条件是r≤√2／2． 证这里我们只证明格式(15)的稳定性条件，其他格式类似． 由Fourier分析法可知格式(15)的特征方程为 m2+132+C一0， 其中 A=10+2cos口+4r2(1一cos疗)， B=一4cos口一20+16r2(1一COS口)， C=10+2cos口+4r2(1一COS口)． 显然，由引理可知A，B，C满足 A>0，A—C一0， A+B+C=24r2(1一COS口)≥0， A—B+C=40+8cos口一8r2(1一COS口)． 所以，要使 A—B+C=40+8cos口一8r2(1一COS口)≥0， 则，．≤厄．因此，格式(15)稳定的必要条件是r≤√虿． 5数值实验与结果分析 例[3]在初值问题(1‘1)一(13)中取 L一1，fl(z)一sin艘+1，f2(z)一0，gl(￡)=gz(￡)一1， 下面，用本文格式(14)、格式(15)和格式(17)进行计算，并与文[3]中的精细时程积分法、经典的GN 方法n1和文[5]中的格式(I)进行比较，计算结果列于表1—4．表中的精确值由解析解 U(z，￡)=1+sin7口：COSa't 算出．由于解在z∈[o，1]内关于z=0．5对称，所以，表中只列出0≤z≤0．5范围内的结果． 万方数据 第1期 齐远节，等：求解二阶波动方程的三次样条差分方法 63 表1取r=h=0．05计算100步时的数值结果 文[3]方法 工 精确解 本文格式(15) GN格式‘4] N≥30 N一20 O．1 O．69098 O．69102 O．69104 O．69121 O．69133 O．2 O．41221 O．41228 0．41232 O．41264 O．41288 O．3 O．19098 O．19108 O．19114 O．19153 O．19190 0．4 0．04894 O．04906 0．04913 0．04951 O．05002 0．5 0．00000 0．00012 0．00018 0．00062 0．00113 表2取r=h=0．05计算1000步时的数值结果 Z 精确解 本文格式(15) 文E3]方法(N=30)C-N格式‘4’ 0．1 1．30901 1．30501 1．30507 1．27384 0．2 1-58778 1．58017 1．58019 1．52088 O．3 1．80901 1．79854 1．79858 1．71693 O．4 1．95105 1．93874 1．93879 1．84280 O．5 2．00000 1．98705 1．98708 1．88618 表3取r=h=0．01计算1000步时的误差比较 Z 本文格式(14) 本文格式(15) 本文格式(17) C-N格式。4、 文Es]格式(I) O．1 O．12434F一130．25733P—06O．12434e—13O．231547P—05O．064302e一04 0．2 0．02664e一130．48948e—06O．02664P一13O．440430e—050．122311e一04 0．3 O．10658e一13O．67372e—06O．10658e一130．606200e—05O．168347e一04 O．4 0．24646e一130．79200P—06O．24646P—13O．712631P—05O．197903e一04 O．5 O．26645e一13O．83276e—06O．26645e一13O．749304e—050．208088e一04 表4取r一0．01。h=0．05计算1000步时的误差比较 Z 本文格式(14) 本文格式(15) 本文格式(17) C～格式143 文Es]格式(I) O．1 0．14829le一30．273446e—06O．244503e—09O．187407P—03O．065097P一04 0．2 0．282067e一30．520126P一060．465145P一09O．356470e—03O．123823e一04 O．3 0．388232F一3O．715893e—060．640193e—09O．490639e—03O．170428e一04 0．4 0．456395e一3O．841582e—060．752516e—09O．576781e～030．200350e一04 0．5 0．479882e一30．884892e～06O．791223e—090．606463e一03O．210661e—04 通过前面的理论分析与上面的表1—4可以得到如下结论： (i)本文提出了一维二阶波动方程的四种三次样条隐式差分格式，当网格比，．一1时，本文格式 (14)与格式(17)为同一个三层的显式差分格式．四种格式都是条件稳定的，其中格式(15)的稳定性条 件最弱，稳定的条件是r≤√2．本文的四种差分格式在每一个时间层上只有3个网格点，所以差分方 程是三对角型，因而可直接采用追赶法进行求解，无需迭代，计算简便． (ii)文[4]的重叠区域分解并行算法是在经典的C—N格式基础上设计的，因此其误差阶仍然仅为 O(r2+^2)，虽然适合并行计算，但是其精度太低，不适合实际问题的计算．上面的数据结果验证了此 结论． (iii)当网格比取r=l时，本文格式(14)和(17)为同一个具有四阶时间和空间精度的差分格式，与 万方数据 64 大 学数 学 第27卷 文[3]的精细时程积分方法相比，精度有了很大的提高．另外文[3]的方法在实际计算中h的值不能太 小，否则计算量将增大许多，而在实际计算中，我们希望h的值越小越好．显然，本文给出的四种三次 样条差分方法能满足此

要求

对教师党员的评价套管和固井爆破片与爆破装置仓库管理基本要求三甲医院都需要复审吗

． (iv)本文给出的四种格式包含了文[5]的高精度紧致差分格式(II)作为特例，由表3和表4结果 可看出，本文的格式(15)和格式(17)比文[5]的格式(I)的精度更高． [1] [2] [3] [4] [5] [6] [7] [8] [参考文献] 黄明游．发展方程数值计算方法[M]．北京：科学出版社，2004． 孙狄．波动方程的一类显式辛格式[J1．计算数学，1997，19(1)：1—10． 金承日，吕万金．二阶双曲型方程的精细积分法口]．计算力学学报，2003，20(1)：113—115． 田敏。羊丹平．波动方程的重叠型区域分解并行有限差分算法[J]．山东大学学报(理学版)，2007，42(2)： 28—38． 葛永斌，朱琳，田振夫．求解波动方程的高精度紧致隐式差分方法[J]．宁夏大学学报(自然科学版)，2005， 26(4)：297—299． 王璞．流体力学问题的三次样条配置法[J]．力学进展，1990，20(3)：316—327． RubinSG，KhoslaPK．Higher-ordernumericalsolutionsusingcubicsplines[J]．AIAAJ．，1976，14(7)：851--858． MillerJJ H，Onthelocationofzerosofcertainclassesofpolynomialswithapplicationtonumericalanalysis[J-]． J．Inst．Math．Appls，1971，8(3)：394—406． CubicSplinesDifferenceMethodforSolving theSecondOrderWaveEquation QfYuan-jiel，L，ULi—bin 2 (1．TheNumber0neMiddleSchoolofChizhou，Chizhou，Anhui24700，China； 2．DepartmentofMathematicsandComputerScience，ChizhouCollege，Chizhou，Anhui247000，China) Abstract：Insolvinginitialboundaryvalueproblemofthesecondorderwaveequation，theclassicalfinitedifference methodisgenerallyrestrictedbyitsstabilityandprecision．Inthispaper，fourclassesofthree-levelimplicitschemesare proposedforsolvingtheone-dimensionalwaveequationbyusingcubic—splinefunction．Thosemethodsareoforder O(r2+h2)，O(r2+h4)，O(r4+h2)and0(r4+h4)respectively．StabilityconditionsareobtainedbyFourieranalysis method．Itisshownbynumericalexamplesthatthetwoschemespresentedinthispaperaremuchbetterthantheprecise time—integrationmethod。theclassicalC—Nmethodandthehighaccuracycompactschemes，andtheyhavehighaccuracy evenforalongtimecalculation．Inaddition，thecurrentschemesincludethehighaccuracycompactschemeasaspecial case． Keywords：waveequation；cubicspline；differencemethod；stability 万方数据 求解二阶波动方程的三次样条差分方法 作者： 齐远节， 刘利斌， QI Yuan-jie， LIU Li-bin 作者单位： 齐远节,QI Yuan-jie(池州市第一中学,安徽,池州,247000)， 刘利斌,LIU Li-bin(池州学院 ,数学计算机科学系,安徽,池州,247000) 刊名： 大学数学 英文刊名： COLLEGE MATHEMATICS 年，卷(期)： 2011,27(1) 参考文献(8条) 1.Miller J J H On the location of zeros of certain classes of polynomials with application to numerical analysis[外文期刊] 1971(03) 2.Rubin S G;Khosla P K Higher-order numerical solutions using cubic splines 1976(07) 3.王璞 流体力学问题的三次样条配置法 1990(03) 4.葛永斌;朱琳;田振夫 求解波动方程的高精度紧致隐式差分方法[期刊论文]-宁夏大学学报（自然科学版） 2005(04) 5.田敏;羊丹平 波动方程的重叠型区域分解并行有限差分算法[期刊论文]-山东大学学报（理学版） 2007(02) 6.金承日;吕万金 二阶双曲型方程的精细积分法[期刊论文]-计算力学学报 2003(01) 7.孙狄 波动方程的一类显式辛格式 1997(01) 8.黄明游 发展方程数值计算方法 2004 本文链接：http://d.g.wanfangdata.com.cn/Periodical_dxsx201101015.aspx