第九章常微方程數(shù)值解法

第九章常微方程數(shù)值解法第章9章常微分方程數(shù)值解法8-1第1頁，課件共31頁，創(chuàng)作于2023年2月第9章目錄§1歐拉（Euler）方法

1.1Euler法及其簡單改進(jìn)1.2改進(jìn)的Euler法§2龍格庫塔（Runge-kutta）方法

2.1龍格-庫塔方法的基本思想2.2二階龍格-庫塔公式2.3高階R-K公式2.4變步長R-K法§3線性多步法§4一階方程組與高階方程初值問題§5收斂性與穩(wěn)定性第2頁，課件共31頁，創(chuàng)作于2023年2月第8章序

許多科學(xué)技術(shù)問題，例如天文學(xué)中的星體運(yùn)動，空間技術(shù)中的物體飛行，自動控制中的系統(tǒng)分析，力學(xué)中的振動，工程問題中的電路分析等，都可歸結(jié)為常微分方程的初值問題。

所謂初值問題，是函數(shù)及其必要的導(dǎo)數(shù)在積分的起始點(diǎn)為已知的一類問題，一般形式為：

我們先介紹簡單的一階問題：第3頁，課件共31頁，創(chuàng)作于2023年2月第9章序由常微分方程的理論可知：上述問題的解唯一存在。

常微分方程求解求什么？應(yīng)求一滿足初值問題（8—1）的解函數(shù)y=y(x)，如對下列微分方程：《高等數(shù)學(xué)》中，微分方程求解，如對一階微分方程：y

=f(x,y)是求解解函數(shù)y=y(x)，使?jié)M足上述方程。但能夠求出準(zhǔn)確的解析函數(shù)y(x)的微分方程是很少的，《高數(shù)》中研究微分方程的求解，是分門別類討論，對不同類型的微分方程，求解方法不一樣，因此，要求解微分方程，首先必須認(rèn)清類型。第4頁，課件共31頁，創(chuàng)作于2023年2月微分方程數(shù)值解而常微分方程初值問題的數(shù)值解法，是要尋求解函數(shù)y(x)在一系列點(diǎn)y(xi)（離散點(diǎn)）:上y(xi)的近似值yi（i=1,2,…,n），并且還可由這些（xi,yi)（i=1,2,…,n）構(gòu)造插值函數(shù)作為近似函數(shù)。上述離散點(diǎn)相鄰兩點(diǎn)間的距離hi=xi-1-xi稱為步長，若hi都相等為一定數(shù)h,則稱為定步長，否則為變步長。由于在實(shí)際問題和科學(xué)研究中遇到的微分方程往往很復(fù)雜，絕大多數(shù)很難，甚至不可能求出解析函數(shù)y(x)，因此只能考慮求其數(shù)值解。本章重點(diǎn)討論如下一階微分方程：在此基礎(chǔ)上介紹一階微分方程組與高階微分方程的數(shù)值解法。第5頁，課件共31頁，創(chuàng)作于2023年2月§1歐拉（Euler）法

以Euler法及其改進(jìn)方法為例,說明常微分方程初值問題數(shù)值解法的一般概念，Euler法很簡單，準(zhǔn)確度也不高，介紹此方法的目的，是由于對它的分析討論能夠比較清楚地顯示出方法的一些特點(diǎn)，而這些特點(diǎn)及基本方法反映了其它方法的特點(diǎn)。Euler法用于求解一階微分方程初值問題：第6頁，課件共31頁，創(chuàng)作于2023年2月1.1Euler法及其簡單改進(jìn)

Euler公式為：由x0出發(fā)

x1,x2,…,xN，而利用此式可算出對應(yīng)的y1,y2,…,yN，式（8-2）稱為差分方程（序列{yn}滿足的方程）。下面是Euler公式的推導(dǎo)：一、從幾何意義出發(fā)：y

=f(x,y)的解函數(shù)y=y(x)在xoy平面上是一族解曲線，而初值問題則是其中一條積分曲線假定y=y(x)的曲線如圖8-1從給定的點(diǎn)P0(x0,y0)出發(fā)，以P0為切點(diǎn)，作切線，切線斜率為曲線y(x)的切線斜率

y

=f(x0,y0)，因此可得切線：（點(diǎn)斜式）P1P2y(x)P0

x2x1x0第7頁，課件共31頁，創(chuàng)作于2023年2月Euler公式的推導(dǎo)(續(xù)1）幾何意義：用折線近似解曲線y=y(x)，折線不會偏離太遠(yuǎn)，因?yàn)槊宽?xiàng)以f(x,y)（斜率）修正。切線與x=x1交于P1(x1,y1)，在[x0,x1]上以切線近似曲線，

第8頁，課件共31頁，創(chuàng)作于2023年2月Euler公式的推導(dǎo)（續(xù)2）二、利用Taylor級數(shù)：將y(x)在xn處展開：

第9頁，課件共31頁，創(chuàng)作于2023年2月Euler公式的推導(dǎo)（續(xù)4）第10頁，課件共31頁，創(chuàng)作于2023年2月Euler公式的推導(dǎo)（續(xù)5）四、利用數(shù)值積分公式：在[xn,xn+1]上對y

(x)=f(x,y(x))積分

對右端積分項(xiàng)采用不同的數(shù)值積分公式，便可得到各種不同的求解dE初值問題的計算公式。如，以矩形面積代替曲邊梯形面積1）以左矩形面積代替曲邊梯形面積如圖8-2，亦即以yf(x,y)xnxxn+1圖8-2

第11頁，課件共31頁，創(chuàng)作于2023年2月yf(x,y)xnxxn+1圖8-3yf(x,y)xnxxn+1圖8-43）以梯形公式（面積）代替曲邊形如圖8-4則有

式（8-5）稱為求dE初值問題的梯形公式，梯形公式看作是以（xn,yn）(xn+1,yn+1)構(gòu)造的插值多項(xiàng)式代替被積函數(shù)得到的，而Euler公式則是以左端點(diǎn)函數(shù)值近似被積函數(shù)而得到，還可以用多個點(diǎn)做插值多項(xiàng)式近似被積函數(shù)構(gòu)造另一些精度更高的解微分方程的數(shù)值公式，梯形公式比Euler公式更準(zhǔn)確一些，誤差更小。Euler公式的推導(dǎo)（續(xù)6）2）以右矩形面積代替曲邊梯形(后退的梯形公式）：如圖8-3

第12頁，課件共31頁，創(chuàng)作于2023年2月Euler公式注釋注1：Euler公式為顯式，后退的Euler公式，梯形公式為隱式；注2：Euler公式，梯形，后退的Euler公式為單步法，計算yn+1只用yn，而中點(diǎn)法公式為多步法（還可如上二所述，構(gòu)造多步法）即必須已知yn-1，yn才能計算yn+1，（求y0,y1不能用此公式。y0,y1稱為多步法的開始值，y0給定，而y1必須由其它公式算出，然后才能用中點(diǎn)法）；注3：前面已有Euler法的局部截斷誤差：后退Euler法的局部截斷誤差：

誤差階：如果局部截斷誤差則稱方法為P階的。第13頁，課件共31頁，創(chuàng)作于2023年2月顯然，步長h越小，階數(shù)P越高，局部截斷誤差越小，當(dāng)然計算精度越高；

注4：梯形法是幾階？梯形法精度比Euler法高，階數(shù)肯定比Euler法高，其實(shí)我們可以利用數(shù)值積分公式的誤差估計式，因?yàn)槲覀兪怯锰菪螖?shù)值積分公式計算因此由積分中梯形公式的誤差知此時的局部截斷誤差為：∴梯形法為2階方法！Euler法，后退Euler法為1階方法，而中點(diǎn)法為2階，Euler公式注釋（續(xù)）第14頁，課件共31頁，創(chuàng)作于2023年2月關(guān)于Euler法的整體截斷誤差注釋注5：關(guān)于Euler法的整體截斷誤差：

實(shí)際計算時，yn是y(xn)的近似值，因此，計算過程中除每步所產(chǎn)生的局部截斷誤差外，還有因前面的計算不準(zhǔn)確而引起的誤差。在不考慮舍入誤差的情況下，稱y(xn+1)與yn+1之差為整體截斷誤差，記為：下面討論Euler方法的整體截斷誤差。

為簡便起見，假定函數(shù)f(x,y)充分光滑，問題（8-1）解y(x)在[a,b]上二階連續(xù)可微，于是由式（8-6），局部截斷誤差有界，即存在M>0，使得對任意x

[a,b]，都有|y′(x)|≤M，從而有：（緊接下屏）第15頁，課件共31頁，創(chuàng)作于2023年2月結(jié)論

對于實(shí)際問題來說，由于L，M難以估計，因很難應(yīng)用，而且上述推導(dǎo)過程中一再放大了誤差上限，這樣的估計往往也很保守，遠(yuǎn)遠(yuǎn)大于實(shí)際的誤差，但是，從估計式（8-7）卻可以得到下面很有用處的結(jié)論。

1）當(dāng)h

0時，en

0即，亦即數(shù)值解yn，一致收斂于初值問題（8-1）的真y(xn)，并且，Euler法的整體截斷誤差的階為O(h)與h同階，比局部截斷誤差低一階。2）舍入誤差局部截斷誤差對實(shí)際計算結(jié)果有影響，并且隨h減少而減少或增大。第16頁，課件共31頁，創(chuàng)作于2023年2月3）計算結(jié)果與解法的階數(shù)p，真解的導(dǎo)數(shù)y(p+1)有關(guān)，p越大，hp+1越小，|y(p+1)(ξ)|的上限越大，M也越大，因此為保證精度當(dāng)然應(yīng)選階數(shù)p較高的方法。但如果M很大，當(dāng)f(x,y)是分段連續(xù)的函數(shù)時，則應(yīng)采用低階的方法如用Euler法。結(jié)論（續(xù)）

4）計算結(jié)果還與開始值的精度有關(guān)，為使這種誤差的影響不致于超過局部截斷誤差，對多步法，應(yīng)采用跟多步法同階的方法計算開始值。第17頁，課件共31頁，創(chuàng)作于2023年2月1.2改進(jìn)的Euler法

梯形公式為二階方法，但卻是隱式格式，即若利用梯形公式求yn+1，就要求解方程（8-5）式，計算量較大，通常在實(shí)際計算時，將Euler法與梯形公式合起來使用，即先使用Euler公式,由（xn,yn）算出yn+1，記為yn+1(0)，稱為

預(yù)測值，然后用梯形公式去提高精度，將yn+1(0)

校正為較準(zhǔn)確的值：由于函數(shù)f(x,y)滿足Lipschitz條件，容易得出:第18頁，課件共31頁，創(chuàng)作于2023年2月改進(jìn)的Euler法(續(xù)）第19頁，課件共31頁，創(chuàng)作于2023年2月預(yù)測──校正型公式

實(shí)際經(jīng)驗(yàn)表明，式（8-8）的迭代效果主要體現(xiàn)在第一次，由此構(gòu)成如下的預(yù)測──校正型公式:此式稱為改進(jìn)的Euler公式，為上機(jī)計算編程方便，常將式

（8-9）改寫為:下面分析改進(jìn)的Euler公式的局部截斷誤差：

第20頁，課件共31頁，創(chuàng)作于2023年2月改進(jìn)的Euler公式的局部截斷誤差分析假定yn=y(xn)，y(xn+1)的Taylor展式為:對于改進(jìn)的Euler公式，由于

這說明改進(jìn)的Euler法的局部截斷誤差為O(h3)，比Euler公式高一階，是二階方法。第21頁，課件共31頁，創(chuàng)作于2023年2月改進(jìn)的Euler公式舉例例1

這些結(jié)果在表8-1中，可見計算結(jié)果的精度，Euler法與后退Euler法差不多，與準(zhǔn)確值相比較Euler法偏小，而后退Euler法偏大；中點(diǎn)法與梯形法精度同為2階，但梯形法更好一些，這跟它們局部截斷誤差的符號，階數(shù)和系

數(shù)的大小是完全一致的。表見下屏：第22頁，課件共31頁，創(chuàng)作于2023年2月表格8-1表8-1

y

=

y,y(0)=1的數(shù)值解(h=0.1)

x精確解歐拉法后退歐拉中點(diǎn)法梯形法.1.904837.900000.909091.900000.904762.2.808731.810000.826446.820000.818594.3.740818.729000.751315.736000.740633.4.670320.656100.683013.627800.670096.5.060531.590490.620921.601440.606278.6.548812.531441.654474.552512.548537.7.496585.478298.5131458.490938.496295.8.449329.430467.466507.454324.449029.9.406570.387421.424098.400073.4062641.367879.348679.385543.374310.367573第23頁，課件共31頁，創(chuàng)作于2023年2月表格8-2

而表8-2是分別取了不同的h=0.1,h=0.01，h=0.001,h=0.0001，還是利用這些公式，經(jīng)過若干步的計算（h越小，計算量越大）算到y(tǒng)(1)的近似值，可見：隨著h的減小，y(1)的近似值的精度在提高，0.01比0.001差，即0.001比0.01時的y(1)準(zhǔn)確。Y′=-y,y(0)=1的解y(1)的近似值(y(1)=0.367879)h歐拉法后退歐拉法中點(diǎn)法梯形法0.1.348678.385543.374310.3675730.01.366033.369711.367944.3678770.001.367700.368052.367879.3678760.0001.367800.367800.367881.368020（緊接下屏）第24頁，課件共31頁，創(chuàng)作于2023年2月表8-2計算結(jié)果說明（續(xù)）

但h太小，到h=0.0001時卻又變得誤差大了，這與前面所說h越小，p階越高，應(yīng)該局部截斷誤差越小，因而計算精度更高矛盾了，為什么會產(chǎn)生這種情況呢？這是由于h太小而引起計算量大因而造成了舍入誤差和截斷誤差的積累，這種情況由于初值問題不同可能會影響更大，偏離更嚴(yán)重，如下面的例2。這種問題實(shí)際上是穩(wěn)定性問題，我們將會討論方法的穩(wěn)定性，由此得出對h有一定的要求的穩(wěn)定性制區(qū)域。第25頁，課件共31頁，創(chuàng)作于2023年2月§5收斂性與穩(wěn)定性

通過前面的討論可以看到，微分方

程數(shù)值解法的基本思想是：

通過某種離散化手段，將微分方程

轉(zhuǎn)化為差分格式求解。從理論上說，這

里需要解決兩個問題，一是當(dāng)步長h

0

時，差分解，這就是差分格式的收斂性

問題；二是利用差分格式求解時，初始

誤差及計算過程中的舍入誤差能否得到

控制，這就是差分格式的穩(wěn)定性問題。

本節(jié)對這兩個問題作一簡單介紹。第26頁，課件共31頁，創(chuàng)作于2023年2月5.1收斂性

用某種差分格式解微分方程時，若對求解區(qū)間[a,b]中的任一點(diǎn)x，當(dāng)h

0時差分解yn