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

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

2、五点差分格式；MATLAB；抛物型方程目 录一 引言41.1 背景41.2 问题的提出5二 方法介绍62.1 五点差分方法62.2 二阶边值问题分解82.2.1 当a²4b0时82.2.2 当a²4b<0时92.3 推广到抛物型方程10三 算法程序113.1 程序列表113.2 算法流程图与介绍123.2.1 重调和方程二阶边值程序123.2.2 抛物型方程程序14四 结果与分析164.1 误差对比方法判断重调和方程算法的收敛阶164.1.1 a²4b0的误差对比方法164.1.2 a²4b<0的误差对比方法184.1.3 结论204.2 图

3、像方法判断重调和方程算法的收敛阶214.2.1 a²4b0的图像方法214.2.2 a²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 问题的提出重调和方程的二阶边值问题 （1）其中是Laplace算子，是二维平面上的有限区域，是其光滑边界，a，

5、b非负常数。数值算列：1. ，使问题（1）存在精确解。2. ，问题（1）存在精确解u=ysinx+xsiny分别取不同的a，b值计算，并且包括和两种情形。并将上面的方法推广应用到抛物型方程 数值算列， ，使问题（11）-（14）存在精确解。二 方法介绍由于区域的不同，会导致算法格式的变化，所以这里为了说明更简单，我们固定区域为矩形区域。其他的一些区域，可由矩形区域变化而来，或者进一步讨论即可得到差分格式，我们这里不再讨论。2.1 五点差分方法考虑矩形区域=0<x<a，0<y<b上，二阶线性椭圆型方程Lu=f,其中(x,y) 第一边值问题。假设矩形区域网格剖分均匀：h1=

6、a/M，h2=b/N。于是网域包含(M-1)*(N-1)个内点，且均为正则内点。 设（i，j），并且u(x,y)充分光滑，则沿着x和y方向分别用中心差商代替导数有 这里表示u(,)。而qu和右端项f直接有qu ，f= 于是方程在方程的(i,j)点被表示为：当我们略去截断误差=，就可以得到五点差分格式： 其中 由方程我们可以看到，点(i，j)与相邻的四个点（i-1,j）、（i+1,j）、（i,j+1）、（i,j-1）都有关系，于是称之为五点差分格式。i-1,ji,j+1i,j-1i+1,ji,j图一 五点差分格式(i,j)关系点图 现在我们再进一步简化问题，令p(x,y)=1,q(x,y)=0时

7、，方程变为Poisson方程 Lu 而我们的五点差分格式也相应地变为 又由于我们题目中的求解域刚好是个正方形，所以我们可以把步长h1=h2=h，用正方形网格剖分区域后，五点差分格式简化为 方程就是我们问题中所运用到的最简化五点差分格式。下面的二阶边值问题和推广到抛物型方程都需要利用到它。2.2 二阶边值问题分解 重调和方程的二阶边值问题 （1）其中是Laplace算子，是二维平面上的有限区域，是其光滑边界，a，b非负常数。虽然直接利用差商逼近导数的方法可以得到相应的差分格式，但由于涉及的离散点数较多，从而边界点处的计算较为困难。于是，我们采用另一种方法，适当地将其分解成两个Possion方程求

8、解。 我们使用两种完全不同的方法去逼近问题（1），而具体运用哪一种方法与a²4b有关。当a²4b0，我们立刻可以把问题（1）分解成两个二阶方程，然后运用五点差分格式进行求解。而当a²4b<0时，我们提出了一种方法，把问题（1）简化为顺序迭代二阶方程，然后逐步求解。而这个方法的收敛性也已经被实验证明。2.2.1 当a²4b0时令，那么边值问题（1）可转化为下列两个边值问题 （2） （3）Dirichlet边值问题（2）及（3）可以用有限元法或者有限差分法求解，而程序中用的是五点差分格式逼近，从而求出问题（1）的数值解。2.2.2 当a²4b

9、<0时令 （4）那么边值问题（1）可转化为下列两个边值问题 （5） （6）其中函数未知且满足关联式（4），所以（5）很难求解，但是若已知，那么（5）和（6）就是易求解的Poisson问题。下面我们采用迭代法计算：先估计，再求解Poisson问题（5）和（6）得,。1.给定，如 （7） 2.已知，求解两个Poisson问题 （8） （9）3.计算新的迭代值 （10）其中为迭代参数。假设, 空间步长,最优参数其中， 方法具有二阶精度。2.3 推广到抛物型方程现在将上面的方法推广应用到抛物型方程 将方程（11）时间离散得 （15）其中，。 得到（15）之后，我们可以发现它与问题（1）类似，于是

10、我们可以运用上面的方法求解。从底层算起，然后计算下一个时间层，如此下去，知道计算完全部的时间层。三 算法程序3.1 程序列表表一 第一、二题相关的程序列表文件名标注main1.m第一题头文件，输入main1(a,b,h)可得到近似解main2.m第二题头文件，输入main2(a,b,h)可得到近似解main3.m附加例子头文件，输入main3(a,b,h)可得到近似解Fivepoints1.m当满足条件a²4b0的差分算法Fivepoints2.m当满足条件a²4b<0的差分算法f.m计算右端项f的函数f1.m计算边值g1的函数f2.m计算边值g2的函数go.m计算边

11、值的函数bes.m计算最优参数Iteration.m迭代计算得到满足精度的Correct.m计算精确解Testerr.m计算平均误差和最大误差表二 第三题相关的程序列表文件名标注main4.m第三题头文件，输入main4(a,b,h,t)可得到近似解Parabolic.m计算每一层的近似解，由底层开始，逐步计算更高的层级reFivepoints1.m专门用于第三题，其余类似Fivepoints1.mreFivepoints2.m专门用于第三题，其余类似Fivepoints2.mIteration1.m专门用于第三题，其余类似Iteration1.mCorrect1.m专门用于第三题，其余类似C

12、orrect1.mf.m、f1.m、f2.m、go.m、bes.m、Testerr.m用途如上表表三 试验列表文件名标注test1.m得到a²4b0的差分算法的结果，并用第一种方法计算收敛阶test2.m得到a²4b<0的差分算法的结果，并用第一种方法计算收敛阶test3.m检验如何得到较高精度的结果test4.m用第二种方法计算a²4b0的收敛阶test5.m用第二种方法计算a²4b<0的收敛阶test6.m对第三题中，稳定性判断条件h412/b的测试3.2 算法流程图与介绍3.2.1 重调和方程二阶边值程序图二 重调和方程二阶边值程序图输

13、入a,b,h得到近似解u= A2G2,并与精确解比较，计算误差大小否是<00计算a²4b 计算公式（2），得到A1,G1v=A1G1得到G1，计算v=A1G1计算公式（3），得到A2,G2令Qu1=Qu2,初始Qu1=0u=A2/G2|Qu1-Qu2|达到精度要求？由公式（6），得到A2,G2计算公式（5），得到A1计算公式参考2.2.2的方法介绍，上面就是算法的流程图。 对应第一题、第二题和附加例子的程序。首先输入a,b,h到相应地题目，如main1(a,b,h)就是运算第一题。然后进行条件a²4b的判断，当0时，进入算法1（Fivepoints1）；否则，进入算法

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

15、|满足精度要求，则认为u是近似解；否则把带入G1，重新计算，直至达到精度要求。 下面给出的迭代算法：（其中Qu1就是，Qu2就是）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) %使精度达到要求，并且迭代次数小于1000 Qu1=Qu2; U,Qu2=Iteration(g1,Qu1,A1,a1,x1,x2,y1,y2,x,y,h,m

16、a,b,M,N,l); k=k+1;end%3.2.2 抛物型方程程序图三 抛物型方程程序图是否，计算下一层floor=floor+1输入a,b,h,t得到近似解u= A2G2,并与精确解比较，计算误差大小否是<00计算a²4b 计算公式（2），得到A1,G1v=A1G1得到G1，计算v=A1G1计算公式（3），得到A2,G2令Qu1=Qu2,初始Qu1=0u=A2/G2|Qu1-Qu2|达到精度要求？由公式（6），得到A2,G2计算公式（5），得到A1首先计算第一层的解，floor=1T= floor*t已经是顶层了吗？结束 计算公式参考2.2.3的方法介绍，上面就是算法的流

17、程图。 由于算法和3.2.1大部分相同，所以这里主要介绍时间层的迭代问题。输入main4(a,b,h,t)就会进入程序。算法开始会根据t来决定时间层一共有多少层，设定完毕后，就进入第一层的计算。之后就如同3.2.1计算，不过特别的是，右端项G1与上一层的近似解u有关系（当计算第一层时，G1可由初值得到）。 得到第一层的解u1之后，令floor=floor+1，进入第二层的计算。如此迭代下去，得到全部层数的近似解。下面给出时间层的迭代算法：%if a*a-4*b<0 for floor=1:L %表示第几层 T1=floor*t; %时间平面参数 Gu2=reFivepoints2(x1,

18、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层的误差 endelse for floor=1:L %表示第几层 T1=floor*t; %时间平面参数 Gu2=reFivepoints1(x1,x2,y1,y2,T1,a,b,h,t,Gu1); Gu1=Gu2; Co=Correct1(x1,x2,y1,y2,h,T1); %第几层的精

19、确解 str2=sprintf('层数为:%d',floor); disp(str2) Testerr(Gu1,Co,M,N); endend四 结果与分析4.1 误差对比方法判断重调和方程算法的收敛阶 由其他科学家的研究中，我们可以知道文章中两种逼近重调和方程的方法收敛阶均为。现在，我们来证明这个结论的正确性。因为方法收敛阶均为，所以我们可以令a,b相同，h1=2*h2。假如输入h1的误差是输入h2的误差的四倍，那么就验证了我们想法。 不过这种判断方法有个相当不利的地方，就是需要步长足够小，这样才能使计算出来的数据较为精确。由于一般家用电脑把网格点取成60*60就已经要计算相

20、当长的时间，所以这里只有把别有的结果展示出来。由于有两种方法，一种a²4b0，另外一种a²4b<0，我们分开讨论。4.1.1 a²4b0的误差对比方法 当运行程序test1.m之后，就可以得到如下表格。表四 第一题，b=0、0.25、0.5的误差表格a²4b0 第一题 b=0ah=pi/10的误差E1h=pi/20的误差E2h=pi/40的误差E3E1/E2E2/E32.06.110002e-0031.381526e-0033.285226e-0044.4226474.2052692.35.966975e-0031.349345e-0033.208

21、794e-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.4203714.204728a²4b0 第一题 b=0.25ah=pi/10的误差E1h=pi/20的误差

22、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.4202784.2047073.55.426637e-0031.227771e-0032.920051e-0044.

23、4199114.2046204.05.315428e-0031.202688e-0032.860443e-0044.4196244.204552a²4b0 第一题 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.4199164.2046232.75.508321e-0031.246343e-0032.964277e-0044.4

24、195874.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）结果接近4.20。这说明网格点越小，结果会更加靠近4.00。表五 第二题，b=0、0.25、0.5的误差表格

25、a²4b0 第二题 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.889534e-0032.256301e-0035.377504e-0044.3830734.1958153.5

26、9.741373e-0032.222717e-0035.297593e-0044.3826424.1957114.09.617922e-0032.194722e-0035.230974e-0044.3822974.195627a²4b0 第二题 b=0.25ah=pi/10的误差E1h=pi/20的误差E2h=pi/40的误差E3E1/E2E2/E32.09.987178e-0032.278699e-0035.430954e-0044.3828424.1957622.39.870067e-0032.252149e-0035.367777e-0044.3825124.1956822.7

27、9.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.195405a²4b0 第二题 b=0.5ah=pi/10的误差E1h=pi/20的误差E2h=pi/40的误差E3E1/E2E2/E32.09

28、.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.127861e-0035.072099e-0044.380625 4.1952284.09.236530e-

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

30、误差对比方法 下面给出其他论文中计算出来的结果。当然也可以运行文件中的test2.m程序，可是由于步长过大，不建议使用，所以不再列出，但读者可以自己测试一下。表六 第一题，a=0、0.5、1的误差表格a²4b<0 第一题 a=0b网格点为65*65网格点为129*129网格点为257*257E1/E2E2/E3迭代k次误差E1迭代k次误差E2迭代k次误差E30.353.7378e459.3527e552.3477e53.99653.98380.763.4125e468.4859e562.0771e54.02144.08551.073.2179e478.0793e572.0550

31、e53.98293.93151.582.9118e487.2101e581.7337e54.03854.15882.092.6904e496.8215e591.8013e53.94403.78702.5102.4561e4106.0247e5101.3909e54.07674.3315a²4b<0 第一题 a=0.5b网格点为65*65网格点为129*129网格点为257*257E1/E2E2/E3迭代k次误差E1迭代k次误差E2迭代k次误差E30.353.7378e458.5291e552.1354e53.99903.99410.763.4125e467.9097e561.9

32、648e54.00694.02571.073.2179e477.5415e571.8938e53.99593.98221.582.9118e486.9316e581.7184e54.00864.03382.092.6904e496.4782e591.6374e53.98923.95642.592.4561e496.1521e591.6340e53.93783.7651a²4b<0 第一题 a=1b网格点为65*65网格点为129*129网格点为257*257E1/E2E2/E3迭代k次误差E1迭代k次误差E2迭代k次误差E30.353.1879e457.9702e551.993

33、8e53.99983.99750.762.9968e467.4868e561.8672e54.00284.00961.062.8644e467.1267e561.7480e54.01934.07711.572.6823e476.7405e571.7204e53.97943.91802.082.5061e486.2346e581.5283e54.01974.07942.592.3658e495.9385e591.5090e53.98383.9354分析：从表六可以看出计算结果与目标结果相当接近。表七 第二题，a=0、0.5、1的误差表格a²4b<0 第二题 a=0b网格点为65*

34、65网格点为129*129网格点为257*257E1/E2E2/E3迭代k次误差E1迭代k次误差E2迭代k次误差E30.365.4815e461.3703e463.4247e54.00024.00120.775.0136e471.2555e473.1525e53.99333.98261.084.7101e481.1763e482.9257e54.00424.02061.594.2909e491.0770e492.7341e53.98413.93912.0103.9191e4109.7264e5102.3582e54.02934.12452.5113.6420e4119.2141e5112.40

35、78e53.95263.8268a²4b<0 第二题 a=0.5b网格点为65*65网格点为129*129网格点为257*257E1/E2E2/E3迭代k次误差E1迭代k次误差E2迭代k次误差E30.355.1875e451.2985e453.2589e53.99503.98450.774.8209e471.2061e473.0184e53.99713.99581.074.5846e471.1495e472.9067e53.98833.95471.584.2213e481.0498e482.5687e54.02114.08692.093.9362e499.9099e592.54

36、58e53.97203.89262.5103.6555e4109.0709e5102.1945e54.02994.1335a²4b<0 第二题 a=1b网格点为65*65网格点为129*129网格点为257*257E1/E2E2/E3迭代k次误差E1迭代k次误差E2迭代k次误差E30.354.9843e451.2472e453.1232e53.99643.99330.764.6839e461.1696e462.9067e54.00474.02381.074.4871e471.1229e472.8174e53.99603.98561.584.1852e481.0449e482.5

37、970e54.00544.02352.093.9278e499.8369e592.4760e53.99293.97292.593.7068e499.3636e592.4338e53.95873.8473分析：从表六、表七可以看到E1/E2和E2/E3确实在4.00附近变动，这证明了我们猜想方法收敛阶为是正确的。而且还可以看到，1.当a，b固定的时候，迭代次数k与我们的网格点数无关，说明步长对迭代次数k无影响；2.当a和网格点数固定，即a、h固定时，b越大，迭代次数k也越大。4.1.3 结论 从4.1的讨论当中我们可以知道一下几个结论：1.步长的大小对误差影响十分大。步长越小，方法精度越高；2.

38、两种方法（a²4b0和a²4b<0）收敛阶为；3.误差对比方法确实可以得知收敛阶，但是前提是步长足够小，这就需要有较好的设备才能得出结论；4.对于a²4b<0的算法，我们发现.当a，b固定的时候，步长对迭代次数k无影响；.当a、h固定时，b越大，迭代次数k也越大。所以在满足a²4b<0时，我们可以适当减小b使迭代次数减小，更快得到近似解。4.2 图像方法判断重调和方程算法的收敛阶 现在介绍另外一种方法判断算法的收敛阶，这也是一种相当常用的方法来判断收敛阶。假设误差e=， 我们可令e=ch²，其中c为常数 然后对上式两边自然对数

39、，有lne = lnc + 2lnh这条公式说明了算法的误差对数与步长对数呈线性相关。如果我们计算出来的误差对数线与y=2lnh相互平行，则证明了算法收敛阶为e=。同样的，需要对a²4b0和a²4b<0两种算法分析。4.2.1 a²4b0的图像方法 运行程序test4.m即可得到下面两个图像：图四 第一题，a²4b0收敛阶判断图分析：可以看到两条曲线是互相平行的。图五 第二题，a²4b0收敛阶判断图分析：图三和图四中的两条曲线都是互相平行的，说明了a²4b0收敛阶确实是。4.2.2 a²4b<0的图像方法运行程序

40、test5.m，有如下图像：图六 第一题，a²4b<0收敛阶判断图分析：在a²4b<0的方法中，同样看到两条曲线平行。图七 第二题，a²4b<0收敛阶判断图分析：由图六和图七可以看出，a²4b<0方法的收敛阶也是。4.2.3 结论 由上面的分析，可以得知：1.a²4b0和a²4b<0两种算法收敛阶是；2.如果觉得用看图来判断两条直线是否平行不太严密的话，可以计算两条直线斜率是否相等。如果相等，则证明结论无误；否则说明收敛阶不是。不过这样的话，又会和4.1的方法一样，出现因为步长不够小而导致结论不正确的情况

41、，所以这里就不给出来了。4.3 特例（大步长，高精度）4.3.1 特殊例子在某一些特殊例子中，尽管步长较大，我们仍然可以取得高精度的结果。例如精确解为，在-2,2*-2,2的区域里面，而右端项f=8-2a(x²+y²-8)+b(x²-4)(y²-4)。这个时候我们能取到相当高精度的解，如表八和表九（运行test3.m可得到）。表八 a²4b0，特殊例子的误差表格网格数为10*10ab=0.0时的误差b=0.25时的误差b=0.5时的误差2.02.116277e-0157.779785e-0151.989629e-0142.32.812565e-

42、0157.352144e-0151.744010e-0142.72.050486e-0152.675500e-0152.532953e-0153.02.412336e-0155.104285e-0153.234724e-015表九 a²4b<0，特殊例子的误差表格网格数为10*10b迭代次数ka=0.0时的误差迭代次数ka=0.5时的误差迭代次数ka=1.0时的误差0.3122.253198e-008101.825354e-00891.282883e-0080.7243.223164e-008172.416017e-008142.030519e-0081.0423.718790

43、e-008251.455840e-008191.322050e-0081.53013.218638e-008482.140648e-008301.504249e-008分析：上面两个表都是在网格点数为10*10的时候得到的，但是精度都分别达到了和，与前面两题得出的结果相比准确很多（因为前面的两题尽管网格数为65*65，精度也只达到）。这是因为对于任何步长大小，二次函数中心差分格式的逼近误差为零，所以使近似解有相当高的精度。 对于这个例子，下面给出一幅在a=0.5，b=1,网格点为20*20的函数图，图八。可以见到x和y越趋向0，函数值就越大。图八 a=0.5，b=1,h=0.2的函数图4.3.

44、2 结论 由上面可以知道：在某一些特例之中，我们的两种算法可以在步长较小的情况下，得到高精度的解。这些特例指的是二次函数。4.4 抛物型稳定性条件测试 当推广到抛物型方程的时候，我们需要对其稳定性进行分析。在对算法拆分和用傅里叶方法计算稳定性之后，可以得到稳定性条件：下面我们来验证其是否正确。4.4.1 验证稳定性条件 运行程序test6.m，可以得到如下的结果：表十和表十一是不满足稳定性条件（即）的测试例子。表十 不稳定，且满足a*a-4b0的表a=10 b=0.5 h=pi/20 t=0.2有（满足a*a-4b0）层数平均误差大小最大误差大小19.556834e-0022.136926e-

45、00121.329571e-0012.972946e-00131.651503e-0013.692793e-00142.021831e-0014.520853e-00152.470265e-0015.523560e-001表十一 不稳定，且满足a*a-4b<0的表a=0.5 b=1 h=pi/20 t=0.2有（满足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-4b0的表a=100 b=2000 h=pi/10 t=0.2有（满足a*a-4b0）层数平均误差大小最大误差大小11.851405e-