重调和方程的二阶边值问题的求解

1、微分方程数值解课程设计第9页孝甫理N/tSouthChinaUniversityofTechnology微分方程数值解课程设计重调和方程的二阶边值问题的求解学院、系：理学院数学系专业：信息与计算科学姓名：陈剑宇学号：201030470140任课教师：黄凤辉提交日期：2012/12/27总评成绩：摘要本次课程设计，主要讨论重调和方程二阶边值问题的求解。文章分成四部分：第一部分，介绍如何将重调和方程的二阶边值问题分解，进而用经典的五点差分格式逼近；第二部分，列出对应的MATLAB算法和流程图；第三部分，通过上述方法求解具体的问题，并分析方法的精度和收敛阶；第四部分，将方法推广应用到抛物型方程，然后

2、分析稳定性。关键词：重调和方程的二阶边值问题；五点差分格式；MATLAB;抛物型方程目录引言1.1 背景41.2 问题的提出5二方法介绍2.1 五点差分方法62.2 二阶边值问题分解82.2.1 当a2-4bR0时82.2.2 当a2-4b<0时911162.3 推广到抛物型方程10三算法程序3.1 程序列表113.2 算法流程图与介绍123.2.1 重调和方程二阶边值程序123.2.2 抛物型方程程序14四结果与分析4.1 误差对比方法判断重调和方程算法的收敛阶164.1.1 a2-4bR0的误差对比方法164.1.2 a2-4b<0的误差对比方法184.1.3 结论204.2

3、图像方法判断重调和方程算法的收敛阶214.2.1 a2-4bR0的图像方法214.2.2 a2-4b<0的图像方法234.2.3 结论244.3 特例（大步长，高精度）254.3.1 特殊例子254.3.2 结论264.4 抛物型稳定性条件测试264.4.1 验证稳定性条件264.4.2 结论28五总结29六参考文献30引言1.1 背景微分方程是构造力学等领域的数学模型的主要方法。一般来讲，无论是物体运动轨迹或是流体速度测定模型都需要用到数值方法去求解。通过对椭圆型、抛物型和双曲型方程的研究和探讨，可以和现实中许多实际问题相互挂钩，并且得到解决问题的方法。许多的数学家都致力于其中，如欧拉

4、、柯西、贝努利、拉格朗日等人都为之做出了重要贡献。有限差分法、有限元法等等，为我们提供了求解的方法和手段。通过对微分方程数值解的学习和研究，我们可以得到许多实际问题的求解方法。虽然大部分问题求出来的是近似解，但是只要精度够高，能满足现实中的需求，这就足以体现它的重要性。在偏微分方程的求解问题中，Possion方程第一边值问题占据重要位置，五点差分格式等方法为我们提供了求解的方法。而在实际问题中，重调和方程的二阶边值问题也相当重要。于是，我们提出如下问题。1.2 问题的提出重调和方程的二阶边值问题2lu-atubu=fin11*u=g1onr(1)u=g2onT其中是Laplace算子，c是二维

5、平面上的有限区域，F是其光滑边界，a,b非负常数。数值算列：1. Q=0,nx0,冗,f=(4十2a+b)sinxsiny,使问题(1)存在精确解u=sinxsiny。2. Q=0,nx0,n,f=(1+a+b)(xsiny十ysinx),问题(1)存在精确解u=ysinx+xsiny分别取不同的a,b值计算，并且包括a2-4b之0和a2-4b<0两种情形。并将上面的方法推广应用到抛物型方程11121314u2.u-aubu=fin11:tIonu=gonu=g2u|J数值算列，Q=(xy,t):(xy)%。产产0产。0T,f=（5+2a+b）sinxsinyet,使问题（11）-（14

6、）存在精确解u=sinxsinyet二方法介绍由于区域的不同，会导致算法格式的变化，所以这里为了说明更简单，我们固定区域为矩形区域。其他的一些区域，可由矩形区域变化而来，或者进一步讨论即可得到差分格式，我们这里不再讨论。2.1 五点差分方法考虑矩形区域Q=0<x<a,0<y<b上，二阶线性椭圆型方程Lu三一三（p型）+工（p史）+qu=f,其中（x,y）CQ二x二xcycy第一边值问题。假设矩形区域。网格剖分均匀：h1=a/M,h2=b/N。于是网域包含（M-1）*（N-1）个内点，且均为正则内点。设（i,j）CQ,并且u（x,y）充分光滑，则沿着x和y方向分别用中心差

7、商代替导数有u u(p ) x x ijui 1,j -ui,jh1uij uij 2Pi.2,j)0/1)L,L,u(p )y y jX(ui,j 1ui,jh2 Pi,j 2 h2ui,jui,jh22)0(h2)这里uj表示u（xyj。而qu和右端项f直接有"qjulj，仁fj于是方程在方程的（i,j）点被表示为:2LUj=Lhuij(h12h2)当我们略去截断误差Ri,j(u)=o(hi2+h22),就可以得到五点差分格式:11Lhuij2、x(Pij、Juj2、y(Pijyuij)qjuijhlh2其中6x(Pij6xuj)=P+j(u*j-u，-Pi.2，j(ui,j-u

8、i,j)y(Pjyuj=Pi,j：(ui,ji-ui,j)-PhjJ.(ui,j-ui,jJ)22由方程我们可以看到，点(i,j)与相邻的四个点(i-1,j)、(i+1,j)、(i,j+1)、(ij-1)都有关系，于是称之为五点差分格式。i,j+1i-1ji,ji+1,ji,j-1图一五点差分格式(i,j)关系点图现在我们再进一步简化问题，令P(x,y)=1,q(x,y)=0时，方程变为Poisson-2-2uu/、方程Lu三一一2一一2=f(x,y)x二y而我们的五点差分格式也相应地变为1122一小5和+U»j)2(5书+Ui,j')(一+2)Ui,j=f.h12h22h1

9、2h22i,j又由于我们题目中的求解域刚好是个正方形，所以我们可以把步长h1=h2=h,用正方形网格剖分区域后，五点差分格式简化为1-;j(Ui+,j+U5+Ui,j书+U“-4Ui,j)=fijhi,j方程就是我们问题中所运用到的最简化五点差分格式。下面的二阶边值问题和推广到抛物型方程都需要利用到它。2.2二阶边值问题分解重调和方程的二阶边值问题2.2UuaAu+bu=finQ*u=g1onr（1）u=g2onT其中是Laplace算子，c是二维平面上的有限区域，F是其光滑边界，a,b非负常数。虽然直接利用差商逼近导数的方法可以得到相应的差分格式，但由于涉及的离散数数较多，从而边界点处的计算

10、较为困难。于是，我们采用另一种方法，适当地将其分解成两个Possion方程求解。我们使用两种完全不同的方法去逼近问题（1）,而具体运用哪一种方法与a2-4b有关。当a2-4b>0,我们立刻可以把问题（1）分解成两个二阶方程，然后运用五点差分格式进行求解。而当a2-4b<0时，我们提出了一种方法，把问题（1）简化为顺序迭代二阶方程，然后逐步求解。而这个方法的收敛性也已经被实验证明。2.2.1 当a2-4bA0时令 R =1(a + Ja2 -4b), 2那么边值问题（1）可转化为下列两个边值问题v -bv = f1v ="2 _g1in 1 on -in j on -(2)

11、(3)Dirichlet边值问题（2）及（3）可以用有限元法或者有限差分法求解，而程序中用的是五点差分格式逼近，从而求出问题（1）的数值解。2.2.2当a2-4b<0时令u=v-bu那么边值问题（1）可转化为下列两个边值问题(4)：v-av=finL'v=g2on-u=vin11u=gion(5)(6)其中函数u,vW未知且满足关联式（4）,所以（5）很难求解，但是若中已知和（6）就是易求解的Poisson问题。下面我们采用迭代法计算：先估计那么（5）,再求解Poisson问题(5)和(6)得v,u。1.给定中(0),如;:(0)-0in'J2.已知中，求解两个Poiss

12、on问题Rk)-avk)=f,：:(k)in1Jk)=g2on(8)(k)(k)J&u=vinC(k)u=g1on-3.计算新的迭代值3k4)(k1)(k)(k)呼=(1。)中-bouinQ(10)其中切为迭代参数。假设G=0,1im0,12,空间步长h,=l1/N,h2=l2/M,最优参数切=一11b(。C2)其中C1=23(3+“c2一22(p2+a)'1=gsin2-h1-4sinh122l1h222h22l2:2=-4-cos2h2二h142二cos21ih2h22l2方法具有二阶精度微分方程数值解课程设计2.3推广到抛物型方程现在将上面的方法推广应用到抛物型方程1)：

13、u2一u-aubu=fin11:t第ii页u=gU=g2ulJonon121314将方程(11)时间离散得n1nu-u.2n1,n1n1n1u-aubu=fA2n1.n11n1n:;11nu-a-u(b-)ufuTT.n1n1n1-aubuf(15)其中b1n-1=b+一,fTn-11nfuT得到（15）之后，我们可以发现它与问题（1）类似，于是我们可以运用上面的方法求解。从底层算起，然后计算下一个时间层，如此下去，知道计算完全部的时间层。三算法程序3.1 程序列表表一第一、二题相关的程序列表文件名标注mainl.m第一题头文彳%输入main1(a,b,h)可得到近似解main2.m第二题头文

14、彳%输入main2(a,b,h)可得到近似解main3.m附加例子头文件，输入main3(a,b,h)可得到近似解Fivepointsl.m当满足条件a2-4b>0的差分算法Fivepoints2.m当满足条件a2-4b<0的差分算法f.m计算右端项f的函数fl.m计算边值g1的函数f2.m计算边值g2的函数go.m1计算边值1g2g1的函数bes.m计算最优参数0Iteration.m迭代计算得到满足精度的中Correct.m计算精确解Testerr.m计算平均误差和最大误差表二第三题相关的程序列表文件名标注main4.m第三题头文彳输入main4(a,b,h,t)可得到近似解P

15、arabolic.m计算每一层的近似解，由底层开始，逐步计算更高的层级reFivepoints1.m专门用于第三题，其余类似Fivepoints1.mreFivepoints2.m专门用于第三题，其余类似Fivepoints2.mIteration1.m专门用于第三题，其余类似Iteration1.mCorrect1.m专门用于第三题，其余类似Correct1.mf.m、f1.m、f2.m、go.m、bes.m、Testerr.m用途如上表表二试验列表文件名标注testl.m得到a2-4b>0的差分算法的结果，并用第一种方法计算收敛阶test2.m得到a2-4b<0的差分算法的结果

16、，并用A种方法计算收敛阶test3.m检验如何得到较高精度的结果test4.m用第二种方法计算a2-4b>0的收敛阶test5.m用第二种方法计算a2-4b<0的收敛阶test6.m对第三题中，稳定性判断条件hA4>12/b的测试3.2 算法流程图与介绍3.2.1 重调和方程二阶边值程序图二重调和方程二阶边值程序图微分方程数值解课程设计计算公式参考2.2.2的方法介绍，上面就是算法的流程图。对应第一题、第二题和附加例子的程序。首先输入a,b,h到相应地题目，如main1（a,b,h）就是运算第一题。然后进行条件a2-4b的判断，当方0时，进入算法1（Fivepointsl）;

17、否贝进入算法2（Fivepoints2）。当a2-4b0,即算法1中，先对公式（2）进行五点差分算法，构造矩阵A1和右边项G1,得到（2）的解v。然后对公式（3）进行五点差分算法，构造矩阵A2和右边项G2（其中G2需要通过v的计算得到）。最后得到公式（3）的解u,这也是我们需要求的近似解。当a2-4b<0,即算法2中，类似地对公式（5）进行五点差分算法，构造矩阵A1和右边项G1。但是由于右边项G1与中有关，无法直接得到，所以我们运用迭代的方法。先令中（0）=0,然后得到G1,然后计算出公式（5）的解v。再根据公式（6）进行五点差分，得到（6）的解uo计算到这里，还要计算中（与u和於0）有

18、关），彳贸.（1）（1）若|中3（0）|满足精度要求，则认为u是近似解；否则把邛带入G1,重新计算，直至达到精度要求。下面给出邛的迭代算法：（其中Qu1就是中"*）,Qu2就是（k由）Qu1=zeros（M-1,N-1）;%令Qu为0,然后开始迭代%U,Qu2=Iteration（g1,Qu1,A1,a1,x1,x2,y1,y2,x,y,h,ma,b,M,N,l）;k=1;%迭代次数while（norm（Qu2-Qu1）>1.0e-6）&&（k<1000）%使精度达到要求，并且迭代次数小于1000Qu1=Qu2;U,Qu2=Iteration（g1,Qu1

19、,A1,a1,x1,x2,y1,y2,x,y,h,ma,b,M,N,l）;k=k+1;end%3.2.2抛物型方程程序图三抛物型方程程序图计算公式参考2.2.3的方法介绍，上面就是算法的流程图。由于算法和3.2.1大部分相同，所以这里主要介绍时间层的迭代问题输入main4(a,b,h,t)就会进入程序。算法开始会根据t来决定时间层一共有多少层，设定完毕后，就进入第一层的计算。之后就如同3.2.1计算，不过特别的是,右端项G1与上一层的近似解u有关系(当计算第一层时，G1可由初值得到)。得到第一层的解u1之后，令floor=floor+1,进入第二层的计算。如此迭代下去，得到全部层数的近似解。下

20、面给出时间层的迭代算法：%ifa*a-4*b<0forfloor=1:L%表示第几层T1=floor*t;%时间平面参数Gu2=reFivepoints2(x1,x2,y1,y2,T1,a,b,h,t,Gu1);Gu1=Gu2;Co=Correct1(x1,x2,y1,y2,h,T1);%第几层的精确解str2=sprintf('层数为：d',floor);disp(str2)Testerr(Gu1,Co,M,N);%计算floor层的误差endelseforfloor=1:L%表示第几层T1=floor*t;%时间平面参数Gu2=reFivepoints1(x1,x2,

21、y1,y2,T1,a,b,h,t,Gu1);Gu1=Gu2;Co=Correct1(x1,x2,y1,y2,h,T1);%第几层的精确解str2=sprintf('层数为：d',floor);disp(str2)Testerr(Gu1,Co,M,N);endend四结果与分析4.1 误差对比方法判断重调和方程算法的收敛阶由其他科学家的研究中，我们可以知道文章中两种逼近重调和方程的方法收敛阶均为0%2）。现在，我们来证明这个结论的正确性。因为方法收敛阶均为o（h2）,所以我们可以令a,b相同，h1=2*h2。假如输入hl的误差是输入h2的误差的四倍，那么就验证了我们想法。不过这种

22、判断方法有个相当不利的地方，就是需要步长足够小，这样才能使计算出来的数据较为精确。由于一般家用电脑把网格点取成60*60就已经要计算相当长的时间，所以这里只有把别有的结果展示出来。由于有两种方法，一种a2-4b>0,另外一种a2-4b<0,我们分开讨论。4.1.1 a2-4bA0的误差对比方法当运行程序test1.m之后，就可以得到如下表格。表四第一题，b=0、0.25、0.5的误差表格a2-4b>0第一>b=0ah=pi/10的误差E1h=pi/20的误差E2h=pi/40的误差E3E1/E2E2/E32.06.110002e-0031.381526e-0033.28

23、5226e-0044.4226474.2052692.35.966975e-0031.349345e-0033.208794e-0044.4221274.2051462.75.804773e-0031.312833e-0033.122066e-0044.4215644.2050123.05.700209e-0031.289285e-0033.066129e-0044.4212164.2049293.55.551360e-0031.255753e-0032.986462e-0044.4207434.2048174.05.427388e-0031.227813e-0032.920076e-0044

24、.4203714.204728a2-4b>0第一题b=0.25ah=pi/10的误差E1h=pi/20的误差E2h=pi/40的误差E3E1/E2E2/E32.05.922623e-0031.339536e-0033.185600e-0044.4213984.2049742.35.796431e-0031.311116e-0033.118084e-0044.4209904.2048772.75.652663e-0031.278724e-0033.041126e-0044.4205494.2047723.05.559609e-0031.257751e-0032.991293e-0044.4

25、202784.2047073.55.426637e-0031.227771e-0032.920051e-0044.4199114.2046204.05.315428e-0031.202688e-0032.860443e-0044.4196244.204552a2-4b>0第一题b=0.5ah=pi/10的误差E1h=pi/20的误差E2h=pi/40的误差E3E1/E2E2/E32.05.746394e-0031.300024e-0033.091838e-0044.4202234.2046962.35.635365e-0031.274994e-0033.032362e-0044.4199

26、164.2046232.75.508321e-0031.246343e-0032.964277e-0044.4195874.2045443.05.425778e-0031.227722e-0032.920023e-0044.4193854.2044963.55.307394e-0031.201008e-0032.856529e-0044.4191164.2044324.05.207995e-0031.178571e-0032.803196e-0044.4189074.204382分析：从表四最后两列可以看出，网格点较小的两次比较（E1/E2）结果接近4.40,而网格点较大的两次比较（E2/E3

27、）结果接近4.20。这说明网格点越小，结果会更加靠近4.00。表五第二题，b=0、0.25、0.5的误差表格a2-4b>0第二题b=0ah=pi/10的误差E1h=pi/20的误差E2h=pi/40的误差E3E1/E2E2/E32.01.029717e-0022.348624e-0035.597133e-0044.3843414.1961202.31.015493e-0022.316423e-0035.520536e-0044.3838864.1960112.79.993581e-0032.279877e-0035.433596e-0044.3833864.1958903.09.8895

28、34e-0032.256301e-0035.377504e-0044.3830734.1958153.59.741373e-0032.222717e-0035.297593e-0044.3826424.1957114.09.617922e-0032.194722e-0035.230974e-0044.3822974.195627a2-4b>0第二题b=0.25ah=pi/10的误差E1h=pi/20的误差E2h=pi/40的误差E3E1/E2E2/E32.09.987178e-0032.278699e-0035.430954e-0044.3828424.1957622.39.870067

29、e-0032.252149e-0035.367777e-0044.3825124.1956822.79.736602e-0032.221878e-0035.295739e-0044.3821504.1955953.09.650187e-0032.202271e-0035.249076e-0044.3819254.1955403.59.526657e-0032.174234e-0035.182343e-0044.3816164.1954654.09.423299e-0032.150766e-0035.126480e-0044.3813684.195405a2-4b>0第二题b=0.5ah=

30、pi/10的误差E1h=pi/20的误差E2h=pi/40的误差E3E1/E2E2/E32.09.695610e-0032.212894e-0035.274545e-0044.3814164.1954212.39.601012e-0032.191411e-0035.223405e-0044.3812014.1953692.79.492728e-0032.166811e-0035.164840e-0044.3809674.1953123.09.422346e-0032.150817e-0035.126759e-0044.3808224.1952773.59.321363e-0032.127861

31、e-0035.072099e-0044.3806254.1952284.09.236530e-0032.108570e-0035.026162e-0044.3804704.195190分析：从表四和表五我们可以看到误差虽然和a,b有关系，但是影响误差最大的因素还是步长。当步长减小为一半，误差可以变为1/4。这说明要有高精度的结果，要求步长足够小。由于表中E1/E2约为4.40,E2/E3约为4.20,和我们的目标结果4.00相比,显然不太理想。这是由于步长太小，使算法的其他影响因素过大，导致结果产生偏差。但是如同前文所说的，当网格点过多，家用电脑计算速度无法达到要求。所以在a2-4b<0

32、的方法中，将会引用其他论文中计算的结果。4.1.2a2-4b<0的误差对比方法下面给出其他论文中计算出来的结果。当然也可以运行文件中的test2.m程序,可是由于步长过大，不建议使用，所以不再列出，但读者可以自己测试一下。表六第一题，a=0、0.5、1的误差表格a2-4b<0第一题a=0b网格点为65*65网格点为129*129网格点为257*257E1/E2E2/E3迭代k次误差E1迭代k次误差E2迭代k次误差E30.353.7378e-459.3527e-552.3477e-53.99653.98380.763.4125e-468.4859e-562.0771e-54.0214

33、4.08551.073.2179e-478.0793e-572.0550e-53.98293.93151.582.9118e-487.2101e-581.7337e-54.03854.15882.092.6904e-496.8215e-591.8013e-53.94403.78702.5102.4561e-4106.0247e-5101.3909e-54.07674.3315a2-4b<0第一题a=0.5b网格点为65*65网格点为129*129网格点为257*257E1/E2E2/E3迭代k次误差E1迭代k次误差E2迭代k次误差E30.353.7378e-458.5291e-552.1

34、354e-53.99903.99410.763.4125e-467.9097e-561.9648e-54.00694.02571.073.2179e-477.5415e-571.8938e-53.99593.98221.582.9118e-486.9316e-581.7184e-54.00864.03382.092.6904e-496.4782e-591.6374e-53.98923.95642.592.4561e-496.1521e-591.6340e-53.93783.7651a2-4b<0第意a=1b网格点为65*65网格点为129*129网格点为257*257E1/E2E2/E3

35、迭代k次误差E1迭代k次误差E2迭代k次误差E30.353.1879e-457.9702e-551.9938e-53.99983.99750.762.9968e-467.4868e-561.8672e-54.00284.00961.062.8644e-467.1267e-561.7480e-54.01934.07711.572.6823e-476.7405e-571.7204e-53.97943.91802.082.5061e-486.2346e-581.5283e-54.01974.07942.592.3658e-495.9385e-591.5090e-53.98383.9354分析：从表六

36、可以看出计算结果与目标结果相当接近表七第二题，a=0、0.5、1的误差表格a2-4b<0第二题a=0b网格点为65*65网格点为129*129网格点为257*257E1/E2E2/E3迭代k次误差E1迭代k次误差E2迭代k次误差E30.365.4815e-461.3703e-463.4247e-54.00024.00120.775.0136e-471.2555e-473.1525e-53.99333.98261.084.7101e-481.1763e-482.9257e-54.00424.02061.594.2909e-491.0770e-492.7341e-53.98413.93912

37、.0103.9191e-4109.7264e-5102.3582e-54.02934.12452.5113.6420e-4119.2141e-5112.4078e-53.95263.8268a2-4b<0第二题a=0.5b网格点为65*65网格点为129*129网格点为257*257E1/E2E2/E3迭代k次误差E1迭代k次误差E2迭代k次误差E30.355.1875e-451.2985e-453.2589e-53.99503.98450.774.8209e-471.2061e-473.0184e-53.99713.99581.074.5846e-471.1495e-472.9067e

38、-53.98833.95471.584.2213e-481.0498e-482.5687e-54.02114.08692.093.9362e-499.9099e-592.5458e-53.97203.89262.5103.6555e-4109.0709e-5102.1945e-54.02994.1335a2-4b<0第二题a=1b网格点为65*65网格点为129*129网格点为257*257E1/E2E2/E3迭代k次误差E1迭代k次误差E2迭代k次误差E30.354.9843e-451.2472e-453.1232e-53.99643.99330.764.6839e-461.1696e

39、-462.9067e-54.00474.02381.074.4871e-471.1229e-472.8174e-53.99603.98561.584.1852e-481.0449e-482.5970e-54.00544.02352.093.9278e-499.8369e-592.4760e-53.99293.97292.593.7068e-499.3636e-592.4338e-53.95873.8473分析：从表六、表七可以看到E1/E2和E2/E3确实在4.00附近变动，这证明了我们猜想方法收敛阶为0,2）是正确的。而且还可以看到，1.当a,b固定的时候，迭代次数k与我们的网格点数无关，说

40、明步长对迭代次数k无影响；2.当a和网格点数固定,即a、h固定时，b越大，迭代次数k也越大。4.1.3结论从4.1的讨论当中我们可以知道一下几个结论:1 .步长的大小对误差影响十分大。步长越小，方法精度越高；2 .两种方法（a2-4b>0和a2-4b<0）收敛阶为02）；3 .误差对比方法确实可以得知收敛阶，但是前提是步长足够小，这就需要有较好的设备才能得出结论；4 .对于a2-4b<0的算法，我们发现I.当a,b固定的时候，步长对迭代次数k无影响；H.当a、h固定时，b越大，迭代次数k也越大。所以在满足a2-4b<0时，我们可以适当减小b使迭代次数减小，更快得到近似解

41、。4.2图像方法判断重调和方程算法的收敛阶现在介绍另外一种方法判断算法的收敛阶，这也是一种相当常用的方法来判断收敛阶。假设误差e=0（h2）,我们可令e=0（h2ch2,其中c为常数然后对上式两边自然对数，有lne=lnc+2lnh这条公式说明了算法的误差对数与步长对数呈线性相关。如果我们计算出来的误差对数线与y=2lnh相互平行，则证明了算法收敛阶为e=0（h2）。同样的，需要对a2-4b>0和a2-4b<0两种算法分析。4.2.1 a2-4b>0的图像方法运行程序test4.m即可得到下面两个图像：第33页13间区Eil* EdiL 近e Insert L&ols

42、 Oesklep JLiitdo* Halp卢P昌4 qq白®后咀_s_g第一题，足脸口收敛阶判断图+一误差对数曲畿产2*曲畿76.中2.52-1.5步长时数图四 第一题，a2- 4b>0收敛阶判断图分析：可以看到两条曲线是互相平行的。Figure 2E®区）Eilt Edit 於" Lrebt% 工 *q1e Resktop 心 nd。理 H*lp口方口鸟fe虱氢争® /口亘口第二题，辞比此却收敛阶判断图-6.54-1.5-1误差对数曲线 - o ¥乂曲线2 5 3 5 4 5 5 5 -2.-3.h4 £步长对数图五第二题，

43、a2-4b>0收敛阶判断图分析：图三和图四中的两条曲线都是互相平行的，说明了a2-4b>0收敛阶确实是叫"。4.2.2 a2-4b<0的图像方法运行程序test5.m,有如下图像：nltx。FigureIElleEditJfiewIjtEWFtZoelsHesktopWindowHelpBP昌玲或1磬®要口目L.n第一题，才去必<口收敛阶判断图_ _ _ _ _ _ _ _ _ 325354 6 556 512 3 5 6十一误差对数曲统k2)(曲线-1.0-17-1615-1J-1.3-12-11步长对数图六第一题，a2-4b<0收敛阶判断图

44、分析：在a2-4b<0的方法中，同样看到两条曲线平行。9-1.0-1.7-1 6-1 5 九4 d 3-1.2-1 1步长对数?Fi£ure2匚|回|瓦EileEditewLnsartToolsHeskioplindnHelp不屈昌玲/改已冥则qMu口第二题，a、处如收敛阶判断图误差对数曲线y=4曲线一。*图七第二题，a2-4b<0收敛阶判断图分析：由图六和图七可以看出，a2-4b<0方法的收敛阶也是0%2）。4.2.3结论由上面的分析，可以得知：1 .a2-4b>0和a2-4b<0两种算法收敛阶是0+2）；2 .如果觉得用看图来判断两条直线是否平行不太

45、严密的话，可以计算两条直线斜率是否相等。如果相等，则证明结论无误；否则说明收敛阶不是0（h2）。不过这样的话，又会和4.1的方法一样，出现因为步长不够小而导致结论不正确的情况，所以这里就不给出来了。4.3 特例(大步长，高精度)4.3.1 特殊例子在某一些特殊例子中，尽管步长较大，我们仍然可以取得高精度的结果。例如精确解为u=(x241y2-4；,在-2,2*-2,2的区域里面，而右端项f=8-2a(x2+y2-8)+b(x2-4)(y2-4)。这个时候我们能取到相当高精度的解，如表八和表九(运行test3.m可得到)。表八a2-4b>0,特殊例子的误差表格网格数为10*10ab=0.0

46、时的误差b=0.25时的误差b=0.5时的误差2.02.116277e-0157.779785e-0151.989629e-0142.32.812565e-0157.352144e-0151.744010e-0142.72.050486e-0152.675500e-0152.532953e-0153.02.412336e-0155.104285e-0153.234724e-015表九a2-4b<0,特殊例子的误差表格网格数为10*10b迭代次数ka=0.0时的误差迭代次数ka=0.5时的误差迭代次数ka=1.0时的误差0.3122.253198e-008101.825354e-00891

47、.282883e-0080.7243.223164e-008172.416017e-008142.030519e-0081.0423.718790e-008251.455840e-008191.322050e-0081.53013.218638e-008482.140648e-008301.504249e-008分析：上面两个表都是在网格点数为10*10的时候得到的，但是精度都分别达到了10*和10、与前面两题得出的结果相比准确很多(因为前面的两题尽管网格数为465*65,精度也只达到10)。这是因为对于任何步长大小，二次函数中心差分格式的逼近误差为零，所以使近似解有相当高的精度。对于这个例子

48、，下面给出一幅在a=0.5,b=1,网格点为20*20的函数图，图八可以见到x和y越趋向0,函数值就越大。图八a=0.5,b=1,h=0.2的函数图4.3.2 结论由上面可以知道：在某一些特例之中，我们的两种算法可以在步长较小的情况下,得到高精度的解。这些特例指的是二次函数。4.4 抛物型稳定性条件测试当推广到抛物型方程的时候，我们需要对其稳定性进行分析。在对算法拆分和用傅里叶方法计算稳定性之后，可以得到稳定性条件：412h-T卜面我们来验证其是否正确4.4.1 验证稳定性条件运行程序test6.m,可以得到如下的结果:表十和表4是不满足稳定性条件（即h4<12）的测试例子。hb层数平均

49、误差大小最大误差大小19.556834e-0022.136926e-00121.329571e-0012.972946e-00131.651503e-0013.692793e-00142.021831e-0014.520853e-00152.470265e-0015.523560e-001不稳定，且满足a*a-4b<0的表a=0.5b=1h=pi/20t=0.2有Q黄足a*a-4b<0)层数平均误差大小最大误差大小12.542859e-0015.685882e-00124.263645e-0019.533591e-00135.734781e-0011.282308e+00047.244495e-0011.619883e+00058.957729e-0012.002965e+000表十不稳定，且满足a*a-4b>0的表有Q黄足a*a-4b>0)a=10b=0.5h=pi/20t=0.2分析：可见数据误差大，而且每下一层，误差增速更快，两层之间的误差增长较