常微分方程數(shù)值解

1、1 1、常微分方程與解、常微分方程與解為為n n階常微分方程階常微分方程。0 ), , ,()(nyyyyxF如果函數(shù)如果函數(shù) 在區(qū)間在區(qū)間a,ba,b內(nèi)內(nèi)n n階可導(dǎo)，稱方程階可導(dǎo)，稱方程)(xyy )(xyy 滿足方程的函數(shù)滿足方程的函數(shù)稱為微分方程的稱為微分方程的解解。則則如如為任意常數(shù)）為任意常數(shù)）xy2 CCxy(, 2一般稱為方程的一般稱為方程的通解通解。為方程的解。為方程的解。12 xy如果如果則有則有10 )(y為方程滿足定解條件的解。為方程滿足定解條件的解。第第10章章 常微分方程的數(shù)值解常微分方程的數(shù)值解10.1引引 言言 科學(xué)研究和工程技術(shù)中的問題往往歸結(jié)為求某個(gè)常微分科

2、學(xué)研究和工程技術(shù)中的問題往往歸結(jié)為求某個(gè)常微分方程的定解問題方程的定解問題. 常微分方程的理論指出，除少數(shù)簡單情況能獲得初值問常微分方程的理論指出，除少數(shù)簡單情況能獲得初值問題的初等解（用初等函數(shù)表示的解）外，絕大多數(shù)情況下是題的初等解（用初等函數(shù)表示的解）外，絕大多數(shù)情況下是求不出初等解的求不出初等解的.有些初值問題即便有初等解，也往往由于形有些初值問題即便有初等解，也往往由于形式過于復(fù)雜而不便處理。式過于復(fù)雜而不便處理。 常微分方程的數(shù)值解法常用來求近似解，由于它提供的常微分方程的數(shù)值解法常用來求近似解，由于它提供的算法能通過計(jì)算機(jī)便捷地實(shí)現(xiàn)，因此近年來得到迅速的發(fā)展算法能通過計(jì)算機(jī)便捷地

3、實(shí)現(xiàn)，因此近年來得到迅速的發(fā)展和廣泛的應(yīng)用。和廣泛的應(yīng)用。10.2 初值問題解法的基本概念 科學(xué)技術(shù)中常常需要求解常微分方程的定解問科學(xué)技術(shù)中常常需要求解常微分方程的定解問題題. 這類問題最簡單的形式，是本章將著重考察的這類問題最簡單的形式，是本章將著重考察的一一階方程的初值問題階方程的初值問題00( , ),(2.1)().yf x yy xy 我們知道，只有我們知道，只有f(x, y)適當(dāng)光滑適當(dāng)光滑譬如關(guān)于譬如關(guān)于y滿足滿足利普希茨利普希茨(Lipschitz)條件條件理論上就可以保證初值問題的解理論上就可以保證初值問題的解yf(x)存在并且唯一存在并且唯一.我們以下的討論，都在滿足上述

4、條件下進(jìn)行。我們以下的討論，都在滿足上述條件下進(jìn)行。( , )( , ).(2.2)f x yf x yL yy 雖然求解常微分方程有各種各樣的解析方法，但雖然求解常微分方程有各種各樣的解析方法，但解析方法只能用來求解一些特殊類型的方程，實(shí)際問解析方法只能用來求解一些特殊類型的方程，實(shí)際問題中歸結(jié)出來的微分方程主要靠數(shù)值解法題中歸結(jié)出來的微分方程主要靠數(shù)值解法. 所謂所謂數(shù)值解法數(shù)值解法, 就是尋求解就是尋求解y(x)在一系列離散節(jié)在一系列離散節(jié)點(diǎn)點(diǎn) 121nnxxxx上的近似值上的近似值 y1,y2,yn,yn+1,. 相鄰兩個(gè)節(jié)點(diǎn)的間距相鄰兩個(gè)節(jié)點(diǎn)的間距hn=xn+1- -xn稱為稱為步長

5、步長. 今后如不特別說明，總是假定今后如不特別說明，總是假定 hi=h(i=1,2,)為為定數(shù)定數(shù), 這時(shí)節(jié)點(diǎn)為這時(shí)節(jié)點(diǎn)為xn=x0+nh(i=0,1,2,) (等距節(jié)點(diǎn)等距節(jié)點(diǎn)). 常微分方程數(shù)值解是一組離散的函數(shù)值數(shù)據(jù)，它常微分方程數(shù)值解是一組離散的函數(shù)值數(shù)據(jù)，它的精確表達(dá)式很難求解得到，但可以進(jìn)行插值計(jì)算后的精確表達(dá)式很難求解得到，但可以進(jìn)行插值計(jì)算后用插值函數(shù)逼近用插值函數(shù)逼近 y(x) 初值問題的初值問題的數(shù)值解法的基本特點(diǎn)：數(shù)值解法的基本特點(diǎn)： 都采取都采取“步進(jìn)式步進(jìn)式”，即求解過程順著節(jié)點(diǎn)排列的，即求解過程順著節(jié)點(diǎn)排列的次序一步一步地向前推進(jìn)次序一步一步地向前推進(jìn). 首先，要對

6、微分方程離散化，建立求解數(shù)值解的首先，要對微分方程離散化，建立求解數(shù)值解的遞推公式遞推公式. 一類是計(jì)算一類是計(jì)算yn+1時(shí)時(shí)只用到只用到xn+1，xn和和yn，即，即前一步的值。因此，有了初值以后就可以逐步往下計(jì)前一步的值。因此，有了初值以后就可以逐步往下計(jì)算，其代表是龍格庫塔法算，其代表是龍格庫塔法稱為稱為單步法單步法. 另一類是用另一類是用到到y(tǒng)n+1前面前面 k 點(diǎn)的值點(diǎn)的值yn,yn-1, yn-k+1，稱為，稱為多步法多步法. 其次，要研究公式的其次，要研究公式的局部截?cái)嗾`差局部截?cái)嗾`差和和階階，數(shù)值解，數(shù)值解yn與精確解與精確解y(xn)的的誤差估計(jì)誤差估計(jì)及及收斂性收斂性，還有

7、遞推公式，還有遞推公式的的計(jì)算穩(wěn)定性計(jì)算穩(wěn)定性等問題等問題.數(shù)值解的思想數(shù)值解的思想（1 1）將連續(xù)變量）將連續(xù)變量 離散為離散為,bax bxxxxank 10nkxyykk,)(21 （2 2）用代數(shù)的方法求出解函數(shù)）用代數(shù)的方法求出解函數(shù) 在在 點(diǎn)的近似值點(diǎn)的近似值)(xyy kx)(kxy* *ky)(xyy 數(shù)學(xué)界關(guān)注數(shù)學(xué)界關(guān)注工程師關(guān)注工程師關(guān)注如果找不到解函數(shù)如果找不到解函數(shù)數(shù)學(xué)界還關(guān)注：數(shù)學(xué)界還關(guān)注：解的存在性解的存在性解的唯一性解的唯一性解的光滑性解的光滑性解的振動(dòng)性解的振動(dòng)性解的周期性解的周期性解的穩(wěn)定性解的穩(wěn)定性解的混沌性解的混沌性第一步：連續(xù)變量離散化第一步：連續(xù)變量離

8、散化,nkxxxxx10第二步：用直線步進(jìn)第二步：用直線步進(jìn)1 1、EulerEuler格式格式10.3 簡單單步法簡單單步法10.3.1 歐拉歐拉(Euler)方法方法00(,)xy000()(,)y xf xy0000(,)()yyf xyxx1xx100010(,)()yyf xyxx11()y xy11( ,)x y111( )( ,)y xf x y1111( ,)()yyf x yxx2xx211121( ,)()yyf x yxx22()y xy(,)nnxy()(,)nnny xf xy(,)()nnnnyyf xyxx過過 做以做以 為切線斜率的方程為切線斜率的方程 當(dāng)當(dāng)時(shí)，

9、得時(shí)，得，取取 當(dāng)當(dāng)時(shí)，得時(shí)，得，取取 過過做以做以為切線斜率的方程為切線斜率的方程一般地，過一般地，過做以做以為切線斜率的方程為切線斜率的方程kp0p1p1npnpkx0 x1x1nxnx),(),(nnnnnnnnnnyxhfyyyxfxxyy 111EulerEuler格式格式nxx11(,)()nnnnnnyyf xyxx()nny xy當(dāng)當(dāng)時(shí)，得時(shí)，得，取取 (3.1) 例例1 用歐拉公式求解初值問題用歐拉公式求解初值問題2(01),(0)1.xyyxyy 解解 取步長取步長h=0.1，歐拉公式的具體形式為，歐拉公式的具體形式為)2(1nnnnnyxyhyy 其中其中xn=nh=0.

10、1n (n=0,1,10), 已知已知y0 =1, 由此式可得由此式可得191818. 1)1 . 12 . 01 . 1 ( 1 . 01 . 1)2(1 . 11 . 01)2(1111200001 yxyhyyyxyhyy依次計(jì)算下去，依次計(jì)算下去，部分計(jì)算結(jié)果部分計(jì)算結(jié)果見下表見下表. xy21 與準(zhǔn)確解與準(zhǔn)確解 相比，可看出歐拉公式的計(jì)算結(jié)相比，可看出歐拉公式的計(jì)算結(jié)果精度很差果精度很差. xn 歐拉公式數(shù)值解歐拉公式數(shù)值解yn準(zhǔn)確解準(zhǔn)確解y(xn) 誤差誤差 0.2 0.4 0.6 0.8 1.0 1.191818 1.358213 1.508966 1.649783 1.7847

11、70 1.183216 1.341641 1.483240 1.612452 1.732051 0.008602 0.016572 0.025726 0.037331 0.052719 歐拉公式具有明顯的幾何意義歐拉公式具有明顯的幾何意義, , 就是就是用折線近似用折線近似代替方程的解曲線代替方程的解曲線，因而常稱公式，因而常稱公式(3.1)為為歐拉折線法歐拉折線法. .( )yy xxynx1nxnp1np1np x 還可以通過幾何直觀來考察歐拉方法的精度還可以通過幾何直觀來考察歐拉方法的精度. .假假設(shè)設(shè)yn=y(xn), ,即頂點(diǎn)即頂點(diǎn)Pn落在積分曲線落在積分曲線y=y(x)上，那么，上

12、，那么，按歐拉方法做出的折線按歐拉方法做出的折線PnPn+1便是便是y=y(x)過點(diǎn)過點(diǎn)Pn的切線的切線. .從圖形上看從圖形上看, ,這這樣定出的頂點(diǎn)樣定出的頂點(diǎn)Pn+1顯著顯著地偏離了原來的積分曲地偏離了原來的積分曲線，可見歐拉方法是線，可見歐拉方法是相相當(dāng)粗糙當(dāng)粗糙的的. .12 方法一化導(dǎo)數(shù)為差商的方法10()()()()()limnnnnnhy xhy xy xy xy xhh 由于在逐步求解的過程中，由于在逐步求解的過程中，y(xn) 的準(zhǔn)確值無法求解的準(zhǔn)確值無法求解出來，因此用其近似值代替。出來，因此用其近似值代替。為避免混淆，以下學(xué)習(xí)簡記：為避免混淆，以下學(xué)習(xí)簡記：y(xn)：

13、待求函數(shù)：待求函數(shù) y(x) 在在 xn 處的精確函數(shù)值處的精確函數(shù)值yn ：待求函數(shù)：待求函數(shù) y(x) 在在 xn 處的近似函數(shù)值處的近似函數(shù)值歐拉歐拉(Euler)方法（幾種推導(dǎo)法）方法（幾種推導(dǎo)法）13 代入初值問題表達(dá)式可得：根據(jù)根據(jù) y0 可以一步步計(jì)算出函數(shù)可以一步步計(jì)算出函數(shù) y y(x) 在在 x1, x2, x3 x4, 上的近似值上的近似值 y1, y2, y3, y4 , 10()()()()()limnnnnnhy xhy xy xy xy xhh 100(,)0,1,2,()nnnnyyhf xynyy x 1(,)nnnnyyf xyh 為了分析計(jì)算公式的精度，通

14、?？捎锰├照归_為了分析計(jì)算公式的精度，通?？捎锰├照归_將將y(xn+1)在在xn處展開，則有處展開，則有).,()(2),()()(2)()()()(1221 nnnnnnnnnnnnxxyhyxhfxyyhxyhxyhxyxy 在在yn=y(xn)的前提下，的前提下，f(xn,yn )=f(xn,y(xn)=y ( (xn n) ). .于是于是可得歐拉法可得歐拉法(3.1)的的公式誤差公式誤差為為2211()( )( ),(3.2)22nnnnhhy xyyy x稱為此方法的稱為此方法的局部截?cái)嗾`差局部截?cái)嗾`差. . 方法二方法二泰勒級數(shù)展開法泰勒級數(shù)展開法15 方法三數(shù)值積分法數(shù)值積分法

15、1111d()() , ( )d() , (), ()nnnnxxnnxxnnnnnnyxy xy xf x y xxxxf xy xhf xy x 同樣以近似值同樣以近似值 yn 代替精確值代替精確值 y(xn) 可得可得：100(,)()nnnnyyhf xyyy x 將微分方程將微分方程 y f (x, y) 在區(qū)間在區(qū)間 xn, xn+1 上積分：上積分：162.隱式歐拉法隱式歐拉法(后退后退) 在數(shù)值積分法推導(dǎo)中，積分的近似值取為積分區(qū)間寬在數(shù)值積分法推導(dǎo)中，積分的近似值取為積分區(qū)間寬度與右端點(diǎn)處的函數(shù)值乘積，即：度與右端點(diǎn)處的函數(shù)值乘積，即：這樣便得到了隱式歐拉法：這樣便得到了隱式

16、歐拉法：11100(,)()nnnnyyh f xyyy x 含有未知含有未知的函數(shù)值的函數(shù)值11111111d()() , ( )d() , (), ()nnnnxxnnxxnnnnnnyxy xy xf x y xxxxf xy xhf xy x （3.3） 隱式歐拉公式與歐拉公式有著本質(zhì)的區(qū)別隱式歐拉公式與歐拉公式有著本質(zhì)的區(qū)別, 后者后者是關(guān)于是關(guān)于yn+1的一個(gè)直接計(jì)算公式，這類公式稱作是的一個(gè)直接計(jì)算公式，這類公式稱作是顯顯式的式的；前者公式的；前者公式的右端含有未知的右端含有未知的yn+1，它實(shí)際上是，它實(shí)際上是關(guān)于關(guān)于yn+1的一個(gè)函數(shù)方程的一個(gè)函數(shù)方程, ,這類方程稱作是這類

17、方程稱作是隱式的隱式的. . 顯式顯式與與隱式隱式兩類方法各有特點(diǎn)，考慮到數(shù)值穩(wěn)兩類方法各有特點(diǎn)，考慮到數(shù)值穩(wěn)定性等其他因素，人們有時(shí)需要選用定性等其他因素，人們有時(shí)需要選用隱式隱式方法，但方法，但使用使用顯式顯式算法遠(yuǎn)比算法遠(yuǎn)比隱式隱式方便方便. . 隱式方程通常用迭代法求解，而迭代過程的實(shí)隱式方程通常用迭代法求解，而迭代過程的實(shí)質(zhì)是質(zhì)是逐步逐步顯式化顯式化. . 設(shè)用歐拉公式設(shè)用歐拉公式),() 0(1nnnnyxhfyy 給出迭代初值給出迭代初值 ，用它代入，用它代入(3.1)式的式的右端，使之轉(zhuǎn)右端，使之轉(zhuǎn)化為顯式，直接計(jì)算得化為顯式，直接計(jì)算得) 0(1 ny),() 0(11) 1

18、 (1 nnnnyxhfyy然后再用然后再用 代入代入(3.1)式，又有式，又有) 1 (1 ny).,() 1 (11) 2(1 nnnnyxhfyy如此反復(fù)進(jìn)行，得如此反復(fù)進(jìn)行，得(1)( )111(,) (0,1, ).(3.4)kknnnnyyhf xyk由于由于f(x, y)對對y滿足滿足Lipschitz條件條件(2.1). 由由(3.4)減減(3.3)得得),(),(11)(111) 1(1 nnknnnknyxfyxfhyy.1)(1 nknyyhL由此可知，只要由此可知，只要hL n) 上產(chǎn)生的擾動(dòng)為上產(chǎn)生的擾動(dòng)為 ，如果：，如果：md d0nd d (1,2,)mnmnnd

19、ddd定義定義：設(shè)在節(jié)點(diǎn)：設(shè)在節(jié)點(diǎn) xn 處用數(shù)值算法得到的理想數(shù)值解為處用數(shù)值算法得到的理想數(shù)值解為 yn，而實(shí)際計(jì)算得到的近似解為，而實(shí)際計(jì)算得到的近似解為 ，稱差值：，稱差值：ny nnnyyd d 為第為第 n 步的數(shù)值解的步的數(shù)值解的擾動(dòng)擾動(dòng)。則稱該數(shù)值方法是則稱該數(shù)值方法是穩(wěn)定穩(wěn)定的。的。 下面以歐拉法為例考察計(jì)算穩(wěn)定性下面以歐拉法為例考察計(jì)算穩(wěn)定性. 例例4 用歐拉公式求解初值問題用歐拉公式求解初值問題 . 1)0(,100yyy 解解 用歐拉法解方程用歐拉法解方程y=- -100y 得得其準(zhǔn)確解其準(zhǔn)確解 是一個(gè)按指數(shù)曲線衰減很快的是一個(gè)按指數(shù)曲線衰減很快的函數(shù)函數(shù).xexy10

20、0)( .)1001(1nnyhy 若取步長若取步長h=0.025，則歐拉公式的具體形式為，則歐拉公式的具體形式為.5 . 11nnyy 節(jié)點(diǎn)節(jié)點(diǎn)xn歐拉方法歐拉方法yn后退歐拉方法后退歐拉方法yn0.0250.0500.0750.100- -1.5 2.25 - -3.375 5.06250.28570.08160.02330.0067計(jì)算結(jié)果見表計(jì)算結(jié)果見表, , 明顯計(jì)算過程不穩(wěn)定明顯計(jì)算過程不穩(wěn)定, , 但取但取h=0.005, , yn+1=- -1.5yn, , 則計(jì)算過程穩(wěn)定則計(jì)算過程穩(wěn)定. .對后退的歐拉公式，取對后退的歐拉公式，取h=0.025時(shí)，則計(jì)算公式時(shí)，則計(jì)算公式為為

21、yn+1=- -(1/3.5)yn . .計(jì)算結(jié)果見表計(jì)算結(jié)果見表, , 這時(shí)計(jì)算過程是這時(shí)計(jì)算過程是穩(wěn)定的穩(wěn)定的. . 例題表明穩(wěn)定性不但例題表明穩(wěn)定性不但與方法有關(guān)與方法有關(guān)，也，也與步長與步長h有有關(guān)，當(dāng)然關(guān)，當(dāng)然與方程中的與方程中的f(x, y)有關(guān)有關(guān). 為了只考察數(shù)值方為了只考察數(shù)值方法本身，通常只檢驗(yàn)數(shù)值方法用于解法本身，通常只檢驗(yàn)數(shù)值方法用于解模型方程模型方程的穩(wěn)的穩(wěn)定性，定性，模型方程模型方程為為,(4.4)yy 其中其中為已知實(shí)數(shù)或復(fù)數(shù)為已知實(shí)數(shù)或復(fù)數(shù)(Re()0) ，這個(gè)方程分析，這個(gè)方程分析較簡單，對一般方程可以通過局部線性化化為這種較簡單，對一般方程可以通過局部線性化

22、化為這種形式。形式。50 定義定義6 6 單步法（單步法（4.24.2）用于解模型方程（）用于解模型方程（4.44.4），若得），若得到的解到的解 ，滿足，滿足 ，則稱方法，則稱方法（4.14.1）是）是絕對穩(wěn)定絕對穩(wěn)定的的. . 1()nnyE hy1)(hE 在在 的平面上，使的平面上，使 的變量圍成的的變量圍成的區(qū)域，稱為區(qū)域，稱為絕對穩(wěn)定域絕對穩(wěn)定域，h1)(hE它與實(shí)軸的交稱為它與實(shí)軸的交稱為絕對穩(wěn)定區(qū)間絕對穩(wěn)定區(qū)間. . 51歐拉法：歐拉法：1(,)nnnnyyhf xy 考察模型方程：考察模型方程：(0)yy 1(1)nnyhy 即：即：假設(shè)在節(jié)點(diǎn)值假設(shè)在節(jié)點(diǎn)值 yn 上有擾動(dòng)上有擾動(dòng) d dn，在節(jié)點(diǎn)值，在節(jié)點(diǎn)值 yn 1 上有上有擾動(dòng)擾動(dòng) d dn 1，且，且 d dn 1 僅由僅由 d dn 引起（即：計(jì)算過程中引起（即：計(jì)算過程中不再引起新的誤差）不再引起新的誤差）52111(1)(1)(1)()(1)nnnnnnnnyyhyhyhyyhd d d d 歐拉法穩(wěn)定歐拉法穩(wěn)定1nndddd 11h 即：即：111h 歐拉法穩(wěn)定的條件：歐拉法穩(wěn)定的條件：20h 0 1(1)nnyhy 針對模型方程：針對模型方程：的顯式歐拉法：的顯式歐拉法：1(,)nnnnyyhf xy 化簡得