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

1、常微分方程的數(shù)值解法第1頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三 實(shí)際中，很多問題的數(shù)學(xué)模型都是微分方程。我們可以研究它們的一些性質(zhì)。但是，只有極少數(shù)特殊的方程有解析解。對于絕大部分的微分方程是沒有解析解的。 常微分方程作為微分方程的基本類型之一，在自然界與工程界有很廣泛的應(yīng)用。很多問題的數(shù)學(xué)表述都可以歸結(jié)為常微分方程的定解問題。很多偏微分方程問題，也可以化為常微分方程問題來近似求解。 本章討論常微分方程的數(shù)值解法 引 言第2頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三 本章討論一階常微分方程的初值問題9.19.2雖然求解常微分方程有各種各樣的解析方法，但解析方法只

2、能用來求解一些特殊類型的方程，大量從實(shí)際問題當(dāng)中歸結(jié)出來的微分方程主要靠數(shù)值解法。第3頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三定義：所謂數(shù)值解法，就是尋求初值問題上的近似值相鄰兩個(gè)節(jié)點(diǎn)間的距離稱為步長。 今后如不特別申明，總假定步長h為定數(shù)。第4頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三一、幾何解釋：圖91 歐拉法的幾何解釋9.1 歐拉（Euler）方法第5頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三二、計(jì)算格式：1、公式推導(dǎo)：第6頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三2、幾何意義：第7頁，共40頁，2022年，5月20日，5點(diǎn)22分

3、，星期三3、計(jì)算格式：9.3 歐拉（Euler）法（也叫歐拉折線法）是最古老的一種數(shù)值解法，它體現(xiàn)了數(shù)值方法的基本思想，但精度很低，單獨(dú)用它來作計(jì)算往往不能滿足精度要求。第8頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三9.2 改進(jìn)的歐拉方法 同一種計(jì)算格式往往可以通過多種途徑構(gòu)造出來，本節(jié)與下一節(jié)就會(huì)看到這一點(diǎn)。一、計(jì)算格式：1、公式推導(dǎo)：將方程（9.1）的兩端同時(shí)積分，9.4第9頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三 選擇不同的近似方法計(jì)算這個(gè)積分項(xiàng)會(huì)得到不同的計(jì)算格式。例如：用矩形公式計(jì)算積分項(xiàng)代入（9.4）得第10頁，共40頁，2022年，5月20日，5點(diǎn)2

4、2分，星期三這樣建立起來的格式就是歐拉法的計(jì)算格式（9.3）。用矩形公式求積分值很粗糙，故歐拉格式精度也很低。為了改進(jìn)精度，我們改用梯形法計(jì)算左端積分第11頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三將其代入（9.4）得9.5（9.5）式被稱為解常微分方程的梯形法則。 格式（9.3）與 （9.5）有本質(zhì)上的區(qū)別，歐拉格式（9.3）是個(gè)直接的計(jì)算公式，這類格式稱作顯式的。第12頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三這個(gè)方程可以用迭代法求解（參看第五章），不過計(jì)算量比較大。2、預(yù)報(bào)校正系統(tǒng)：綜合使用上述兩種格式，先用歐拉格式，求得一個(gè)稱為預(yù)報(bào)值。這樣建立起來的預(yù)報(bào)校正

5、系統(tǒng)稱為改進(jìn)的歐拉格式。第13頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三3、改進(jìn)的歐拉格式：9.6格式（9.6）的每一步需要兩次調(diào)用函數(shù)f，它可以改寫成下列形式：第14頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三二、算法與流程圖：1、算法分析：歐拉法每一步只需對f調(diào)用一次，而改進(jìn)的歐拉法則不然，需對f調(diào)用兩次，其計(jì)算量比歐拉法增加一倍，付出這種代價(jià)的目的是為了提高精度。由此可見，它比歐拉格式的截?cái)嗾`差提高了一倍。第15頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三2、流程圖：（略）3、C源程序：#include #include #define H 0.1#

6、define N 10float f(x,y)float x,y; return(y-2*x/y); 第16頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三main() float x0=0; float y0=1; float x1,y1; float yp,yc; float h=H; int i; for(i=1;i=N;i+) x1=x0+h; yp=y0+h*f(x0,y0); yc=y0+h*f(x1,yp); y1=(yp+yc)/2; printf(x=%f,y=%fn,x1,y1); x0=x1; y0=y1; 第17頁，共40頁，2022年，5月20日，5點(diǎn)22分，

7、星期三例解解初值問題我們分別用兩種格式進(jìn)行計(jì)算，這里歐拉格式的具體形式是 第18頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三而改進(jìn)的歐拉格式是計(jì)算結(jié)果見下表：第19頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三同準(zhǔn)確解比較，第二列歐拉格式的結(jié)果大致只有兩位有效數(shù)字，而第三列改進(jìn)的歐拉格式的結(jié)果則有三位有效數(shù)字。第20頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三結(jié)點(diǎn)歐拉法改進(jìn)歐拉法準(zhǔn)確解00.10.20.30.40.50.60.70.80.91.011.11.1918181.2774381.3582131.4351331.5089661.5803381.6497

8、831.7177791.7847711.0959091.1840971.2662011.343361.4164021.4859561.5525141.6164751.6781661.73786711.0954451.1832161.2649111.3416411.4142141.4832401.5491931.6124521.6733201.732051第21頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三9.3 龍格庫塔(Runge-kutta)方法 最常用的是四階的龍格庫塔格式，但推導(dǎo)極為繁瑣，我們以二階為例，說明其思想方法。一、二階龍格庫塔法：1、基本思想：設(shè)初值問題：第22頁，

9、共40頁，2022年，5月20日，5點(diǎn)22分，星期三對差商9.7第23頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三平均斜率。由此得知，只要對平均斜率提供一種算法，由（9.7）式便相應(yīng)地得到一種計(jì)算格式。 歐拉格式： 改進(jìn)的歐拉格式：由于僅取一個(gè)點(diǎn)，所以精度很低。第24頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三 這就是龍格庫塔方法的基本思想1、二階龍格庫塔法：第25頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三（同改進(jìn)歐拉法）這樣設(shè)計(jì)出的計(jì)算格式：第26頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三9.8我們希望適當(dāng)選擇參數(shù)的值，第27頁，共40頁，

10、2022年，5月20日，5點(diǎn)22分，星期三代入（9.8）知和二階泰勒展開式第28頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三比較系數(shù)即可發(fā)現(xiàn)，要使（9.8）的截?cái)嗾`差為只要成立下列條件：9.9這里共有三個(gè)參數(shù)，但滿足兩個(gè)條件，因此有一個(gè)自由度。滿足條件（9.9）的一族格式（9.8）統(tǒng)稱二階龍格庫塔格式。這時(shí)龍格庫塔格式稱作變形的歐拉格式，其形式是：第29頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三9.10二、三階龍格庫塔格式：第30頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三仍用（9.8）所給的形式可以使上述格式（9.10）的截?cái)嗾`差為這類格式統(tǒng)稱為三階龍格

11、庫塔格式。第31頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三三、四階龍格庫塔法：繼續(xù)上述過程，可以進(jìn)一步討論四階龍格庫塔格式。一種最常用的經(jīng)典龍格庫塔格式為9.11第32頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三 若使四階龍格庫塔格式與改進(jìn)歐拉格式的總體計(jì)算量相同，可以取較大的步長，但計(jì)算精度比改進(jìn)歐拉法高很多。四、經(jīng)典龍格庫塔格式算法與流程圖： 四階龍格庫塔格式（9.11）的每一步需要四次調(diào)用函數(shù)2、流程圖：（略）3、C-源程序：1、算法分析：第33頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三#include #include #define H 0.2

12、#define N 5float f(x,y)float x,y; return(y-2*x/y); main() float x0=0; float y0=1; float x1,y1,k1,k2,k3,k4; float h=H; int i; for(i=1;i=N;i+)第34頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三 x1=x0+h; k1=f(x0,y0); k2=f(x0+h/2,y0+(h/2)*k1); k3=f(x0+h/2,y0+(h/2)*k2); k4=f(x1,y0+h*k3); y1=y0+(h/6)*(k1+2*k2+2*k3+k4); print

13、f(x=%f,y=%fn,x1,y1); x0=x1; y0=y1; 第35頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三例用四階經(jīng)典龍格庫塔格式解初值問題4、Mathematica求解函數(shù)：NDSolveeqns, y, x, xmin, xmax函數(shù)功能：對常微分方程或方程組eqns，求函數(shù)y關(guān)于x在xmin, xmax范圍內(nèi)的數(shù)值解。第36頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三解由公式（9.11），有第37頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三例用四階經(jīng)典龍格庫塔格式解初值問題第38頁，共40頁，2022年，5月20日，5點(diǎn)22分，星期三解計(jì)算結(jié)果如下表：結(jié) 點(diǎn)改進(jìn)歐拉法龍格庫塔法準(zhǔn)確解01110.21.134096 1.1832291.1832160.41.343360 1.3416671.3416410.6 1.485956 1.4832811.4832400.11.6164761.612513