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

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

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

3、法214.2.2 a4b0的圖像方法234.2.3 結(jié)論244.3 特例（大步長(zhǎng)，高精度）254.3.1 特殊例子254.3.2 結(jié)論264.4 拋物型穩(wěn)定性條件測(cè)試264.4.1 驗(yàn)證穩(wěn)定性條件264.4.2 結(jié)論28五 總結(jié)29六 參考文獻(xiàn)30一 引言1.1 背景 微分方程是構(gòu)造力學(xué)等領(lǐng)域的數(shù)學(xué)模型的主要方法。一般來(lái)講，無(wú)論是物體運(yùn)動(dòng)軌跡或是流體速度測(cè)定模型都需要用到數(shù)值方法去求解。通過(guò)對(duì)橢圓型、拋物型和雙曲型方程的研究和探討，可以和現(xiàn)實(shí)中許多實(shí)際問(wèn)題相互掛鉤，并且得到解決問(wèn)題的方法。許多的數(shù)學(xué)家都致力于其中，如歐拉、柯西、貝努利、拉格朗日等人都為之做出了重要貢獻(xiàn)。有限差分法、有限元法等等

4、，為我們提供了求解的方法和手段。通過(guò)對(duì)微分方程數(shù)值解的學(xué)習(xí)和研究，我們可以得到許多實(shí)際問(wèn)題的求解方法。雖然大部分問(wèn)題求出來(lái)的是近似解，但是只要精度夠高，能滿足現(xiàn)實(shí)中的需求，這就足以體現(xiàn)它的重要性。 在偏微分方程的求解問(wèn)題中，Possion方程第一邊值問(wèn)題占據(jù)重要位置，五點(diǎn)差分格式等方法為我們提供了求解的方法。而在實(shí)際問(wèn)題中，重調(diào)和方程的二階邊值問(wèn)題也相當(dāng)重要。于是，我們提出如下問(wèn)題。1.2 問(wèn)題的提出重調(diào)和方程的二階邊值問(wèn)題 （1）其中是Laplace算子，是二維平面上的有限區(qū)域，是其光滑邊界，a，b非負(fù)常數(shù)。數(shù)值算列：1. ，使問(wèn)題（1）存在精確解。2. ，問(wèn)題（1）存在精確解u=ysinx

5、+xsiny分別取不同的a，b值計(jì)算，并且包括和兩種情形。并將上面的方法推廣應(yīng)用到拋物型方程 數(shù)值算列， ，使問(wèn)題（11）-（14）存在精確解。二 方法介紹由于區(qū)域的不同，會(huì)導(dǎo)致算法格式的變化，所以這里為了說(shuō)明更簡(jiǎn)單，我們固定區(qū)域?yàn)榫匦螀^(qū)域。其他的一些區(qū)域，可由矩形區(qū)域變化而來(lái)，或者進(jìn)一步討論即可得到差分格式，我們這里不再討論。2.1 五點(diǎn)差分方法考慮矩形區(qū)域=0xa，0yb上，二階線性橢圓型方程Lu=f,其中(x,y) 第一邊值問(wèn)題。假設(shè)矩形區(qū)域網(wǎng)格剖分均勻：h1=a/M，h2=b/N。于是網(wǎng)域包含(M-1)*(N-1)個(gè)內(nèi)點(diǎn)，且均為正則內(nèi)點(diǎn)。 設(shè)（i，j），并且u(x,y)充分光滑，則沿著

6、x和y方向分別用中心差商代替導(dǎo)數(shù)有 這里表示u(,)。而qu和右端項(xiàng)f直接有qu ，f= 于是方程在方程的(i,j)點(diǎn)被表示為：當(dāng)我們略去截?cái)嗾`差=，就可以得到五點(diǎn)差分格式： 其中 由方程我們可以看到，點(diǎn)(i，j)與相鄰的四個(gè)點(diǎn)（i-1,j）、（i+1,j）、（i,j+1）、（i,j-1）都有關(guān)系，于是稱之為五點(diǎn)差分格式。i-1,ji,j+1i,j-1i+1,ji,j圖一 五點(diǎn)差分格式(i,j)關(guān)系點(diǎn)圖 現(xiàn)在我們?cè)龠M(jìn)一步簡(jiǎn)化問(wèn)題，令p(x,y)=1,q(x,y)=0時(shí)，方程變?yōu)镻oisson方程 Lu 而我們的五點(diǎn)差分格式也相應(yīng)地變?yōu)?又由于我們題目中的求解域剛好是個(gè)正方形，所以我們可以把步長(zhǎng)

7、h1=h2=h，用正方形網(wǎng)格剖分區(qū)域后，五點(diǎn)差分格式簡(jiǎn)化為 方程就是我們問(wèn)題中所運(yùn)用到的最簡(jiǎn)化五點(diǎn)差分格式。下面的二階邊值問(wèn)題和推廣到拋物型方程都需要利用到它。2.2 二階邊值問(wèn)題分解 重調(diào)和方程的二階邊值問(wèn)題 （1）其中是Laplace算子，是二維平面上的有限區(qū)域，是其光滑邊界，a，b非負(fù)常數(shù)。雖然直接利用差商逼近導(dǎo)數(shù)的方法可以得到相應(yīng)的差分格式，但由于涉及的離散點(diǎn)數(shù)較多，從而邊界點(diǎn)處的計(jì)算較為困難。于是，我們采用另一種方法，適當(dāng)?shù)貙⑵浞纸獬蓛蓚€(gè)Possion方程求解。 我們使用兩種完全不同的方法去逼近問(wèn)題（1），而具體運(yùn)用哪一種方法與a4b有關(guān)。當(dāng)a4b0，我們立刻可以把問(wèn)題（1）分解成兩

8、個(gè)二階方程，然后運(yùn)用五點(diǎn)差分格式進(jìn)行求解。而當(dāng)a4b0時(shí)，我們提出了一種方法，把問(wèn)題（1）簡(jiǎn)化為順序迭代二階方程，然后逐步求解。而這個(gè)方法的收斂性也已經(jīng)被實(shí)驗(yàn)證明。2.2.1 當(dāng)a4b0時(shí)令，那么邊值問(wèn)題（1）可轉(zhuǎn)化為下列兩個(gè)邊值問(wèn)題 （2） （3）Dirichlet邊值問(wèn)題（2）及（3）可以用有限元法或者有限差分法求解，而程序中用的是五點(diǎn)差分格式逼近，從而求出問(wèn)題（1）的數(shù)值解。2.2.2 當(dāng)a4b0時(shí)令 （4）那么邊值問(wèn)題（1）可轉(zhuǎn)化為下列兩個(gè)邊值問(wèn)題 （5） （6）其中函數(shù)未知且滿足關(guān)聯(lián)式（4），所以（5）很難求解，但是若已知，那么（5）和（6）就是易求解的Poisson問(wèn)題。下面我們采

9、用迭代法計(jì)算：先估計(jì)，再求解Poisson問(wèn)題（5）和（6）得,。1.給定，如 （7） 2.已知，求解兩個(gè)Poisson問(wèn)題 （8） （9）3.計(jì)算新的迭代值 （10）其中為迭代參數(shù)。假設(shè), 空間步長(zhǎng),最優(yōu)參數(shù)其中， 方法具有二階精度。2.3 推廣到拋物型方程現(xiàn)在將上面的方法推廣應(yīng)用到拋物型方程 將方程（11）時(shí)間離散得 （15）其中，。 得到（15）之后，我們可以發(fā)現(xiàn)它與問(wèn)題（1）類似，于是我們可以運(yùn)用上面的方法求解。從底層算起，然后計(jì)算下一個(gè)時(shí)間層，如此下去，知道計(jì)算完全部的時(shí)間層。三 算法程序3.1 程序列表表一 第一、二題相關(guān)的程序列表文件名標(biāo)注main1.m第一題頭文件，輸入main

10、1(a,b,h)可得到近似解main2.m第二題頭文件，輸入main2(a,b,h)可得到近似解main3.m附加例子頭文件，輸入main3(a,b,h)可得到近似解Fivepoints1.m當(dāng)滿足條件a4b0的差分算法Fivepoints2.m當(dāng)滿足條件a4b0的差分算法f.m計(jì)算右端項(xiàng)f的函數(shù)f1.m計(jì)算邊值g1的函數(shù)f2.m計(jì)算邊值g2的函數(shù)go.m計(jì)算邊值的函數(shù)bes.m計(jì)算最優(yōu)參數(shù)Iteration.m迭代計(jì)算得到滿足精度的Correct.m計(jì)算精確解Testerr.m計(jì)算平均誤差和最大誤差表二 第三題相關(guān)的程序列表文件名標(biāo)注main4.m第三題頭文件，輸入main4(a,b,h,t

11、)可得到近似解Parabolic.m計(jì)算每一層的近似解，由底層開(kāi)始，逐步計(jì)算更高的層級(jí)reFivepoints1.m專門用于第三題，其余類似Fivepoints1.mreFivepoints2.m專門用于第三題，其余類似Fivepoints2.mIteration1.m專門用于第三題，其余類似Iteration1.mCorrect1.m專門用于第三題，其余類似Correct1.mf.m、f1.m、f2.m、go.m、bes.m、Testerr.m用途如上表表三 試驗(yàn)列表文件名標(biāo)注test1.m得到a4b0的差分算法的結(jié)果，并用第一種方法計(jì)算收斂階test2.m得到a4b0的差分算法的結(jié)果，并用

12、第一種方法計(jì)算收斂階test3.m檢驗(yàn)如何得到較高精度的結(jié)果test4.m用第二種方法計(jì)算a4b0的收斂階test5.m用第二種方法計(jì)算a4b0的收斂階test6.m對(duì)第三題中，穩(wěn)定性判斷條件h412/b的測(cè)試3.2 算法流程圖與介紹3.2.1 重調(diào)和方程二階邊值程序圖二 重調(diào)和方程二階邊值程序圖輸入a,b,h得到近似解u= A2G2,并與精確解比較，計(jì)算誤差大小否是00計(jì)算a4b 計(jì)算公式（2），得到A1,G1v=A1G1得到G1，計(jì)算v=A1G1計(jì)算公式（3），得到A2,G2令Qu1=Qu2,初始Qu1=0u=A2/G2|Qu1-Qu2|達(dá)到精度要求？由公式（6），得到A2,G2計(jì)算公式（

13、5），得到A1計(jì)算公式參考2.2.2的方法介紹，上面就是算法的流程圖。 對(duì)應(yīng)第一題、第二題和附加例子的程序。首先輸入a,b,h到相應(yīng)地題目，如main1(a,b,h)就是運(yùn)算第一題。然后進(jìn)行條件a4b的判斷，當(dāng)0時(shí)，進(jìn)入算法1（Fivepoints1）；否則，進(jìn)入算法2（Fivepoints2）。 當(dāng)a4b0，即算法1中，先對(duì)公式（2）進(jìn)行五點(diǎn)差分算法，構(gòu)造矩陣A1和右邊項(xiàng)G1，得到（2）的解v。然后對(duì)公式（3）進(jìn)行五點(diǎn)差分算法，構(gòu)造矩陣A2和右邊項(xiàng)G2（其中G2需要通過(guò)v的計(jì)算得到）。最后得到公式（3）的解u，這也是我們需要求的近似解。 當(dāng)a4b1.0e-6)&(k1000) %使精度達(dá)到要

14、求，并且迭代次數(shù)小于1000 Qu1=Qu2; U,Qu2=Iteration(g1,Qu1,A1,a1,x1,x2,y1,y2,x,y,h,ma,b,M,N,l); k=k+1;end%3.2.2 拋物型方程程序圖三 拋物型方程程序圖是否，計(jì)算下一層floor=floor+1輸入a,b,h,t得到近似解u= A2G2,并與精確解比較，計(jì)算誤差大小否是00計(jì)算a4b 計(jì)算公式（2），得到A1,G1v=A1G1得到G1，計(jì)算v=A1G1計(jì)算公式（3），得到A2,G2令Qu1=Qu2,初始Qu1=0u=A2/G2|Qu1-Qu2|達(dá)到精度要求？由公式（6），得到A2,G2計(jì)算公式（5），得到A1首

15、先計(jì)算第一層的解，floor=1T= floor*t已經(jīng)是頂層了嗎？結(jié)束 計(jì)算公式參考2.2.3的方法介紹，上面就是算法的流程圖。 由于算法和3.2.1大部分相同，所以這里主要介紹時(shí)間層的迭代問(wèn)題。輸入main4(a,b,h,t)就會(huì)進(jìn)入程序。算法開(kāi)始會(huì)根據(jù)t來(lái)決定時(shí)間層一共有多少層，設(shè)定完畢后，就進(jìn)入第一層的計(jì)算。之后就如同3.2.1計(jì)算，不過(guò)特別的是，右端項(xiàng)G1與上一層的近似解u有關(guān)系（當(dāng)計(jì)算第一層時(shí)，G1可由初值得到）。 得到第一層的解u1之后，令floor=floor+1，進(jìn)入第二層的計(jì)算。如此迭代下去，得到全部層數(shù)的近似解。下面給出時(shí)間層的迭代算法：%if a*a-4*b0 for

16、floor=1:L %表示第幾層 T1=floor*t; %時(shí)間平面參數(shù) Gu2=reFivepoints2(x1,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); %計(jì)算floor層的誤差 endelse for floor=1:L %表示第幾層 T1=floor*t; %時(shí)間平面參數(shù) Gu2=reFivepoints1(x1,x2,y1,y2,T1,a,b,h,t,Gu1)

17、; Gu1=Gu2; Co=Correct1(x1,x2,y1,y2,h,T1); %第幾層的精確解 str2=sprintf(層數(shù)為:%d,floor); disp(str2) Testerr(Gu1,Co,M,N); endend四 結(jié)果與分析4.1 誤差對(duì)比方法判斷重調(diào)和方程算法的收斂階 由其他科學(xué)家的研究中，我們可以知道文章中兩種逼近重調(diào)和方程的方法收斂階均為。現(xiàn)在，我們來(lái)證明這個(gè)結(jié)論的正確性。因?yàn)榉椒ㄊ諗侩A均為，所以我們可以令a,b相同，h1=2*h2。假如輸入h1的誤差是輸入h2的誤差的四倍，那么就驗(yàn)證了我們想法。 不過(guò)這種判斷方法有個(gè)相當(dāng)不利的地方，就是需要步長(zhǎng)足夠小，這樣才能使

18、計(jì)算出來(lái)的數(shù)據(jù)較為精確。由于一般家用電腦把網(wǎng)格點(diǎn)取成60*60就已經(jīng)要計(jì)算相當(dāng)長(zhǎng)的時(shí)間，所以這里只有把別有的結(jié)果展示出來(lái)。由于有兩種方法，一種a4b0，另外一種a4b0，我們分開(kāi)討論。4.1.1 a4b0的誤差對(duì)比方法 當(dāng)運(yùn)行程序test1.m之后，就可以得到如下表格。表四 第一題，b=0、0.25、0.5的誤差表格a4b0 第一題 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.349345

19、e-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.4203714.204728a4b0 第一題 b=0.25ah=pi/10的誤差E1h=pi/2

20、0的誤差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-0

21、044.4199114.2046204.05.315428e-0031.202688e-0032.860443e-0044.4196244.204552a4b0 第一題 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.419

22、5874.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)格點(diǎn)較小的兩次比較（E1/E2）結(jié)果接近4.40，而網(wǎng)格點(diǎn)較大的兩次比較（E2/E3）結(jié)果接近4.20。這說(shuō)明網(wǎng)格點(diǎn)越小，結(jié)果會(huì)更加靠近4.00。表五 第二題，b=0、0.25、0.5的誤差表格a4

23、b0 第二題 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.59.741373

24、e-0032.222717e-0035.297593e-0044.3826424.1957114.09.617922e-0032.194722e-0035.230974e-0044.3822974.195627a4b0 第二題 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.79.736602e-0032

25、.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.195405a4b0 第二題 b=0.5ah=pi/10的誤差E1h=pi/20的誤差E2h=pi/40的誤差E3E1/E2E2/E32.09.695610e-0032.212894

26、e-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-0032.108570e-0035.02

27、6162e-0044.380470 4.195190分析：從表四和表五我們可以看到誤差雖然和a,b有關(guān)系，但是影響誤差最大的因素還是步長(zhǎng)。當(dāng)步長(zhǎng)減小為一半，誤差可以變?yōu)?/4。這說(shuō)明要有高精度的結(jié)果，要求步長(zhǎng)足夠小。 由于表中E1/E2約為4.40，E2/E3約為4.20，和我們的目標(biāo)結(jié)果4.00相比，顯然不太理想。這是由于步長(zhǎng)太小，使算法的其他影響因素過(guò)大，導(dǎo)致結(jié)果產(chǎn)生偏差。但是如同前文所說(shuō)的，當(dāng)網(wǎng)格點(diǎn)過(guò)多，家用電腦計(jì)算速度無(wú)法達(dá)到要求。所以在a4b0的方法中，將會(huì)引用其他論文中計(jì)算的結(jié)果。4.1.2 a4b0的誤差對(duì)比方法 下面給出其他論文中計(jì)算出來(lái)的結(jié)果。當(dāng)然也可以運(yùn)行文件中的test2

28、.m程序，可是由于步長(zhǎng)過(guò)大，不建議使用，所以不再列出，但讀者可以自己測(cè)試一下。表六 第一題，a=0、0.5、1的誤差表格a4b0 第一題 a=0b網(wǎng)格點(diǎn)為65*65網(wǎng)格點(diǎn)為129*129網(wǎng)格點(diǎn)為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.0550e53.98293.93151.582.9118e487.2101e581.7337e54.0385

29、4.15882.092.6904e496.8215e591.8013e53.94403.78702.5102.4561e4106.0247e5101.3909e54.07674.3315a4b0 第一題 a=0.5b網(wǎng)格點(diǎn)為65*65網(wǎng)格點(diǎn)為129*129網(wǎng)格點(diǎn)為257*257E1/E2E2/E3迭代k次誤差E1迭代k次誤差E2迭代k次誤差E30.353.7378e458.5291e552.1354e53.99903.99410.763.4125e467.9097e561.9648e54.00694.02571.073.2179e477.5415e571.8938e53.99593.98221

30、.582.9118e486.9316e581.7184e54.00864.03382.092.6904e496.4782e591.6374e53.98923.95642.592.4561e496.1521e591.6340e53.93783.7651a4b0 第一題 a=1b網(wǎng)格點(diǎn)為65*65網(wǎng)格點(diǎn)為129*129網(wǎng)格點(diǎn)為257*257E1/E2E2/E3迭代k次誤差E1迭代k次誤差E2迭代k次誤差E30.353.1879e457.9702e551.9938e53.99983.99750.762.9968e467.4868e561.8672e54.00284.00961.062.8644e46

31、7.1267e561.7480e54.01934.07711.572.6823e476.7405e571.7204e53.97943.91802.082.5061e486.2346e581.5283e54.01974.07942.592.3658e495.9385e591.5090e53.98383.9354分析：從表六可以看出計(jì)算結(jié)果與目標(biāo)結(jié)果相當(dāng)接近。表七 第二題，a=0、0.5、1的誤差表格a4b0 第二題 a=0b網(wǎng)格點(diǎn)為65*65網(wǎng)格點(diǎn)為129*129網(wǎng)格點(diǎn)為257*257E1/E2E2/E3迭代k次誤差E1迭代k次誤差E2迭代k次誤差E30.365.4815e461.3703e46

32、3.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.4078e53.95263.8268a4b0 第二題 a=0.5b網(wǎng)格點(diǎn)為65*65網(wǎng)格點(diǎn)為129*129網(wǎng)格點(diǎn)為257*257E1/E2E2/E3迭代k次誤差E1

33、迭代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.5458e53.97203.89262.5103.6555e4109.0709e5102.1945e54.02994.1335a4b0 第二題 a=1b網(wǎng)格點(diǎn)為65*65網(wǎng)格點(diǎn)為129

34、*129網(wǎng)格點(diǎn)為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.5970e54.00544.02352.093.9278e499.8369e592.4760e53.99293.97292.593.7068e499.3636e592.4338e53.95873.847

35、3分析：從表六、表七可以看到E1/E2和E2/E3確實(shí)在4.00附近變動(dòng)，這證明了我們猜想方法收斂階為是正確的。而且還可以看到，1.當(dāng)a，b固定的時(shí)候，迭代次數(shù)k與我們的網(wǎng)格點(diǎn)數(shù)無(wú)關(guān)，說(shuō)明步長(zhǎng)對(duì)迭代次數(shù)k無(wú)影響；2.當(dāng)a和網(wǎng)格點(diǎn)數(shù)固定，即a、h固定時(shí)，b越大，迭代次數(shù)k也越大。4.1.3 結(jié)論 從4.1的討論當(dāng)中我們可以知道一下幾個(gè)結(jié)論：1.步長(zhǎng)的大小對(duì)誤差影響十分大。步長(zhǎng)越小，方法精度越高；2.兩種方法（a4b0和a4b0）收斂階為；3.誤差對(duì)比方法確實(shí)可以得知收斂階，但是前提是步長(zhǎng)足夠小，這就需要有較好的設(shè)備才能得出結(jié)論；4.對(duì)于a4b0的算法，我們發(fā)現(xiàn).當(dāng)a，b固定的時(shí)候，步長(zhǎng)對(duì)迭代次數(shù)

36、k無(wú)影響；.當(dāng)a、h固定時(shí)，b越大，迭代次數(shù)k也越大。所以在滿足a4b0時(shí)，我們可以適當(dāng)減小b使迭代次數(shù)減小，更快得到近似解。4.2 圖像方法判斷重調(diào)和方程算法的收斂階 現(xiàn)在介紹另外一種方法判斷算法的收斂階，這也是一種相當(dāng)常用的方法來(lái)判斷收斂階。假設(shè)誤差e=， 我們可令e=ch，其中c為常數(shù) 然后對(duì)上式兩邊自然對(duì)數(shù)，有l(wèi)ne = lnc + 2lnh這條公式說(shuō)明了算法的誤差對(duì)數(shù)與步長(zhǎng)對(duì)數(shù)呈線性相關(guān)。如果我們計(jì)算出來(lái)的誤差對(duì)數(shù)線與y=2lnh相互平行，則證明了算法收斂階為e=。同樣的，需要對(duì)a4b0和a4b0兩種算法分析。4.2.1 a4b0的圖像方法 運(yùn)行程序test4.m即可得到下面兩個(gè)圖像

37、：圖四 第一題，a4b0收斂階判斷圖分析：可以看到兩條曲線是互相平行的。圖五 第二題，a4b0收斂階判斷圖分析：圖三和圖四中的兩條曲線都是互相平行的，說(shuō)明了a4b0收斂階確實(shí)是。4.2.2 a4b0的圖像方法運(yùn)行程序test5.m，有如下圖像：圖六 第一題，a4b0收斂階判斷圖分析：在a4b0的方法中，同樣看到兩條曲線平行。圖七 第二題，a4b0收斂階判斷圖分析：由圖六和圖七可以看出，a4b0方法的收斂階也是。4.2.3 結(jié)論 由上面的分析，可以得知：1.a4b0和a4b0兩種算法收斂階是；2.如果覺(jué)得用看圖來(lái)判斷兩條直線是否平行不太嚴(yán)密的話，可以計(jì)算兩條直線斜率是否相等。如果相等，則證明結(jié)論

38、無(wú)誤；否則說(shuō)明收斂階不是。不過(guò)這樣的話，又會(huì)和4.1的方法一樣，出現(xiàn)因?yàn)椴介L(zhǎng)不夠小而導(dǎo)致結(jié)論不正確的情況，所以這里就不給出來(lái)了。4.3 特例（大步長(zhǎng)，高精度）4.3.1 特殊例子在某一些特殊例子中，盡管步長(zhǎng)較大，我們?nèi)匀豢梢匀〉酶呔鹊慕Y(jié)果。例如精確解為，在-2,2*-2,2的區(qū)域里面，而右端項(xiàng)f=8-2a(x+y-8)+b(x-4)(y-4)。這個(gè)時(shí)候我們能取到相當(dāng)高精度的解，如表八和表九（運(yùn)行test3.m可得到）。表八 a4b0，特殊例子的誤差表格網(wǎng)格數(shù)為10*10ab=0.0時(shí)的誤差b=0.25時(shí)的誤差b=0.5時(shí)的誤差2.02.116277e-0157.779785e-0151.98

39、9629e-0142.32.812565e-0157.352144e-0151.744010e-0142.72.050486e-0152.675500e-0152.532953e-0153.02.412336e-0155.104285e-0153.234724e-015表九 a4b0，特殊例子的誤差表格網(wǎng)格數(shù)為10*10b迭代次數(shù)ka=0.0時(shí)的誤差迭代次數(shù)ka=0.5時(shí)的誤差迭代次數(shù)ka=1.0時(shí)的誤差0.3122.253198e-008101.825354e-00891.282883e-0080.7243.223164e-008172.416017e-008142.030519e-0081

40、.0423.718790e-008251.455840e-008191.322050e-0081.53013.218638e-008482.140648e-008301.504249e-008分析：上面兩個(gè)表都是在網(wǎng)格點(diǎn)數(shù)為10*10的時(shí)候得到的，但是精度都分別達(dá)到了和，與前面兩題得出的結(jié)果相比準(zhǔn)確很多（因?yàn)榍懊娴膬深}盡管網(wǎng)格數(shù)為65*65，精度也只達(dá)到）。這是因?yàn)閷?duì)于任何步長(zhǎng)大小，二次函數(shù)中心差分格式的逼近誤差為零，所以使近似解有相當(dāng)高的精度。 對(duì)于這個(gè)例子，下面給出一幅在a=0.5，b=1,網(wǎng)格點(diǎn)為20*20的函數(shù)圖，圖八?？梢砸?jiàn)到x和y越趨向0，函數(shù)值就越大。圖八 a=0.5，b=1,h

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

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