第七章 微分方程

第七章微分方程第1頁，共37頁，2023年，2月20日，星期三§1引言常微分方程復(fù)習(xí)：一般的，凡表示未知函數(shù)、未知函數(shù)的導(dǎo)數(shù)與自變量之間關(guān)系的方程叫做微分方程。微分方程的解就是一個未知函數(shù)，將該未知函數(shù)代入微分方程中，能使方程成為恒等。

一階常微分方程的初值問題Initial-ValueProblem該微分方程的特解y(x)為一條曲線，稱為微分方程的積分曲線。第2頁，共37頁，2023年，2月20日，星期三

考慮一階常微分方程的初值問題/*Initial-ValueProblem*/:只要f(x,y)在[a,b]R1上連續(xù)，且關(guān)于y滿足Lipschitz

條件，即存在與x,y無關(guān)的常數(shù)L使對任意定義在[a,b]上的(x,y1)和(x,y2)都成立，則上述IVP存在唯一解y=y(x)。要計算出解函數(shù)y(x)在一系列節(jié)點a=x0<x1<…<xn=b

處的近似值節(jié)點間距為步長，通常采用等距節(jié)點，即取hi=h

(常數(shù))。第3頁，共37頁，2023年，2月20日，星期三數(shù)值解法的基本特點（思路）：求解過程順著節(jié)點排列次序一步步向前推進，即按遞推公式由已知的y0,y1,…yi，求出yi+1。yox0y(x)p0x1p1p2x2xxnPn+1Pnxn+1第4頁，共37頁，2023年，2月20日，星期三§2歐拉方法

/*Euler’sMethod*/將兩端在上積分已知依次求出y1,y2,…yn第5頁，共37頁，2023年，2月20日，星期三y(x)xyox0p0x1p1Pi+1xi+1xnPn+1Pnxn+1歐拉公式:yi+1

=yi+hf(xi,yi)歐拉公式的幾何意義近似代替y(x)，用該直線與x=xi+1的交點Pi＋1(xi＋1,yi＋1)的縱坐標(biāo)yi+1作為y(xi＋1)的近似值。xiPi(xi,yi)在區(qū)間用過點Pi(xi,yi)，以f(xi,yi)為斜率的直線y

=yi+f(xi,yi

)(x-xi)第6頁，共37頁，2023年，2月20日，星期三

歐拉公式：x0x1向前差商近似導(dǎo)數(shù)記為定義在假設(shè)yi=y(xi)，即第i步計算是精確的前提下，考慮的截斷誤差Ri=y(xi+1)

yi+1稱為局部截斷誤差/*localtruncationerror*/。定義若某算法的局部截斷誤差為O(hp+1)，則稱該算法有p階精度。

歐拉法的局部截斷誤差：歐拉法具有1階精度。Ri的主項/*leadingterm*/第7頁，共37頁，2023年，2月20日，星期三

歐拉公式的改進：隱式歐拉法向后差商近似導(dǎo)數(shù)x0x1))(,()(1101xyxfhyxy+)1,...,0(),(111-=+=+++niyxfhyyiiii由于未知數(shù)yi+1

同時出現(xiàn)在等式的兩邊，不能直接得到，故稱為隱式/*implicit*/

歐拉公式，而前者稱為顯式/*explicit*/歐拉公式。一般先用顯式計算一個初值，再迭代求解。隱式歐拉法的局部截斷誤差：即隱式歐拉公式具有1階精度。第8頁，共37頁，2023年，2月20日，星期三梯形公式—顯、隱式兩種算法的平均注：局部截斷誤差，即梯形公式具有2階精度，比歐拉方法有了進步。但注意到該公式是隱式公式，計算時不得不用到迭代法，其迭代收斂性與歐拉公式相似。中點歐拉公式/*midpointformula*/中心差商近似導(dǎo)數(shù)x0x2x1假設(shè)，則可以導(dǎo)出即中點公式具有2階精度。需要2個初值y0和y1來啟動遞推過程，這樣的算法稱為雙步法/*double-stepmethod*/，而前面的三種算法都是單步法/*single-stepmethod*/。第9頁，共37頁，2023年，2月20日，星期三方法顯式歐拉隱式歐拉梯形公式中點公式簡單精度低穩(wěn)定精度低,計算量大精度提高計算量大精度提高,顯式多一個初值,可能影響精度第10頁，共37頁，2023年，2月20日，星期三例題：用Euler公式和改進的Euler公式分別求下列初值問題的數(shù)值解(取步長h=0.1計算到y(tǒng)3)：y′=-2xy2

y(0)=1解：由歐拉公式

yn+1=yn+hf(xn,yn)=yn-2hxnyn2計算如下y1=y0-2hx0y02=1-2·0·0·12=1y2=y1-2hx1y12=1-2·0.1·0.1·12=0.98y3=y2-2hx2y22=1-2·0.1·0.2·0.982=0.9416第11頁，共37頁，2023年，2月20日，星期三改進歐拉法的預(yù)-校公式計算如下y1=0.99;y2=0.9614;y3=0.9173精確解y(0.1)=0.99,y(0.2)=0.9614;y(0.3)=0.9173可見改進歐拉公式比歐拉公式精度高第12頁，共37頁，2023年，2月20日，星期三改進歐拉法/*modifiedEuler’smethod*/Step1:先用顯式歐拉公式作預(yù)測，算出),(1iiiiyxfhyy+=+Step2:再將代入隱式梯形公式的右邊作校正，得到1+iy)],(),([2111+++++=iiiiiiyxfyxfhyy注：此法亦稱為預(yù)測-校正法/*predictor-correctormethod*/?？梢宰C明該算法具有2階精度，同時可以看到它是個單步遞推格式，比隱式公式的迭代求解過程簡單。第13頁，共37頁，2023年，2月20日，星期三為便于編制程序上機計算，改進歐拉公式可以改寫為如下形式：第14頁，共37頁，2023年，2月20日，星期三§3

Runge-Kutta方法[一]基本思想Runge-Kutta方法是一種高精度的單步法,簡稱R-K法．得到高精度方法的一個直接想法是利用Taylor展開．假設(shè)式y(tǒng)'=f(x,y)(a≤x≤b)

中的f(x,y)

充分光滑,將y(xn+1)在xn點作Taylor展開:y(xn+1)=y(xn)+hy'(xn)+(h2/2!)y''(xn)+......+(hp/p!)y(p)(xn)+..... 其中 y'(x)=f(x,y(x)) y''(x)=[f(x,y(x))]'x=fx+f·fy ............................... y(p)(x)=[f(x,y(x))](p)x第15頁，共37頁，2023年，2月20日，星期三對照標(biāo)準(zhǔn)形式y(tǒng)n+1=yn+hф(xn,yn;h).若取ф(x,y;h)=y'(x)+(h/2!)y''(x)+......+(hp-1/p!)y(p)(x)并以yn代替y(xn),則得到一個p階近似公式

yn+1=yn+hф(xn,yn;h)(n=0,1,2,......)(**)顯然p=1時,式(**)就是 yn+1=yn+hf(xn,yn)它即為我們熟悉的Euler方法.當(dāng)p≥2時,要利用公式(**)就需要計算f(x,y)的高階微商.這個計算量是很大的.因此，利用式(**)構(gòu)造高階公式是不實用的.第16頁，共37頁，2023年，2月20日，星期三R-K方法不是直接使用Taylor級數(shù),而是利用它的思想,即計算f(x,y)在不同結(jié)點的函數(shù)值,然后作這些函數(shù)值的線性組合,構(gòu)造近似公式,式中有一些可供選擇的參數(shù).將近似公式與Taylor展開式相比較,使前面的若干項密合,從而使近似公式達(dá)到一定的精度.下面以二級二階R-K方法為例說明這一方法的基本思想.第17頁，共37頁，2023年，2月20日，星期三[二]二級二階R-K方法在[xn,xn+1]上,取f(x,y)在兩個點的函數(shù)值作線性組合,記得到二級R-K方法:yn+1=yn+h(c1K1+c2K2)K1=f(xn,yn) （*） K2=f(xn+a2h,yn+b21hK1)其中c1,c2,a2,b21為待定參數(shù).對照式(**)有ф(x,y,h)=c1f(x,y)+c2f(x+a2h,y+b21hf(x,y))若要求式(**)達(dá)到二階精度，則只要局部截斷誤差Tn+1=O(h3).將y(xn+1)及ф(xn,y(xn),h)在xn作Taylor展開, f(xn,y(xn))=y'(xn) f(xn+a2h,y(xn)+b21hK1)=f(xn,y(xn)+a2hfx+b21hf·fy+O(h2))其中 fx=fx(xn,y(xn)),fy=fy(xn,y(xn)),f=f(xn,y(xn))由此得第18頁，共37頁，2023年，2月20日，星期三

ф(xn,y(xn);h)=(c1+c2)y'(xn)+c2(a2hfx+b21hfyf)+O(h2) 因為y(xn+1)在xn處的Taylor展開為 y(xn+1)=y(xn)+hy'(xn)+(h2/2!)y''(xn)+O(h3)由顯式單步法在xn+1的局部截斷誤差定義Tn+1=y(xn+1)-y(xn)-hф(xn,y(xn),h)=h(1-c1-c2)y'(xn)+h2[(1/2-a2c2)fx+(1/2-c2b21)fyf]+O(h3)顯然，若要求Tn+1=O(h3),則應(yīng)有c1+c2=1c2a2=1/2,c2b21=1/2上方程組含有3個方程,4個未知數(shù),其解是不唯一的.若取c2=α為自由參數(shù),則得其的一族解c1=1-α,c2=α,a2=b21=1/（2α） （***） 滿足條件（***)的(*)式稱為二級二階R-K方法.特別當(dāng)α=1/2時,公式即是前面介紹的改進的Euler方法.第19頁，共37頁，2023年，2月20日，星期三

當(dāng)α

=1時,c1=0,c2=1，得yn+1=yn+hK2n=0,1,….N-1K1=f(xn,yn)K2=f(xn+h/2,yn+hK1/2)這就是變形的歐拉方法或中點方法。當(dāng)α

=3/4時，c1=1/4c2=3/4，得到Heun方法。yn+1=yn+h(K1+3K2)/4,n=0,1,….N-1K1=f(xn,yn)K2=f(xn+2h/3,yn+2hK1/3)

第20頁，共37頁，2023年，2月20日，星期三二級R-K方法是顯示單步式,每前進一步需要計算兩個函數(shù)值.由上面的討論可知,適當(dāng)選擇四個參數(shù)c1,c2,a2,b21,可使每步計算兩次函數(shù)值的二階R-K方法達(dá)到二階精度.能否在計算函數(shù)值次數(shù)不變的情況下,通過選擇四個參數(shù),使得二階R-K方法的精度再提高呢?我們說,答案是否定的.無論四個參數(shù)怎樣選擇,都不能使公式(**)提高到三階.這說明每一步計算兩個函數(shù)值的二階R-K方法最高階為二階.若要獲得更高階得數(shù)值方法,就必須增加計算函數(shù)值的次數(shù).第21頁，共37頁，2023年，2月20日，星期三

m級顯式Runge-Kutta方法

仿照二級R-K方法,在[xn,xn+1]上,取f在m個點的函數(shù)值做線性組合,即得到m級R-K方法

m

yn+1=yn+h∑crKr

r=1 K1=f(xn,y

n) r-1 Kr=f(xn+har,yn+h∑brsKs)，r=2,3,....,m

s=1

第22頁，共37頁，2023年，2月20日，星期三使用不同的方法確定參數(shù)cr,ar,brs可使上式成為不同階的R-K方法.在m級R-K方法中,最著名的是經(jīng)典R-K方法:

yn+1=yn+(K1+2K2+2K3+K4)h/6 K1=f(xn,yn) K2=f(xn+h/2,yn+hK1/2） K3=f(xn+h/2,yn+hK2/2) K4=f(xn+h,yn+hK3)

由于它每前進一步需要計算四個點的函數(shù)值,因此稱為四級公式.按定義可直接驗證它的局部截斷誤差為O(h5),故它是四階方法.第23頁，共37頁，2023年，2月20日，星期三

前面已經(jīng)看到,二階、四階R-K方法可分別達(dá)到最高階數(shù)二階、四階，但是N階R-K方法的最高階卻不一定是N階。R-K方法的級數(shù)表示公式中計算函數(shù)值f的次數(shù)。Butcher給出了R-K方法計算函數(shù)值f的次數(shù)與階數(shù)之間的關(guān)系表，如下：

計算f的次數(shù) 1 23 4567 方法的最高階數(shù)12 3 4 4 5 6由表可見，四級以下R-K的方法其最高階數(shù)與計算f的次數(shù)一致，對m階R-K公式，當(dāng)m>4，雖然計算f的次數(shù)增加，但是方法階數(shù)不一定增加。因此四級四階R-K公式是應(yīng)用最為廣泛的公式。 第24頁，共37頁，2023年，2月20日，星期三

絕對穩(wěn)定性與絕對穩(wěn)定域求解初值問題的數(shù)值方法，當(dāng)給定不同步長計算時結(jié)果的舍入誤差影響差別很大，如果舍入誤差不增長算法就是數(shù)值穩(wěn)定的，若舍入誤差增長很快算法就不穩(wěn)定。定義：用一個數(shù)值方法求解微分方程y′=y是復(fù)數(shù)(1.5)對給定的步長h,在計算yn時引起的誤差n，若這個誤差在計算后面的yn+k中所引起的誤差n+k滿足:|n+k|≤|n|(k=1,2,…)就說這個數(shù)值方法對步長h和復(fù)數(shù)是絕對穩(wěn)定的,使得數(shù)值方法是絕對穩(wěn)定的H=h在復(fù)平面上的允許范圍稱為數(shù)值方法的絕對穩(wěn)定域.實驗方程第25頁，共37頁，2023年，2月20日，星期三特拉法加爾（Trafalgar）海戰(zhàn)和納爾森（Nelson）秘訣19世紀(jì)中葉，法國拿破倫統(tǒng)帥大軍要與英國爭奪海上霸主地位，而實施這一戰(zhàn)略的最主要的關(guān)鍵是消滅英國的艦隊。英國海軍統(tǒng)帥、海軍中將納爾森親自制定了周密的戰(zhàn)術(shù)方案。第26頁，共37頁，2023年，2月20日，星期三1805年10月21日，這場海上大戰(zhàn)爆發(fā)了。英國是納爾森親自統(tǒng)帥的地中海艦隊，由27艘戰(zhàn)艦組成；另外一方是由費倫紐夫（Villenuve）率領(lǐng)的法國——西班牙聯(lián)合艦隊，共有33艘戰(zhàn)艦。Trafalgar大海戰(zhàn)的概況是：費倫紐夫（Villenuve）率領(lǐng)的法國——西班牙聯(lián)合艦隊采用常規(guī)的一字橫列，以利炮火充分展開，而納爾森的戰(zhàn)術(shù)使費倫紐夫大出意外。第27頁，共37頁，2023年，2月20日，星期三英國的艦隊分成兩個縱列：前衛(wèi)上風(fēng)縱列由12艘戰(zhàn)艦組成，由納爾森親自指揮，攔腰將法國——西班牙聯(lián)合艦隊切為兩段；后衛(wèi)下風(fēng)縱列由英國海軍中將科林伍德(Collingwood)指揮，由15艘戰(zhàn)艦組成。在一場海戰(zhàn)后，法國——西班牙聯(lián)合艦隊以慘敗告終：聯(lián)合艦隊司令費倫紐夫連同12艘戰(zhàn)艦被俘，8艘沉沒，僅13艘逃走，人員傷亡7000人。而英國戰(zhàn)艦沒有沉沒，人員傷亡1663人，但是，作為統(tǒng)帥的納爾森陣亡。第28頁，共37頁，2023年，2月20日，星期三秘密備忘錄中的納爾森（Nelson）秘訣：預(yù)期參加戰(zhàn)斗的英國艦隊：40艘。法國—西班牙聯(lián)合艦隊：46艘。預(yù)計聯(lián)合艦隊?wèi)?zhàn)斗隊形一字橫列。英國艦隊的戰(zhàn)斗隊形與任務(wù)：分成兩個主縱列及一個小縱列。第29頁，共37頁，2023年，2月20日，星期三主縱列1：16艘，由納爾森親自指揮，攔腰將法國——西班牙聯(lián)合艦隊切為兩段，并攻擊聯(lián)合艦隊的中間部分。主縱列2：16艘，由英國海軍中將科林伍德指揮，從聯(lián)合艦隊后半部再切斷，分割并攻擊后部12艘。小縱列：8艘，在中心部分附近攻擊其先頭部分的3-4艘。第