計算方法第八章

計算方法第八章第一頁，共五十八頁，編輯于2023年，星期五本章主要研究常微分方程初值問題的數(shù)值求解：通常，假設(shè)函數(shù)f關(guān)于第二個變量滿足李普希茨條件（L條件），即為存在常數(shù)L>0，使得第二頁，共五十八頁，編輯于2023年，星期五第一節(jié)一般概念1.1歐拉法及其簡單改進(jìn)方法：選擇適當(dāng)?shù)墓?jié)點，用差分近似微分，將方程離散化，從而求在這些節(jié)點上的函數(shù)值。歐拉方法第三頁，共五十八頁，編輯于2023年，星期五例子：下面我們分別取步長為0.1與0.01進(jìn)行計算，計算結(jié)果顯示在下面的圖中。第四頁，共五十八頁，編輯于2023年，星期五步長為0.1的計算結(jié)果。第五頁，共五十八頁，編輯于2023年，星期五步長為0.01的計算結(jié)果第六頁，共五十八頁，編輯于2023年，星期五0.01 0.99005 0.990.1 0.90484 0.904380.2 0.81873 0.817910.3 0.74082 0.73970.40.67032 0.668970.41 0.66365 0.662280.59 0.55433 0.552680.6 0.54881 0.547160.90.40657 0.404730.91 0.40252 0.400680.99 0.37158 0.3697310.36788 0.36603第七頁，共五十八頁，編輯于2023年，星期五

DOUBLEPRECISIONh,y(0：100)OPEN(20,FILE='OUTPUT.DAT',STATUS="UNKNOWN") h=1.0/100 y(0)=1.0 do10i=1,100 y(i)=y(i-1)*(1.0-h) write(20,*)i*h,y(i)10continue END第八頁，共五十八頁，編輯于2023年，星期五精度更高的歐拉公式：方法：選擇計算導(dǎo)數(shù)值的精度更好的差分格式。歐拉中點公式第九頁，共五十八頁，編輯于2023年，星期五利用中點公式求解微分方程時，有一個問題，就是計算時需要兩個迭代初值！對于這個問題，我們可以先用歐拉公式，通過給定的初值計算出第一個點的值，然后在利用這兩個（第一和第二個點）的值進(jìn)行計算，直到計算出全部節(jié)點上的值。下面，我們用中點公式求解剛才的例子，計算的步長取0.01，可以看到，計算的精度比較高。此時，計算公式為第十頁，共五十八頁，編輯于2023年，星期五計算結(jié)果第十一頁，共五十八頁，編輯于2023年，星期五0.02 0.9802 0.98020.1 0.90484 0.904840.2 0.81873 0.818740.3 0.74082 0.740840.4 0.67032 0.670350.5 0.60653 0.606560.6 0.54881 0.548850.7 0.49659 0.496630.8 0.44933 0.449380.9 0.40657 0.40663第十二頁，共五十八頁，編輯于2023年，星期五1.2歐拉方法的其他改進(jìn)微分方程數(shù)值解的關(guān)鍵在于對導(dǎo)數(shù)的處理，可以用差分來近似導(dǎo)數(shù)，也可以通過積分，將導(dǎo)數(shù)項化掉。對于方程：首先，作出劃分設(shè)已經(jīng)求出第n個節(jié)點的函數(shù)值，在區(qū)間上對方程兩邊積分容易看出，要求第n+1個節(jié)點的函數(shù)值，關(guān)鍵在于選擇適當(dāng)?shù)姆e分公式計算積分！第十三頁，共五十八頁，編輯于2023年，星期五（1）如選擇下矩形公式，則這正是前面的歐拉公式。（2）如選擇上矩形公式，則這是所謂的后退歐拉公式。（3）如選擇梯形公式，則這是所謂的歐拉梯形公式。第十四頁，共五十八頁，編輯于2023年，星期五直接利用已經(jīng)求得的已知節(jié)點上的值計算未知節(jié)點上的函數(shù)值的算法稱為顯式法。例如：歐拉公式、歐拉中點公式計算未知節(jié)點上的函數(shù)值時，用到了未知節(jié)點上的函數(shù)值，這種算法稱為隱式法。例如：后退歐拉法、歐拉梯形公式顯然，利用隱式法求微分方程的數(shù)值解是，需要從表達(dá)式中反解未知節(jié)點上的函數(shù)值。第十五頁，共五十八頁，編輯于2023年，星期五1.3隱式法的具體計算：例如歐拉梯形公式用迭代法計算n+1步的值。（1）簡單迭代收斂的條件：第十六頁，共五十八頁，編輯于2023年，星期五（2）牛頓迭代必須指出，在真正計算中，常用的是簡單迭代，而且只迭代一步，迭代初值稱為預(yù)測，迭代稱為校正，這樣組成的一組計算公式稱為預(yù)測－－校正公式。第十七頁，共五十八頁，編輯于2023年，星期五預(yù)測－校正公式也稱為改進(jìn)的歐拉法，將上面的組合公式改寫為：注意到，將上式進(jìn)一步改寫為：這是我們最終使用的計算格式。第十八頁，共五十八頁，編輯于2023年，星期五例子：取步長為0.1計算，結(jié)果如圖。第十九頁，共五十八頁，編輯于2023年，星期五圖：第二十頁，共五十八頁，編輯于2023年，星期五

DOUBLEPRECISIONh,y(0:10),ak1,ak2OPEN(20,FILE='OUTPUT1.DAT',STATUS="UNKNOWN") h=1.0/10 y(0)=1.0 do10i=1,10

ak1=-h*y(i-1) ak2=-h*(y(i-1)+ak1) y(i)=y(i-1)+(ak1+ak2)/2.010continuedo20i=0,10 write(20,*)i*h,y(i),exp(-i*h)20continue END第二十一頁，共五十八頁，編輯于2023年，星期五同理，對于后退歐拉公式有預(yù)測－校正公式或改寫為：第二十二頁，共五十八頁，編輯于2023年，星期五用此法解前面的例子步長0.1步長0.01第二十三頁，共五十八頁，編輯于2023年，星期五1.4誤差估計定義：利用第n個節(jié)點的函數(shù)值，通過數(shù)值方法計算第n+1個節(jié)點的函數(shù)近似值，所引起的誤差，稱為n+1個節(jié)點上的局部截斷誤差。我們記為n+1個節(jié)點上函數(shù)的精確值，為數(shù)值解，則局部截斷誤差為：如局部截斷誤差為，稱為具有p階局部截斷誤差。第二十四頁，共五十八頁，編輯于2023年，星期五歐拉方法的誤差分析：省略余項得公式：第二十五頁，共五十八頁，編輯于2023年，星期五完全類似的可以得到后退歐拉公式的局部截斷誤差為：歐拉中點公式的局部截斷誤差為：歐拉梯形公式的局部截斷誤差為：第二十六頁，共五十八頁，編輯于2023年，星期五定義：由初值引起的第n個節(jié)點的近似值與精確值之間的誤差稱為第n個節(jié)點整體誤差。定理：設(shè)下面求解微分方程的數(shù)值計算方法局部截斷誤差為p+1階，且函數(shù)關(guān)于y

滿足利普希茨條件，同時初值是準(zhǔn)確的，則整體截斷誤差為p階。歐拉公式、后退歐拉公式的整體截斷誤差為1階。歐拉中點公式、歐拉梯形公式的整體截斷誤差為2階。第二十七頁，共五十八頁，編輯于2023年，星期五微分方程數(shù)值解法的進(jìn)一步改進(jìn)。再回到恒等式如果取作為節(jié)點，將被積函數(shù)用插值多項式來近似，用插值多項式帶到積分中去求出積分，則可以得到所謂的亞當(dāng)斯(Adams)顯式公式局部截斷誤差：第二十八頁，共五十八頁，編輯于2023年，星期五類似地，如果取作為節(jié)點，可得亞當(dāng)斯(Adams)隱式公式局部截斷誤差：第二十九頁，共五十八頁，編輯于2023年，星期五進(jìn)一步，如果我們將恒等式中的積分區(qū)間改為，并在此區(qū)間上用辛普森公式，可得第三十頁，共五十八頁，編輯于2023年，星期五1.5絕對穩(wěn)定性一個常微分方程數(shù)值解法稱為絕對穩(wěn)定，是指在某一步（如第一步）產(chǎn)生的誤差（如計算機(jī)的存儲誤差），在計算中會逐步減小。稱方程為試驗方程，設(shè)計算步長為h，設(shè)在計算開始時產(chǎn)生誤差（存儲誤差），此誤差在以后會逐步減弱，我們稱該算法相對于是絕對穩(wěn)定的，的全體稱為該算法的絕對穩(wěn)定域。第三十一頁，共五十八頁，編輯于2023年，星期五歐拉法的絕對穩(wěn)定區(qū)域后退歐拉法的絕對穩(wěn)定區(qū)域第三十二頁，共五十八頁，編輯于2023年，星期五1.6局部截斷誤差的實用估計（1）用兩種階數(shù)相同的算法求解，計算出n+1步的近似值，從而得到局部截斷誤差估計。（2）用同樣的公式，用不同步長計算出n+1步的近似值，從而得到局部截斷誤差估計。1.7隱式法隱式法具有較好的絕對穩(wěn)定性！只不過在使用隱式法的時候，需要進(jìn)行迭代，或者使用預(yù)測－校正計算格式。第三十三頁，共五十八頁，編輯于2023年，星期五第二節(jié)泰勒級數(shù)法與龍格－庫塔法對于方程：取計算步長為h，則，將函數(shù)進(jìn)行泰勒展開如函數(shù)y(x)有p+1階導(dǎo)數(shù)，容易得到p階泰勒級數(shù)展開法：第三十四頁，共五十八頁，編輯于2023年，星期五公式中的導(dǎo)數(shù)用下面公式計算：例子：第三十五頁，共五十八頁，編輯于2023年，星期五步長0.1步長0.01第三十六頁，共五十八頁，編輯于2023年，星期五龍格－庫塔法：對于常微分方程的數(shù)值解法，一個關(guān)鍵在于選擇精度高的算法計算下面公式中的積分。要高精度的計算積分，常用的方法是適當(dāng)增加計算節(jié)點，不妨設(shè)用m

個節(jié)點計算上面積分，節(jié)點為則積分為第三十七頁，共五十八頁，編輯于2023年，星期五將積分改寫為：則得公式：這樣的公式稱為顯式龍格－庫塔公式。第三十八頁，共五十八頁，編輯于2023年，星期五確定二階R－K法:系數(shù)滿足：此為四個未知數(shù)的三個方程，任意取，得第三十九頁，共五十八頁，編輯于2023年，星期五特別，取，得到通常也稱為變形歐拉法，也常寫為它具有二階精度，也稱為二階二級R－K方法。第四十頁，共五十八頁，編輯于2023年，星期五例子第四十一頁，共五十八頁，編輯于2023年，星期五三階三級庫塔法局部截斷誤差為4階，整體截斷誤差為3階第四十二頁，共五十八頁，編輯于2023年，星期五最常用的四級四階公式，稱為龍格－庫塔公式：局部截斷誤差為5階，整體截斷誤差為4階。第四十三頁，共五十八頁，編輯于2023年，星期五隱式龍格－庫塔法：常用的有二階R－K法：第四十四頁，共五十八頁，編輯于2023年，星期五例子：用隱式二階R－K法：用顯式二階R－K法：第四十五頁，共五十八頁，編輯于2023年，星期五

DOUBLEPRECISIONh,y(0:100),yy(0:100) doubleprecisionakOPEN(20,FILE='OUTPUT5.DAT',STATUS="UNKNOWN") h=1.0/100 y(0)=1.0 yy(0)=1.0 do10i=1,100

y(i)=y(i-1)-(h/(1-h/2))*y(i-1) yy(i)=yy(i-1)*(1-h+h*h/2)10continuedo20i=0,100 write(20,100)i*h,y(i),yy(i),exp(-i*h)100format(1x,f4.2,f8.5,f8.5,f8.5)20continue END第四十六頁，共五十八頁，編輯于2023年，星期五h=0.1第四十七頁，共五十八頁，編輯于2023年，星期五h=0.01第四十八頁，共五十八頁，編輯于2023年，星期五半隱式龍格－庫塔法：第四十九頁，共五十八頁，編輯于2023年，星期五最常用的二級三階半隱式龍格－庫塔公式：第五十頁，共五十八頁，編輯于2023年，星期五微分方程組一階方程組第五十一頁，共五十八頁，編輯于2023年，星期五龍格－庫塔公式第五十二頁，共五十八頁，編輯于2023年，星期五寫成分量形式：第五十三頁，共五十八頁，編輯于2023年，星期五例子：第五十四頁，共五十八頁，編輯于2023年，星期五

DOUBLEPRECISIONh,y1(0:100),y2(0:100) doubleprecisionak1,ak2,ak3,ak4,bk1,bk2,bk3,bk4OPEN(20,FILE='OUTPUT6.DAT',STATUS="UNKNOWN") h=1.0/10 y1(0)=0.0 y2(0)=1.0 do10i=1,10

ak1=h*(-y1(i-1)+y2(i-1))

bk1=-h*y2(i-1)

ak2=h*(-(y1(i-1)+ak1/2.0)+y2(i-1)+bk1/2.0)

bk2=-h*(y2(i-1)+bk1/2.0)

ak3=h*(-(y1(i-1)+ak2/2.0)+y2(i-1)+bk2/2.0)

bk3=-h*(y2(i-1)+bk2/2.0)

ak4=h*(-(y1(i-1)+ak3)+y2(i-1)+bk3)

bk4=-h*(y2(i-1)+bk3)

y1(i)=y1(i-1)+(ak1+2*ak2+2*ak3+ak4)/6.0

y2(i)=y2(i-1)+(bk1+2*bk2+2*bk3+bk4)/6.010continuedo20i=0,10 write(20,100)i*h,y1(i),i*h*exp(-i*h),y2(i),exp(-i*h)100format(1x,f4.2,f12.5,f12.5,f12.5,f12.5)20continue END第五十五頁，共五十八頁，編輯于2023年，