1、School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 引言引言 包含自變量、未知函數(shù)及未知函數(shù)的導(dǎo)數(shù)或微分的方包含自變量、未知函數(shù)及未知函數(shù)的導(dǎo)數(shù)或微分的方程稱為微分方程。在微分方程中程稱為微分方程。在微分方程中, 自變量的個(gè)數(shù)只有一個(gè)自變量的個(gè)數(shù)只有一個(gè), 稱為常微分方程。自變量的個(gè)數(shù)為兩個(gè)或兩個(gè)以上的微分稱為常微分方程。自變量的個(gè)數(shù)為兩個(gè)或兩個(gè)以上的微分方程叫偏微分方程。微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)方程叫偏微分方程。微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱為微分方程的階數(shù)。如

2、果未知函數(shù)數(shù)的階數(shù)稱為微分方程的階數(shù)。如果未知函數(shù)y及其各階導(dǎo)及其各階導(dǎo)數(shù)數(shù)都是一次的都是一次的,則稱它是線性的則稱它是線性的,否則稱為非線性的。否則稱為非線性的。 )(,nyyy 第三章第三章 常微分方程的差分方法常微分方程的差分方法School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 在高等數(shù)學(xué)中，對(duì)于常微分方程的求解，給出了一些典在高等數(shù)學(xué)中，對(duì)于常微分方程的求解，給出了一些典型方程求解析解的基本方法，如可分離變量法、常系數(shù)齊次型方程求解析解的基本方法，如可分離變量法、常系數(shù)齊次

3、線性方程的解法、常系數(shù)非齊次線性方程的解法等。但能求線性方程的解法、常系數(shù)非齊次線性方程的解法等。但能求解的常微分方程仍然是有限的，大多數(shù)的常微分方程是不可解的常微分方程仍然是有限的，大多數(shù)的常微分方程是不可能給出解析解。能給出解析解。 譬如譬如 22yxy 這個(gè)一階微分方程就不能用初等函數(shù)及其積分來(lái)表達(dá)這個(gè)一階微分方程就不能用初等函數(shù)及其積分來(lái)表達(dá)它的解。它的解。 School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 從實(shí)際問題當(dāng)中歸納出來(lái)的微分方程，通常主要依靠數(shù)值解從實(shí)際問題當(dāng)中

4、歸納出來(lái)的微分方程，通常主要依靠數(shù)值解法來(lái)解決。本章主要討論一階常微分方程初值問題法來(lái)解決。本章主要討論一階常微分方程初值問題 00)(),(yxyyxfy ( 3.1 ) 在區(qū)間在區(qū)間a x b上的數(shù)值解法上的數(shù)值解法。 可以證明可以證明,如果函數(shù)在帶形區(qū)域如果函數(shù)在帶形區(qū)域 R=axb,-y內(nèi)連續(xù)，內(nèi)連續(xù)，且關(guān)于且關(guān)于y滿足李普希茲滿足李普希茲(Lipschitz)條件，即存在常數(shù)條件，即存在常數(shù)L(它與它與x,y無(wú)關(guān)無(wú)關(guān))使使 2121),(),(yyLyxfyxf對(duì)對(duì)R內(nèi)任意兩個(gè)內(nèi)任意兩個(gè) 都成立都成立,則方程則方程( 3.1 )的解的解在在 a, b 上存在且唯一。上存在且唯一。 2

5、1, yy)(xyy School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 7.2 數(shù)值方法的基本思想數(shù)值方法的基本思想 對(duì)常微分方程初值問題對(duì)常微分方程初值問題(7.1)式的數(shù)值解法，就是要算出精確式的數(shù)值解法，就是要算出精確解解y(x)在區(qū)間在區(qū)間 a,b 上的一系列離散節(jié)點(diǎn)上的一系列離散節(jié)點(diǎn) 處的函數(shù)值處的函數(shù)值 的近似值的近似值 相鄰兩個(gè)節(jié)點(diǎn)的間距相鄰兩個(gè)節(jié)點(diǎn)的間距 稱為步長(zhǎng)，步長(zhǎng)可以相等，稱為步長(zhǎng)，步長(zhǎng)可以相等，也可以不等。本章總是假定也可以不等。本章總是假定h為定數(shù)，稱為定

6、步長(zhǎng)，這時(shí)節(jié)點(diǎn)為定數(shù)，稱為定步長(zhǎng)，這時(shí)節(jié)點(diǎn)可表示為可表示為數(shù)值解法需要把連續(xù)性的問題加以離散化，從而求出離散節(jié)點(diǎn)數(shù)值解法需要把連續(xù)性的問題加以離散化，從而求出離散節(jié)點(diǎn)的數(shù)值解。的數(shù)值解。 bxxxxann110)(,),(),(10nxyxyxynyyy,10iixxh1niihxxi, 2 , 1 ,0School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 對(duì)常微分方程數(shù)值解法的基本出發(fā)點(diǎn)就是離散化。其數(shù)對(duì)常微分方程數(shù)值解法的基本出發(fā)點(diǎn)就是離散化。其數(shù)值解法有兩個(gè)基本特點(diǎn)，它們都采用值

7、解法有兩個(gè)基本特點(diǎn)，它們都采用“步進(jìn)式步進(jìn)式”，即求解過，即求解過程順著節(jié)點(diǎn)排列的次序一步一步地向前推進(jìn)，描述這類算法，程順著節(jié)點(diǎn)排列的次序一步一步地向前推進(jìn)，描述這類算法，要求給出用已知信息要求給出用已知信息 計(jì)算計(jì)算 的遞推的遞推公式。建立這類遞推公式的基本方法是在這些節(jié)點(diǎn)上用數(shù)值公式。建立這類遞推公式的基本方法是在這些節(jié)點(diǎn)上用數(shù)值積分、數(shù)值微分、泰勒展開等離散化方法，對(duì)初值問題積分、數(shù)值微分、泰勒展開等離散化方法，對(duì)初值問題中的導(dǎo)數(shù)中的導(dǎo)數(shù) 進(jìn)行不同的離散化處理進(jìn)行不同的離散化處理。 021,yyyyiii1iyy00)(),(yxyyxfySchool of Automation En

8、gineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 對(duì)于初值問題對(duì)于初值問題的數(shù)值解法，首先要解決的問題就是如何對(duì)微分方程進(jìn)行離的數(shù)值解法，首先要解決的問題就是如何對(duì)微分方程進(jìn)行離散化，建立求數(shù)值解的遞推公式。遞推公式通常有兩類，一散化，建立求數(shù)值解的遞推公式。遞推公式通常有兩類，一類是計(jì)算類是計(jì)算yi+1時(shí)只用到時(shí)只用到xi+1, xi 和和yi，即前一步的值，因此有了初，即前一步的值，因此有了初值以后就可以逐步往下計(jì)算，此類方法稱為單步法；其代表值以后就可以逐步往下計(jì)算，此類方法稱為單步法；其代表是龍格是龍格庫(kù)塔法。另一類是

9、計(jì)算庫(kù)塔法。另一類是計(jì)算yi+1時(shí)，除用到時(shí)，除用到xi+1,xi和和yi以外，以外，還要用到還要用到 ，即前面，即前面k步的值，此類方法步的值，此類方法稱為多步法；其代表是亞當(dāng)斯法。稱為多步法；其代表是亞當(dāng)斯法。 00)(),(yxyyxfy), 2 , 1( ,kpyxpipiSchool of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.1 歐拉（歐拉（Euler）法）法3.1.1 Euler公式公式 歐拉（歐拉（Euler）方法是解初值問題的最簡(jiǎn)單的數(shù)值方法。）方法是解初值問題的最簡(jiǎn)

10、單的數(shù)值方法。初值問題初值問題的解的解y=y(x)代表通過點(diǎn)代表通過點(diǎn) 的一條稱之為微分方程的積的一條稱之為微分方程的積分曲線。積分曲線上每一點(diǎn)分曲線。積分曲線上每一點(diǎn) 的切線的斜率的切線的斜率 等于等于函數(shù)函數(shù) 在這點(diǎn)的值。在這點(diǎn)的值。 00)(),(yxyyxfy),(00yx),(yx)(xy),(yxfSchool of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 Pi+1 Pn y=y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 Euler法的

11、求解過程是法的求解過程是:從初始點(diǎn)從初始點(diǎn)P0(即點(diǎn)即點(diǎn)(x0,y0)出發(fā)出發(fā),作積分曲線作積分曲線y=y(x)在在P0點(diǎn)上切線點(diǎn)上切線 (其斜率其斜率為為 ),與與x=x1直線直線10PP),()(000yxfxy相交于相交于P1點(diǎn)點(diǎn)(即點(diǎn)即點(diǎn)(x1,y1),得到得到y(tǒng)1作為作為y(x1)的近似值的近似值,如上圖所示。如上圖所示。過點(diǎn)過點(diǎn)(x0,y0),以以f(x0,y0)為斜率的切線方程為為斜率的切線方程為 當(dāng)當(dāng) 時(shí)時(shí),得得 )(,(0000 xxyxfyy1xx )(,(010001xxyxfyy這樣就獲得了這樣就獲得了P1點(diǎn)的坐標(biāo)。點(diǎn)的坐標(biāo)。 School of Automation

12、Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 x0 x1School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 Pi+1 Pn y=y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 同樣同樣, 過點(diǎn)過點(diǎn)P1(x1,y1),作積分曲線作積分曲線y=y(x)的切線的切線交直線交直線x=x2于于P2點(diǎn)點(diǎn),切線切線 的斜率的斜率 =直線方程為直線方程為21PP)(1xy),(11y

13、xf)(,(1111xxyxfyy)(,(121112xxyxfyy當(dāng)當(dāng) 時(shí)時(shí),得得 2xx School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 當(dāng)當(dāng) 時(shí)時(shí),得得 Pi+1 Pn y= y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 由此獲得了由此獲得了P2的坐標(biāo)。重復(fù)以上過程的坐標(biāo)。重復(fù)以上過程,就可獲得一系列的就可獲得一系列的點(diǎn)點(diǎn):P1,P1,Pn。對(duì)已求得點(diǎn)。對(duì)已求得點(diǎn)以以 = 為斜率作直線為斜率作直線 ),(nnnyxP)(nxy),

14、(nnyxf)(,(nnnnxxyxfyy1nxx)(,(11nxnnnnxxyxfyynnyxy)(取取School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 從圖形上看從圖形上看,就獲得了一條近似于曲線就獲得了一條近似于曲線y=y(x)的折線的折線 。 Pi+ 1 Pn y= y(x) P1 Pi Pn Pi+ 1 P0 x0 x1 xi xi+ 1 xn Pi P1 這樣這樣,從從x0 逐個(gè)算出逐個(gè)算出對(duì)應(yīng)的數(shù)值解對(duì)應(yīng)的數(shù)值解 nxxx,21nyyy,21nPPPP321Schoo

15、l of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 通常取通常取 (常數(shù)常數(shù)),則則Euler法的計(jì)算格式法的計(jì)算格式 hhxxiii1)(),(001xyyyxhfyyiiii i=0,1,n ( 3.2 ) 還可用數(shù)值微分、數(shù)值積分法和泰勒展開法推導(dǎo)還可用數(shù)值微分、數(shù)值積分法和泰勒展開法推導(dǎo)Euler格式。格式。以數(shù)值積分為例進(jìn)行推導(dǎo)。以數(shù)值積分為例進(jìn)行推導(dǎo)。將方程將方程 的兩端在區(qū)間的兩端在區(qū)間 上積分得，上積分得， ),(yxfy 1,iixx11),(iiiixxxxdxyxfdxy

16、11)(,)(),()()(1iiiixxixxiidxxyxfxydxyxfxyxy選擇不同的計(jì)算方法計(jì)算上式的積分項(xiàng)選擇不同的計(jì)算方法計(jì)算上式的積分項(xiàng) ,就就會(huì)得到不同的計(jì)算公式。會(huì)得到不同的計(jì)算公式。 1)(,iixxdxxyxf(3.3)School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 用左矩形方法計(jì)算積分項(xiàng)用左矩形方法計(jì)算積分項(xiàng) )(,)()(,11iiiixxxyxfxxdxxyxfii代入代入(3.3)式式,并用并用yi近似代替式中近似代替式中y(xi)即可得到向前歐

17、拉即可得到向前歐拉（Euler）公式）公式 ),(1iiiiyxhfyy 由于數(shù)值積分的矩形方法精度很低，所以歐拉（由于數(shù)值積分的矩形方法精度很低，所以歐拉（Euler）公式當(dāng)然很粗糙。公式當(dāng)然很粗糙。 School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 例例3.1 用歐拉法解初值問題用歐拉法解初值問題 1)0()6 . 00(2yxxyyy取步長(zhǎng)取步長(zhǎng)h=0.2 ,計(jì)算過程保留計(jì)算過程保留4位小數(shù)位小數(shù) 解解: h=0.2, 歐拉迭代格式歐拉迭代格式 2),(xyyyxf21),(

18、iiiiiiiiyhxhyyyxhfyy)2 , 1 , 0()4(2 . 0iyxyiii當(dāng)當(dāng) k=0, x1=0.2時(shí)，已知時(shí)，已知x0=0，y0=1，有，有 y(0.2) y1=0.21(401)0.8當(dāng)當(dāng) k=1, x2=0.4時(shí)，已知時(shí)，已知x1 =0.2， y1 =0.8，有，有 y(0.4) y2 =0.20.8(40.20.8)0.6144當(dāng)當(dāng) k=2, x3 =0.6時(shí)，已知時(shí)，已知x2 =0.4， y2 =0.6144，有，有 y(0.6) y3=0.20.6144(4-0.40.6144)=0.4613 School of Automation Engineering自自

19、 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.1.2 梯形公式梯形公式為了提高精度為了提高精度,對(duì)方程對(duì)方程 的兩端在區(qū)間上的兩端在區(qū)間上積分得，積分得，改用梯形方法計(jì)算其積分項(xiàng)，即改用梯形方法計(jì)算其積分項(xiàng)，即 ),(yxfy 1,iixx1)(,)()(1iixxiidxxyxfxyxy)(,()(,(2)(,1111iiiiiixxxyxfxyxfxxdxxyxfii( 3.4 ) 代入代入(3.4)式式,并用近似代替式中即可得到梯形公式并用近似代替式中即可得到梯形公式 ),(),(2111iiiiiiyxfyxfhyy( 3.5 )

20、由于數(shù)值積分的梯形公式比矩形公式的精度高，因此梯由于數(shù)值積分的梯形公式比矩形公式的精度高，因此梯形公式（形公式（3.5）比歐拉公式）比歐拉公式( 3.2 )的精度高一個(gè)數(shù)值方法。的精度高一個(gè)數(shù)值方法。 School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 ),(),(2111iiiiiiyxfyxfhyy( 3.5 ) (3.5)式的右端含有未知的式的右端含有未知的yi+1,它是一個(gè)關(guān)于它是一個(gè)關(guān)于yi+1的函數(shù)的函數(shù)方程方程,這類數(shù)值方法稱為隱式方法。相反地這類數(shù)值方法稱為隱式方法。

21、相反地,歐拉法是關(guān)于歐拉法是關(guān)于yi+1的一個(gè)直接的計(jì)算公式，的一個(gè)直接的計(jì)算公式， 這類數(shù)值方法稱為顯式方法。這類數(shù)值方法稱為顯式方法。 School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.1.3 兩步歐拉公式兩步歐拉公式 對(duì)方程對(duì)方程 的兩端在區(qū)間上的兩端在區(qū)間上 積分得積分得 ),(yxfy 11,iixx11)(,)()(11iixxiidxxyxfxyxy ( 3.6 ) 改用中矩形公式計(jì)算其積分項(xiàng)，即改用中矩形公式計(jì)算其積分項(xiàng)，即 )(,)(,1111iiiixxxy

22、xfxxdxxyxfii代入上式代入上式,并用并用yi近似代替式中近似代替式中y(xi)即可得到兩步歐拉公式即可得到兩步歐拉公式 ),(211iiiiyxhfyy ( 3.7 ) School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 前面介紹過的數(shù)值方法前面介紹過的數(shù)值方法,無(wú)論是歐拉方法無(wú)論是歐拉方法,還是梯形方還是梯形方法，它們都是單步法法，它們都是單步法,其特點(diǎn)是在計(jì)算其特點(diǎn)是在計(jì)算yi+1時(shí)只用到前一步時(shí)只用到前一步的信息的信息yi;可是公式可是公式(3.7)中除了中除了yi外

23、外,還用到更前一步的信息還用到更前一步的信息yi-1,即調(diào)用了前兩步的信息即調(diào)用了前兩步的信息,故稱其為兩步歐拉公式。故稱其為兩步歐拉公式。School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.1.4 歐拉法的局部截?cái)嗾`差歐拉法的局部截?cái)嗾`差 衡量求解公式好壞的一個(gè)主要標(biāo)準(zhǔn)是求解公式的精度衡量求解公式好壞的一個(gè)主要標(biāo)準(zhǔn)是求解公式的精度, 因因此引入局部截?cái)嗾`差和階數(shù)的概念。此引入局部截?cái)嗾`差和階數(shù)的概念。 定義定義3.1 在在yi準(zhǔn)確的前提下準(zhǔn)確的前提下, 即即 時(shí)時(shí), 用數(shù)值方法

24、計(jì)算用數(shù)值方法計(jì)算yi+1的誤差的誤差 , 稱為該數(shù)值方法計(jì)算時(shí)稱為該數(shù)值方法計(jì)算時(shí)yi+1的局部的局部截?cái)嗾`差。截?cái)嗾`差。 對(duì)于歐拉公式，假定對(duì)于歐拉公式，假定 ，則有，則有)(iixyy 11)(iiiyxyR)(iixyy )()()(,()(1iiiiiixyhxyxyxfhxyy而將真解而將真解y(x)在在xi處按二階泰勒展開：處按二階泰勒展開： ),()(! 2)()()(121 iiiiixxyhxyhxyxy)(!2)(211yhyxyii 因此有因此有 School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常

25、微分方程的差分方法常微分方程的差分方法 定義定義3.2 數(shù)值方法的局部截?cái)嗾`差為數(shù)值方法的局部截?cái)嗾`差為 ,則稱這種數(shù)值方則稱這種數(shù)值方法的階數(shù)是法的階數(shù)是P。步長(zhǎng)。步長(zhǎng)(h N 結(jié)束。結(jié)束。 10,xx)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy11, yx0101,yyxxSchool of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 （2）改進(jìn)歐拉法的流程圖）改進(jìn)歐拉法的流程圖 開 始 輸 入x0, y0,h , N 1 n x0 + h x1 y0+ h

26、f( x0,y0 ) yp y0+ h f( x1,yp) yc ( yp+ yc) /2 y1 輸 出x1, y1 n + 1 n n = N ? x1 x0 y1 y0 結(jié) 束 n y School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 例例3.2 用改進(jìn)歐拉法解初值問題用改進(jìn)歐拉法解初值問題 1)0(2yyxyy區(qū)間為區(qū)間為 0,1 ,取步長(zhǎng)取步長(zhǎng)h=0.1。解解: 改進(jìn)歐拉法的具體形式改進(jìn)歐拉法的具體形式 )(21)2(1.0)2(1.011cpipipiciiiipyyyy

27、xyyyyxyyy本題的精確解為本題的精確解為 ,計(jì)算請(qǐng)對(duì)比計(jì)算請(qǐng)對(duì)比P98列表所示。列表所示。 xxy21)(School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 例例3.3 對(duì)初值問題對(duì)初值問題 1)0(0yyy證明用梯形公式求得的近似解為證明用梯形公式求得的近似解為 nnhhynhx 證明證明: 解初值問題的梯形公式為解初值問題的梯形公式為 ),(),(nnnnnnyxfyxfhyyyyxf),( 211nnnnyyhyy 整理成顯式整理成顯式 nnyhhy反復(fù)迭代反復(fù)迭代,得到

28、得到 yhhyhhyhhyhhynnnnn.10ynnhhy22 School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.2 龍格龍格-庫(kù)塔（庫(kù)塔（Runge-Kutta）法）法3.2.1 龍格龍格-庫(kù)塔庫(kù)塔(Runge-Kutta)法的基本思想法的基本思想 Euler公式可改寫成公式可改寫成 ),(111iiiiyxfKhKyy則則yi+1的表達(dá)式的表達(dá)式y(tǒng)(xi+1)與的與的Taylor展開式的前兩項(xiàng)完全相同展開式的前兩項(xiàng)完全相同,即局部截?cái)嗾`差為即局部截?cái)嗾`差為 。 改進(jìn)的改進(jìn)

29、的Euler公式又可改寫成公式又可改寫成 )(2hO),(),()(21121211hKyxfKyxfKKKhyyiiiiiiSchool of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 上述兩組公式在形式上有一個(gè)共同點(diǎn)上述兩組公式在形式上有一個(gè)共同點(diǎn):都是用都是用f(x,y)在某在某些點(diǎn)上值的線性組合得出些點(diǎn)上值的線性組合得出y(xi+1)的近似值的近似值yi+1,而且增加計(jì)算的而且增加計(jì)算的次數(shù)次數(shù)f(x,y)的次數(shù)的次數(shù),可提高截?cái)嗾`差的階。如歐拉公式可提高截?cái)嗾`差的階。如歐拉公式:每步

30、計(jì)每步計(jì)算一次算一次f(x,y)的值的值,為一階方法。改進(jìn)歐拉公式需計(jì)算兩次為一階方法。改進(jìn)歐拉公式需計(jì)算兩次f(x,y)的值，它是二階方法。它的局部截?cái)嗾`差為的值，它是二階方法。它的局部截?cái)嗾`差為 。)(3hOSchool of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 于是可考慮用函數(shù)于是可考慮用函數(shù)f(x,y)在若干點(diǎn)上的函數(shù)值的線性組合在若干點(diǎn)上的函數(shù)值的線性組合來(lái)構(gòu)造近似公式，構(gòu)造時(shí)要求近似公式在來(lái)構(gòu)造近似公式，構(gòu)造時(shí)要求近似公式在(xi,yi)處的處的Taylor展開展開式與解式與

31、解y(x)在在xi處的處的Taylor展開式的前面幾項(xiàng)重合，從而使近展開式的前面幾項(xiàng)重合，從而使近似公式達(dá)到所需要的階數(shù)。既避免求偏導(dǎo)似公式達(dá)到所需要的階數(shù)。既避免求偏導(dǎo),又提高了計(jì)算方法又提高了計(jì)算方法精度的階數(shù)?；蛘哒f精度的階數(shù)。或者說,在在 這一步內(nèi)多預(yù)報(bào)幾個(gè)點(diǎn)的斜率這一步內(nèi)多預(yù)報(bào)幾個(gè)點(diǎn)的斜率值，然后將其加權(quán)平均作為平均斜率，則可構(gòu)造出更高精度的值，然后將其加權(quán)平均作為平均斜率，則可構(gòu)造出更高精度的計(jì)算格式，這就是龍格計(jì)算格式，這就是龍格庫(kù)塔（庫(kù)塔（Runge-Kutta）法的基本思想。）法的基本思想。 1,iixxSchool of Automation Engineering自自 動(dòng)

32、動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.2.2 二階龍格二階龍格庫(kù)塔法庫(kù)塔法 在在 上取兩點(diǎn)上取兩點(diǎn)xi和和 ,以該兩點(diǎn)處的斜率值以該兩點(diǎn)處的斜率值k1和和k2的加權(quán)平均的加權(quán)平均(或稱為線性組合或稱為線性組合)來(lái)求取平均斜率來(lái)求取平均斜率k*的近似值的近似值K，即，即 1,iixxphxxipi2211kkK式中式中:k1為為xi點(diǎn)處的切線斜率值，點(diǎn)處的切線斜率值， k2為為 點(diǎn)處的切線斜率值點(diǎn)處的切線斜率值,比照改進(jìn)的歐拉法比照改進(jìn)的歐拉法,將將 視為視為 ，即可得，即可得 )(),(1iiixyyxfkpixphxi1ix),(12

33、phkyphxfkii對(duì)常微分方程初值問題對(duì)常微分方程初值問題(3.1)式的解式的解 y=y(x),根據(jù)微分中值定理，根據(jù)微分中值定理，存在點(diǎn)存在點(diǎn) ，使得，使得 ),(1iixx)()()(11iiiixxyxyxySchool of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 )(,()(yfyK式中式中 K可看作是可看作是y=y(x)在區(qū)間在區(qū)間 上的平均斜率。所以可上的平均斜率。所以可得計(jì)算公式為：得計(jì)算公式為： 1,iixxhKxyxyii)()(1)()(2211kkhxyi（3.1

34、4） 將將y(xi)在在x=xi處進(jìn)行二階處進(jìn)行二階Taylor展開：展開： )()(! 2)()()(321hOxyhxyhxyxyiiii （3.15） 也即也即 hKxyxyii)()(1（3.13）School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 將將 在在x=xi處進(jìn)行一階處進(jìn)行一階Taylor展開：展開： ),()(12phkyphxfphxykiii)(),(),(),(),(22hOyxfyxfyxfphyxfkiiyiiiixii)()()(2hOxyphxyii

35、 將以上結(jié)果代入（將以上結(jié)果代入（3.14）得：）得： )()()(22111kkhxyxyii)()()()()(221hOxyphxyxyhxyiiii )()()()()(32221hOxyphxyhxyiii （3.16） 對(duì)式對(duì)式(3.15)和和(3.16)進(jìn)行比較系數(shù)后可知進(jìn)行比較系數(shù)后可知,只要只要 211221p（3.17） 成立成立,格式格式(3.14)的局部截?cái)嗾`差就等于的局部截?cái)嗾`差就等于)(3hO有有2階階精度精度School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分

36、方法 式式(3.17)中具有三個(gè)未知量中具有三個(gè)未知量,但只有兩個(gè)方程但只有兩個(gè)方程,因而有無(wú)窮多解。因而有無(wú)窮多解。若取若取 ,則則p=1，這是無(wú)窮多解中的一個(gè)解，將以上所，這是無(wú)窮多解中的一個(gè)解，將以上所解的值代入式解的值代入式(3.14)并改寫可得并改寫可得 2121),(),()(21121211hkyxfkyxfkkkhyyiiiiii 不難發(fā)現(xiàn)，上面的格式就是改進(jìn)的歐拉格式。凡滿足不難發(fā)現(xiàn)，上面的格式就是改進(jìn)的歐拉格式。凡滿足條件式（條件式（3.17）有一簇形如上式的計(jì)算格式，這些格式統(tǒng)稱）有一簇形如上式的計(jì)算格式，這些格式統(tǒng)稱為二階龍格為二階龍格庫(kù)塔格式。因此改進(jìn)的歐拉格式是眾多

37、的二庫(kù)塔格式。因此改進(jìn)的歐拉格式是眾多的二階龍格階龍格庫(kù)塔法中的一種特殊格式。庫(kù)塔法中的一種特殊格式。 School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 若取若取 ,則則 ，此時(shí)二階龍格，此時(shí)二階龍格-庫(kù)塔庫(kù)塔法的計(jì)算公式為法的計(jì)算公式為 0121, 12p)2,(),(1212121khyxfkyxfkhkyyiiiiii1,2 , 1 , 0ni 此計(jì)算公式稱為變形的二階龍格此計(jì)算公式稱為變形的二階龍格庫(kù)塔法。式中庫(kù)塔法。式中 為區(qū)間為區(qū)間 的中點(diǎn)。的中點(diǎn)。 21ix1,iix

38、xSchool of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.2.3 三階龍格三階龍格-庫(kù)塔法庫(kù)塔法 為了進(jìn)一步提高精度，設(shè)除為了進(jìn)一步提高精度，設(shè)除 外再增加一點(diǎn)外再增加一點(diǎn) pix) 1(qpqhxxiqi并用三個(gè)點(diǎn)并用三個(gè)點(diǎn) , , 的斜率的斜率k1,k2,k3加權(quán)平均加權(quán)平均得出平均斜率得出平均斜率k*的近似值，這時(shí)計(jì)算格式具有形式的近似值，這時(shí)計(jì)算格式具有形式: ixpixqix),(),()1 (1213211phkyphxfkyxfkkkkhyyiiiiii（3.18） 為

39、了預(yù)報(bào)點(diǎn)為了預(yù)報(bào)點(diǎn) 的斜率值的斜率值k3,在區(qū)間在區(qū)間 內(nèi)有兩個(gè)斜率值內(nèi)有兩個(gè)斜率值k1和和k2可以用可以用,可將可將k1,k2加權(quán)平均得出加權(quán)平均得出 上的平均斜率上的平均斜率,從而從而得到得到 的預(yù)報(bào)值的預(yù)報(bào)值 qixqiixx,qiixx,)(qixyqiy21)1 (kkqhyyiqiSchool of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 于是可得于是可得 ),(3qiqiyxfk運(yùn)用運(yùn)用Taylor展開方法選擇參數(shù)展開方法選擇參數(shù) ,可以使格式可以使格式(3.18)的局部截?cái)嗾`

40、差為的局部截?cái)嗾`差為 ,即具有三階精度，這類格式統(tǒng)稱為三即具有三階精度，這類格式統(tǒng)稱為三階龍格階龍格庫(kù)塔方法。下列是其中的一種，稱為庫(kù)塔（庫(kù)塔方法。下列是其中的一種，稱為庫(kù)塔（Kutta）公式。公式。 ,qp)(4hO)4(6)2(,()2,(),(3211211312121kkkhyykkhyxfkkhyxfkyxfkiiiiiiii（3.19） School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.2.4 四階龍格四階龍格庫(kù)塔法庫(kù)塔法 如果需要再提高精度，用類似上述的處理方法，

41、只需在如果需要再提高精度，用類似上述的處理方法，只需在區(qū)間區(qū)間 上用四個(gè)點(diǎn)處的斜率加權(quán)平均作為平均斜率上用四個(gè)點(diǎn)處的斜率加權(quán)平均作為平均斜率k*的近的近似值，構(gòu)成一系列四階龍格似值，構(gòu)成一系列四階龍格庫(kù)塔公式。具有四階精度，即局庫(kù)塔公式。具有四階精度，即局部截?cái)嗾`差是部截?cái)嗾`差是 。 由于推導(dǎo)復(fù)雜，這里從略，只介紹最常用的一種四階經(jīng)典龍由于推導(dǎo)復(fù)雜，這里從略，只介紹最常用的一種四階經(jīng)典龍格格庫(kù)塔公式。庫(kù)塔公式。 qiixx,)(5hO)22(6),()2,()2,(),(43211314221312121kkkkhyyhkyxfkkhyxfkkhyxfkyxfkiiiiiiiiii（3.20）

42、 School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.2.5 四階龍格四階龍格庫(kù)塔法算法實(shí)現(xiàn)庫(kù)塔法算法實(shí)現(xiàn)(1) 計(jì)算步驟計(jì)算步驟 輸入輸入 ,h,N ; 使用龍格使用龍格庫(kù)塔公式（庫(kù)塔公式（3.20）計(jì)算出）計(jì)算出y1; 輸出輸出 ，并使，并使 轉(zhuǎn)到轉(zhuǎn)到 直至直至n N 結(jié)束。結(jié)束。 10,xx11, yx0101,yyxxSchool of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差

43、分方法 (2) 四階龍格四階龍格庫(kù)塔算法流程圖庫(kù)塔算法流程圖 開 始 輸 入 x0, y0,h , N 1 n x0 + h x1 f(x0,y0 ) k1, f(x0+ h /2 ,y0 + h k1/2 ) k2 f(x0+ h /2 ,y0 + h k2/2 ) k3, f(x1,y0+ h k3) k4 y0+ h (k1+ 2 k2+ + 2 k3+ k4)/6 y1 輸 出 x1, y1 n + 1 n n = N ? x1 x0 y1 y0 結(jié) 束 n y School of Automation Engineering自自 動(dòng)動(dòng) 化化 工工 程程 學(xué)學(xué) 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 例例7.4 取步長(zhǎng)取步長(zhǎng)h=0.2，用經(jīng)典格式求解初值問題，用經(jīng)典格式求解初值問題 1)0(2yxyy10 x解解: 由四階龍格由四階龍格-庫(kù)塔公式可得庫(kù)塔公式可得 2 . 0, 1, 0,2),(00hyxxyyxf0),(001yxfk2 . 0) 1 , 1 . 0()2,(10202fkhyxfkh204. 0)02. 1 , 1 . 0()2,(20203fkhyxfkh4