重調(diào)和方程二階邊值問題求解

1、微分方程數(shù)值解課程設(shè)計微分方程數(shù)值解課程設(shè)計重調(diào)和方程的二階邊值問題的求解學(xué)院、系： 理學(xué)院 數(shù)學(xué)系 專 業(yè)： 信息與計算科學(xué) 姓 名： 陳劍宇 學(xué) 號： 201030470140 任課教師： 黃鳳輝 提交日期： 2012/12/27 總評成績： 摘 要本次課程設(shè)計，主要討論重調(diào)和方程二階邊值問題的求解。文章分成四部分：第一部分，介紹如何將重調(diào)和方程的二階邊值問題分解，進而用經(jīng)典的五點差分格式逼近；第二部分，列出對應(yīng)的MATLAB算法和流程圖；第三部分，通過上述方法求解具體的問題，并分析方法的精度和收斂階；第四部分，將方法推廣應(yīng)用到拋物型方程，然后分析穩(wěn)定性。關(guān)鍵詞：重調(diào)和方程的二階邊值問題；

2、五點差分格式；MATLAB；拋物型方程目 錄一 引言41.1 背景41.2 問題的提出5二 方法介紹62.1 五點差分方法62.2 二階邊值問題分解82.2.1 當(dāng)a²4b0時82.2.2 當(dāng)a²4b<0時92.3 推廣到拋物型方程10三 算法程序113.1 程序列表113.2 算法流程圖與介紹123.2.1 重調(diào)和方程二階邊值程序123.2.2 拋物型方程程序14四 結(jié)果與分析164.1 誤差對比方法判斷重調(diào)和方程算法的收斂階164.1.1 a²4b0的誤差對比方法164.1.2 a²4b<0的誤差對比方法184.1.3 結(jié)論204.2 圖

3、像方法判斷重調(diào)和方程算法的收斂階214.2.1 a²4b0的圖像方法214.2.2 a²4b<0的圖像方法234.2.3 結(jié)論244.3 特例（大步長，高精度）254.3.1 特殊例子254.3.2 結(jié)論264.4 拋物型穩(wěn)定性條件測試264.4.1 驗證穩(wěn)定性條件264.4.2 結(jié)論28五 總結(jié)29六 參考文獻30一 引言1.1 背景 微分方程是構(gòu)造力學(xué)等領(lǐng)域的數(shù)學(xué)模型的主要方法。一般來講，無論是物體運動軌跡或是流體速度測定模型都需要用到數(shù)值方法去求解。通過對橢圓型、拋物型和雙曲型方程的研究和探討，可以和現(xiàn)實中許多實際問題相互掛鉤，并且得到解決問題的方法。許多的數(shù)學(xué)

4、家都致力于其中，如歐拉、柯西、貝努利、拉格朗日等人都為之做出了重要貢獻。有限差分法、有限元法等等，為我們提供了求解的方法和手段。通過對微分方程數(shù)值解的學(xué)習(xí)和研究，我們可以得到許多實際問題的求解方法。雖然大部分問題求出來的是近似解，但是只要精度夠高，能滿足現(xiàn)實中的需求，這就足以體現(xiàn)它的重要性。 在偏微分方程的求解問題中，Possion方程第一邊值問題占據(jù)重要位置，五點差分格式等方法為我們提供了求解的方法。而在實際問題中，重調(diào)和方程的二階邊值問題也相當(dāng)重要。于是，我們提出如下問題。1.2 問題的提出重調(diào)和方程的二階邊值問題 （1）其中是Laplace算子，是二維平面上的有限區(qū)域，是其光滑邊界，a，

5、b非負常數(shù)。數(shù)值算列：1. ，使問題（1）存在精確解。2. ，問題（1）存在精確解u=ysinx+xsiny分別取不同的a，b值計算，并且包括和兩種情形。并將上面的方法推廣應(yīng)用到拋物型方程 數(shù)值算列， ，使問題（11）-（14）存在精確解。二 方法介紹由于區(qū)域的不同，會導(dǎo)致算法格式的變化，所以這里為了說明更簡單，我們固定區(qū)域為矩形區(qū)域。其他的一些區(qū)域，可由矩形區(qū)域變化而來，或者進一步討論即可得到差分格式，我們這里不再討論。2.1 五點差分方法考慮矩形區(qū)域=0<x<a，0<y<b上，二階線性橢圓型方程Lu=f,其中(x,y) 第一邊值問題。假設(shè)矩形區(qū)域網(wǎng)格剖分均勻：h1=

6、a/M，h2=b/N。于是網(wǎng)域包含(M-1)*(N-1)個內(nèi)點，且均為正則內(nèi)點。 設(shè)（i，j），并且u(x,y)充分光滑，則沿著x和y方向分別用中心差商代替導(dǎo)數(shù)有 這里表示u(,)。而qu和右端項f直接有qu ，f= 于是方程在方程的(i,j)點被表示為：當(dāng)我們略去截斷誤差=，就可以得到五點差分格式： 其中 由方程我們可以看到，點(i，j)與相鄰的四個點（i-1,j）、（i+1,j）、（i,j+1）、（i,j-1）都有關(guān)系，于是稱之為五點差分格式。i-1,ji,j+1i,j-1i+1,ji,j圖一 五點差分格式(i,j)關(guān)系點圖 現(xiàn)在我們再進一步簡化問題，令p(x,y)=1,q(x,y)=0時

7、，方程變?yōu)镻oisson方程 Lu 而我們的五點差分格式也相應(yīng)地變?yōu)?又由于我們題目中的求解域剛好是個正方形，所以我們可以把步長h1=h2=h，用正方形網(wǎng)格剖分區(qū)域后，五點差分格式簡化為 方程就是我們問題中所運用到的最簡化五點差分格式。下面的二階邊值問題和推廣到拋物型方程都需要利用到它。2.2 二階邊值問題分解 重調(diào)和方程的二階邊值問題 （1）其中是Laplace算子，是二維平面上的有限區(qū)域，是其光滑邊界，a，b非負常數(shù)。雖然直接利用差商逼近導(dǎo)數(shù)的方法可以得到相應(yīng)的差分格式，但由于涉及的離散點數(shù)較多，從而邊界點處的計算較為困難。于是，我們采用另一種方法，適當(dāng)?shù)貙⑵浞纸獬蓛蓚€Possion方程求

8、解。 我們使用兩種完全不同的方法去逼近問題（1），而具體運用哪一種方法與a²4b有關(guān)。當(dāng)a²4b0，我們立刻可以把問題（1）分解成兩個二階方程，然后運用五點差分格式進行求解。而當(dāng)a²4b<0時，我們提出了一種方法，把問題（1）簡化為順序迭代二階方程，然后逐步求解。而這個方法的收斂性也已經(jīng)被實驗證明。2.2.1 當(dāng)a²4b0時令，那么邊值問題（1）可轉(zhuǎn)化為下列兩個邊值問題 （2） （3）Dirichlet邊值問題（2）及（3）可以用有限元法或者有限差分法求解，而程序中用的是五點差分格式逼近，從而求出問題（1）的數(shù)值解。2.2.2 當(dāng)a²4b

9、<0時令 （4）那么邊值問題（1）可轉(zhuǎn)化為下列兩個邊值問題 （5） （6）其中函數(shù)未知且滿足關(guān)聯(lián)式（4），所以（5）很難求解，但是若已知，那么（5）和（6）就是易求解的Poisson問題。下面我們采用迭代法計算：先估計，再求解Poisson問題（5）和（6）得,。1.給定，如 （7） 2.已知，求解兩個Poisson問題 （8） （9）3.計算新的迭代值 （10）其中為迭代參數(shù)。假設(shè), 空間步長,最優(yōu)參數(shù)其中， 方法具有二階精度。2.3 推廣到拋物型方程現(xiàn)在將上面的方法推廣應(yīng)用到拋物型方程 將方程（11）時間離散得 （15）其中，。 得到（15）之后，我們可以發(fā)現(xiàn)它與問題（1）類似，于是

10、我們可以運用上面的方法求解。從底層算起，然后計算下一個時間層，如此下去，知道計算完全部的時間層。三 算法程序3.1 程序列表表一 第一、二題相關(guān)的程序列表文件名標注main1.m第一題頭文件，輸入main1(a,b,h)可得到近似解main2.m第二題頭文件，輸入main2(a,b,h)可得到近似解main3.m附加例子頭文件，輸入main3(a,b,h)可得到近似解Fivepoints1.m當(dāng)滿足條件a²4b0的差分算法Fivepoints2.m當(dāng)滿足條件a²4b<0的差分算法f.m計算右端項f的函數(shù)f1.m計算邊值g1的函數(shù)f2.m計算邊值g2的函數(shù)go.m計算邊

11、值的函數(shù)bes.m計算最優(yōu)參數(shù)Iteration.m迭代計算得到滿足精度的Correct.m計算精確解Testerr.m計算平均誤差和最大誤差表二 第三題相關(guān)的程序列表文件名標注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的差分算法的結(jié)果，并用第一種方法計算收斂階test2.m得到a²4b<0的差分算法的結(jié)果，并用第一種方法計算收斂階test3.m檢驗如何得到較高精度的結(jié)果test4.m用第二種方法計算a²4b0的收斂階test5.m用第二種方法計算a²4b<0的收斂階test6.m對第三題中，穩(wěn)定性判斷條件h412/b的測試3.2 算法流程圖與介紹3.2.1 重調(diào)和方程二階邊值程序圖二 重調(diào)和方程二階邊值程序圖輸

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的方法介紹，上面就是算法的流程圖。 對應(yīng)第一題、第二題和附加例子的程序。首先輸入a,b,h到相應(yīng)地題目，如main1(a,b,h)就是運算第一題。然后進行條件a²4b的判斷，當(dāng)0時，進入算法1（Fivepoints1）；否則，進入算法

14、2（Fivepoints2）。 當(dāng)a²4b0，即算法1中，先對公式（2）進行五點差分算法，構(gòu)造矩陣A1和右邊項G1，得到（2）的解v。然后對公式（3）進行五點差分算法，構(gòu)造矩陣A2和右邊項G2（其中G2需要通過v的計算得到）。最后得到公式（3）的解u，這也是我們需要求的近似解。 當(dāng)a²4b<0，即算法2中，類似地對公式（5）進行五點差分算法，構(gòu)造矩陣A1和右邊項G1。但是由于右邊項G1與有關(guān)，無法直接得到，所以我們運用迭代的方法。先令，然后得到G1，然后計算出公式（5）的解v。再根據(jù)公式（6）進行五點差分，得到（6）的解u。計算到這里，還要計算（與u和有關(guān)），假若|-

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; %迭代次數(shù)while (norm(Qu2-Qu1)>1.0e-6)&&(k<1000) %使精度達到要求，并且迭代次數(shù)小于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已經(jīng)是頂層了嗎？結(jié)束 計算公式參考2.2.3的方法介紹，上面就是算法的流

17、程圖。 由于算法和3.2.1大部分相同，所以這里主要介紹時間層的迭代問題。輸入main4(a,b,h,t)就會進入程序。算法開始會根據(jù)t來決定時間層一共有多少層，設(shè)定完畢后，就進入第一層的計算。之后就如同3.2.1計算，不過特別的是，右端項G1與上一層的近似解u有關(guān)系（當(dāng)計算第一層時，G1可由初值得到）。 得到第一層的解u1之后，令floor=floor+1，進入第二層的計算。如此迭代下去，得到全部層數(shù)的近似解。下面給出時間層的迭代算法：%if a*a-4*b<0 for floor=1:L %表示第幾層 T1=floor*t; %時間平面參數(shù) 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('層數(shù)為:%d',floor); disp(str2) Testerr(Gu1,Co,M,N); %計算floor層的誤差 endelse for floor=1:L %表示第幾層 T1=floor*t; %時間平面參數(shù) 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('層數(shù)為:%d',floor); disp(str2) Testerr(Gu1,Co,M,N); endend四 結(jié)果與分析4.1 誤差對比方法判斷重調(diào)和方程算法的收斂階 由其他科學(xué)家的研究中，我們可以知道文章中兩種逼近重調(diào)和方程的方法收斂階均為?，F(xiàn)在，我們來證明這個結(jié)論的正確性。因為方法收斂階均為，所以我們可以令a,b相同，h1=2*h2。假如輸入h1的誤差是輸入h2的誤差的四倍，那么就驗證了我們想法。 不過這種判斷方法有個相當(dāng)不利的地方，就是需要步長足夠小，這樣才能使計算出來的數(shù)據(jù)較為精確。由于一般家用電腦把網(wǎng)格點取成60*60就已經(jīng)要計算相

20、當(dāng)長的時間，所以這里只有把別有的結(jié)果展示出來。由于有兩種方法，一種a²4b0，另外一種a²4b<0，我們分開討論。4.1.1 a²4b0的誤差對比方法 當(dāng)運行程序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分析：從表四最后兩列可以看出，網(wǎng)格點較小的兩次比較（E1/E2）結(jié)果接近4.40，而網(wǎng)格點較大的兩次比較（E2/E3）結(jié)果接近4.20。這說明網(wǎng)格點越小，結(jié)果會更加靠近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有關(guān)系，但是影響誤差最大的因素還是步長。當(dāng)步長減小為一半，誤差可以變?yōu)?/4。這說明要有高精度的結(jié)果，要求步長足夠小。 由于表中E1/E2約為4.40，E2/E3約為4.20，和我們的目標結(jié)果4.00相比，顯然不太理想。這是由于步長太小，使算法的其他影響因素過大，導(dǎo)致結(jié)果產(chǎn)生偏差。但是如同前文所說的，當(dāng)網(wǎng)格點過多，家用電腦計算速度無法達到要求。所以在a²4b<0的方法中，將會引用其他論文中計算的結(jié)果。4.1.2 a²4b<0的

30、誤差對比方法 下面給出其他論文中計算出來的結(jié)果。當(dāng)然也可以運行文件中的test2.m程序，可是由于步長過大，不建議使用，所以不再列出，但讀者可以自己測試一下。表六 第一題，a=0、0.5、1的誤差表格a²4b<0 第一題 a=0b網(wǎng)格點為65*65網(wǎng)格點為129*129網(wǎng)格點為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網(wǎng)格點為65*65網(wǎng)格點為129*129網(wǎng)格點為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網(wǎng)格點為65*65網(wǎng)格點為129*129網(wǎng)格點為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分析：從表六可以看出計算結(jié)果與目標結(jié)果相當(dāng)接近。表七 第二題，a=0、0.5、1的誤差表格a²4b<0 第二題 a=0b網(wǎng)格點為65*

34、65網(wǎng)格點為129*129網(wǎng)格點為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網(wǎng)格點為65*65網(wǎng)格點為129*129網(wǎng)格點為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網(wǎng)格點為65*65網(wǎng)格點為129*129網(wǎng)格點為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.當(dāng)a，b固定的時候，迭代次數(shù)k與我們的網(wǎng)格點數(shù)無關(guān)，說明步長對迭代次數(shù)k無影響；2.當(dāng)a和網(wǎng)格點數(shù)固定，即a、h固定時，b越大，迭代次數(shù)k也越大。4.1.3 結(jié)論 從4.1的討論當(dāng)中我們可以知道一下幾個結(jié)論：1.步長的大小對誤差影響十分大。步長越小，方法精度越高；2.

38、兩種方法（a²4b0和a²4b<0）收斂階為；3.誤差對比方法確實可以得知收斂階，但是前提是步長足夠小，這就需要有較好的設(shè)備才能得出結(jié)論；4.對于a²4b<0的算法，我們發(fā)現(xiàn).當(dāng)a，b固定的時候，步長對迭代次數(shù)k無影響；.當(dāng)a、h固定時，b越大，迭代次數(shù)k也越大。所以在滿足a²4b<0時，我們可以適當(dāng)減小b使迭代次數(shù)減小，更快得到近似解。4.2 圖像方法判斷重調(diào)和方程算法的收斂階 現(xiàn)在介紹另外一種方法判斷算法的收斂階，這也是一種相當(dāng)常用的方法來判斷收斂階。假設(shè)誤差e=， 我們可令e=ch²，其中c為常數(shù) 然后對上式兩邊自然對數(shù)

39、，有l(wèi)ne = lnc + 2lnh這條公式說明了算法的誤差對數(shù)與步長對數(shù)呈線性相關(guān)。如果我們計算出來的誤差對數(shù)線與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 結(jié)論 由上面的分析，可以得知：1.a²4b0和a²4b<0兩種算法收斂階是；2.如果覺得用看圖來判斷兩條直線是否平行不太嚴密的話，可以計算兩條直線斜率是否相等。如果相等，則證明結(jié)論無誤；否則說明收斂階不是。不過這樣的話，又會和4.1的方法一樣，出現(xiàn)因為步長不夠小而導(dǎo)致結(jié)論不正確的情況

41、，所以這里就不給出來了。4.3 特例（大步長，高精度）4.3.1 特殊例子在某一些特殊例子中，盡管步長較大，我們?nèi)匀豢梢匀〉酶呔鹊慕Y(jié)果。例如精確解為，在-2,2*-2,2的區(qū)域里面，而右端項f=8-2a(x²+y²-8)+b(x²-4)(y²-4)。這個時候我們能取到相當(dāng)高精度的解，如表八和表九（運行test3.m可得到）。表八 a²4b0，特殊例子的誤差表格網(wǎng)格數(shù)為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，特殊例子的誤差表格網(wǎng)格數(shù)為10*10b迭代次數(shù)ka=0.0時的誤差迭代次數(shù)ka=0.5時的誤差迭代次數(shù)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分析：上面兩個表都是在網(wǎng)格點數(shù)為10*10的時候得到的，但是精度都分別達到了和，與前面兩題得出的結(jié)果相比準確很多（因為前面的兩題盡管網(wǎng)格數(shù)為65*65，精度也只達到）。這是因為對于任何步長大小，二次函數(shù)中心差分格式的逼近誤差為零，所以使近似解有相當(dāng)高的精度。 對于這個例子，下面給出一幅在a=0.5，b=1,網(wǎng)格點為20*20的函數(shù)圖，圖八?？梢砸姷絰和y越趨向0，函數(shù)值就越大。圖八 a=0.5，b=1,h=0.2的函數(shù)圖4.3.

44、2 結(jié)論 由上面可以知道：在某一些特例之中，我們的兩種算法可以在步長較小的情況下，得到高精度的解。這些特例指的是二次函數(shù)。4.4 拋物型穩(wěn)定性條件測試 當(dāng)推廣到拋物型方程的時候，我們需要對其穩(wěn)定性進行分析。在對算法拆分和用傅里葉方法計算穩(wěn)定性之后，可以得到穩(wěn)定性條件：下面我們來驗證其是否正確。4.4.1 驗證穩(wěn)定性條件 運行程序test6.m，可以得到如下的結(jié)果：表十和表十一是不滿足穩(wěn)定性條件（即）的測試例子。表十 不穩(wěn)定，且滿足a*a-4b0的表a=10 b=0.5 h=pi/20 t=0.2有（滿足a*a-4b0）層數(shù)平均誤差大小最大誤差大小19.556834e-0022.136926e-

45、00121.329571e-0012.972946e-00131.651503e-0013.692793e-00142.021831e-0014.520853e-00152.470265e-0015.523560e-001表十一 不穩(wěn)定，且滿足a*a-4b<0的表a=0.5 b=1 h=pi/20 t=0.2有（滿足a*a-4b<0）層數(shù)平均誤差大小最大誤差大小12.542859e-0015.685882e-00124.263645e-0019.533591e-00135.734781e-0011.282308e+00047.244495e-0011.619883e+00058.957729e-0012.002965e+000分析：可見數(shù)據(jù)誤差大，而且每下一層，誤差增速更快，兩層之間的誤差增長較大。這說明了初始層產(chǎn)生的誤差不能有效地控制以后每層產(chǎn)生的誤差，與穩(wěn)定含義不符。下面的表十二和表十三測試例子，滿足了穩(wěn)定性條件（）。表十二 穩(wěn)定，且滿足a*a-4b0的表a=100 b=2000 h=pi/10 t=0.2有（滿足a*a-4b0）層數(shù)平均誤差大小最大誤差大小11.851405e-