實(shí)驗(yàn)名稱微分方程數(shù)值解

1、探索實(shí)驗(yàn)8 常微分方程初值問題數(shù)值解一、 實(shí)驗(yàn)?zāi)康?了解常微分方程初值問題的數(shù)值解概念，掌握解常微分方程初值問題的Euler方法及改進(jìn)的Euler方法和Runge-Kutta方法解常微分方程初值問題的算法構(gòu)造和計(jì)算。能用程序?qū)崿F(xiàn)解常微分方程初值問題的Euler方法、改進(jìn)的Euler方法和經(jīng)典的Runge-Kutta方法，學(xué)習(xí)用計(jì)算機(jī)求常微分方程初值問題數(shù)值解的一些科學(xué)計(jì)算方法和簡(jiǎn)單的編程技術(shù)。二、概念與結(jié)論1. 常微分方程初值問題： 常微分方程特解問題稱為初值問題，通常其形式為 2常微分方程初值問題數(shù)值解： 常微分方程初值問題的解y(x)在a,b上的有限個(gè)值y(xk)的近似值yk稱為常微分方程

2、初值問題數(shù)值解，其中 xk= xk-1 + hk-1 ，k=1,2,3,N。xk稱為節(jié)點(diǎn)， hk 稱為步長(zhǎng)。 通常，步長(zhǎng)h不變，取為等距步長(zhǎng) hk=h=（b-a）/N，N為等分區(qū)間a,b分割數(shù)。3.常微分方程初值問題數(shù)值解法： 求常微分方程初值問題數(shù)值解yk的方法稱為微分方程初值問題數(shù)值解法。4.單步法： 在計(jì)算yk+1之時(shí)只用到y(tǒng)k 的方法，其計(jì)算公式有： 顯式單步法計(jì)算公式: yk+1=yk+h(xk,yk;h） 隱式單步法計(jì)算公式: yk+1 = yk +h(xk ,yk,yk+1,h） 式中函數(shù)是連續(xù)函數(shù)，稱為增量函數(shù)。5.多步法： 在計(jì)算yk+1之時(shí)不僅用到y(tǒng)k ，還要用yk-1，y

3、k-2，,一般m 步法要用到y(tǒng)k,yk-1,yk-2， yk-m+1， 多步法也有顯式方法和隱式方法之分。6.數(shù)值解法的局部截?cái)嗾`差： 假設(shè)某常微分方程初值問題數(shù)值解法在x= xk沒有誤差，即y(xk)= yk，稱Tk+1 =y(xk+1)- yk+1為該數(shù)值方法的局部截?cái)嗾`差。 顯式單步法有其局部截?cái)嗾`差為：Tk+1 =y(xk+1)- y(xk)- h(xk,y(xk),h）7.數(shù)值解法的階： 如果某常微分方程初值問題數(shù)值解法的局部截?cái)嗾`差為：Tk+1 =O(h p+1)則稱該數(shù)值方法的階為P,或稱該數(shù)值方法為P階方法。 數(shù)值方法的階越高，方法越好。三、程序中Mathematica語句解釋

4、： 1.k+,+k,k-,-kk+ 表示賦值關(guān)系 k = k+1 , 如: k=1;Table+k,5獲得表2,3,4,5,6+k 表示先處理k的值,再做賦值 k=k+1, 如: k=1;Tablek+,5獲得表1,2,3,4,5k- 表示賦值關(guān)系 k = k-1, 如: k=1;Tablek-,5獲得表1, 0, -1, -2, -3-k 表示先處理k的值,再做賦值 k=k-1,如:k=1;Table-k,4獲得表0,-1,-2,-32. x+=k, x*=kx+=k 表示 x = x + kx*=k 表示 x = x * k3.Forstat,test,incr,body 以stat為初值

5、,重復(fù)計(jì)算incr和body直到test為False終止 。這里start為初始值,test為條件,incr為循環(huán)變量修正式,body為循環(huán)體,通常由incr項(xiàng)控制test的變化。 注意: 上述命令形式中的start可以是由復(fù)合表達(dá)式提供的多個(gè)初值,如果循環(huán)體生成 Break 語句,則退出For循環(huán); 如果循環(huán)體生成Continue 語句,則由incr的增量進(jìn)入For語句的下一次循環(huán)。 四、方法與程序在實(shí)際問題中，常常會(huì)遇到一些常微分方程初值問題。對(duì)這些問題如果只求助于高等數(shù)學(xué)中介紹的用求通解加確定邊界條件的方法去求解往往是無能為力的。因?yàn)橐环矫嫱ń饪赡芮蟛怀鰜?，另一方面所求的解可能是不能用?/p>

6、等函數(shù)表達(dá)的形式。人們處理這類問題主要采用常微分方程初值問題數(shù)值解的方法來求解。這類方法有很多，這里主要介紹解常微分方程初值問題的Euler方法、改進(jìn)的Euler方法和Runge-Kutta方法的構(gòu)造過程和算法程序。1. Euler方法Euler方法是最簡(jiǎn)單的求微分方程數(shù)值解的方法。這個(gè)方法由于精度不高，實(shí)用中較少使用。人們常用Euler方法來說明求微分方程數(shù)值解涉及到的一些問題。1.1 Euler方法的構(gòu)造過程： Euler方法涉及的計(jì)算公式有很多方法，這里介紹在求微分方程數(shù)值解用的較多的Taylor展開法。 設(shè) f(x,y) 充分光滑, xn+1=xn+h,將y(xn+1)在x n點(diǎn)作Ta

7、ylor展開，得:y (x n+1)=y(xn)+hy¢(xn)+(h2/2!)y²(xn)+O(h3)取其關(guān)于h 的線性部分，有：y (x n+1)»y(xn)+hy¢(xn)注意到y(tǒng)¢(xn)= f(xn,y(xn),用yk代替y(xk)并將“»”換為等號(hào)“=”,得到Euler公式：y n+1= yn+h f(xn,yn) 于是，由初始條件y0=y(a),借助Euler公式就可以依次計(jì)算出微分方程初值問題的數(shù)值解：y1 ,y2 ,yk稱由Euler公式求數(shù)值解的方法為Euler方法。顯然Euler方法是單步顯式方法。 因?yàn)閷?duì)Eul

8、er公式有局部截?cái)嗾`差為Tn+1= O(h2)因此Euler方法是一階方法。1.2 Euler方法算法： 1 輸入函數(shù)f(x,y)、初值y0、自變量區(qū)間端點(diǎn)a,b步長(zhǎng)h2 計(jì)算節(jié)點(diǎn)數(shù)n和節(jié)點(diǎn)x k 3 用Euler公式y(tǒng) n+1= yn+h f(xn,yn) 求數(shù)值解1.3 Euler方法程序：Clearx,y,ffx_,y_= Input"函數(shù)f(x,y)="y0=Input"初值y0 ="a=Input"左端點(diǎn)a="b=Input"右端點(diǎn)b="h=Input"步長(zhǎng)h="n=(b-a)/h;F

9、ori=0,i<n,i+,xk=a+i*h;y1=y0+h*fxk,y0;Print"y(",xk+h/N,")=",y1/N;y0=y1說明：本程序用Euler公式求常微方程初值問題數(shù)值解。程序執(zhí)行后，按要求通過鍵盤依次輸入輸入函數(shù)f(x,y)、初值y0、自變量區(qū)間端點(diǎn)a,b步長(zhǎng)h后，計(jì)算機(jī)則給出常微方程初值問題數(shù)值解。程序中變量說明：fx,y: 存放函數(shù)f(x,y)y0: 存放初值y0及數(shù)值解a：存放自變量區(qū)間左端點(diǎn)b: 存放自變量區(qū)間右端點(diǎn)n: 存放節(jié)點(diǎn)個(gè)數(shù)h: 存放節(jié)點(diǎn)步長(zhǎng)xk：存放節(jié)點(diǎn)xiy1: 存放數(shù)值解注：1）語句Print"

10、;y(",xk+h/N,")=",y1/N是將數(shù)值解用6位有效數(shù)字顯示出來，如果要顯示n位數(shù)有效數(shù)字，可以將語句改為Print"y(",xk+h/N,")=",Ny1,n2）Mathematica中有求微分方程初值問題數(shù)值解的命令，形式為:NDSolve yx=fx,yx,ya=y0, y, x, a, b 由NDSolve命令得到的解是以未知函數(shù)名->InterpolatingFunctionrange, <>的形式給出的,其中的InterpolatingFunctionrange, <>是所

11、求的插值函數(shù)表示的數(shù)值解, range就是所求數(shù)值解的自變量范圍。1.4 例題與實(shí)驗(yàn)例1. 用Euler方法求初值問題的數(shù)值解。取步長(zhǎng)h=0.1,并在一個(gè)坐標(biāo)系中畫出數(shù)值解與準(zhǔn)確解的圖形。解：執(zhí)行Euler方法程序后，在輸入的窗口中按提示分別輸入x-2x/y、1、0、1、0.1，每次輸入后用鼠標(biāo)點(diǎn)擊窗口的“OK”按扭，計(jì)算機(jī)在屏幕上給出如下數(shù)值解結(jié)果：y(0.1)=1.1y(0.2)=1.19182y(0.3)=1.27744y(0.4)=1.35821y(0.5)=1.43513y(0.6)=1.50897y(0.7)=1.58034y(0.8)=1.64978y(0.9)=1.71778y

12、(1.)=1.78477本題的準(zhǔn)確解為y(x)= (1 + 2 x)1/2，在一個(gè)坐標(biāo)系中畫出數(shù)值解與準(zhǔn)確解的圖形如下：圖中點(diǎn)圖是數(shù)值解，曲線為準(zhǔn)確解。從圖中可以知道數(shù)值解與準(zhǔn)確解比較接近，但還是有誤差的。2. 改進(jìn)的Euler方法 改進(jìn)的Euler方法比Euler方法精度高，其中把微分方程初值問題數(shù)值隱式解法計(jì)算公式變?yōu)轭A(yù)測(cè)、校正公式的做法很有代表性。2.1改進(jìn)的Euler方法構(gòu)造過程： 將微分方程y¢=f(x,y) 在xk,xk+1積分，得對(duì)右端的定積分用數(shù)值積分的梯形公式，有用yk代替y(xk)并將“»”換為等號(hào)“=”可得求微分方程初值問題的數(shù)值解的梯形公式：這個(gè)公式

13、是單步隱式公式?？梢宰C明它是二階方法，比Euler方法階高。不過，在給定初始條件y0=y(a)后,要求出數(shù)值解，每一次計(jì)算yk的值都要進(jìn)行非線性方程求根的迭代解法來完成，因此計(jì)算量很大。為減少計(jì)算量，通常采用迭代一兩次就轉(zhuǎn)入下一步計(jì)算的方法，特別，如果先用Euler公式計(jì)算一次以對(duì)要計(jì)算的下一步值解進(jìn)行預(yù)測(cè)，然后再用梯形公式對(duì)其進(jìn)行校正的方法得到下一步的值，它的計(jì)算格式為：這個(gè)計(jì)算格式稱為改進(jìn)的Euler公式，而稱由如上計(jì)算格式求數(shù)值解的方法為改進(jìn)的Euler方法。2.2 改進(jìn)的Euler方法算法：1.輸入函數(shù)f(x,y)、初值y0、自變量區(qū)間端點(diǎn)a,b步長(zhǎng)h2.計(jì)算節(jié)點(diǎn)數(shù)n和節(jié)點(diǎn)x k 3.

14、用改進(jìn)的Euler公式求數(shù)值解2.3 改進(jìn)的Euler方法程序：Clearx,y,ffx_,y_= Input"函數(shù)f(x,y)="y0=Input"初值y0 ="a=Input"左端點(diǎn)a="b=Input"右端點(diǎn)b="h=Input"步長(zhǎng)h="n=(b-a)/h;Fori=0,i<n,i+,xk=a+i*h;y1=y0+h*fxk,y0;y1=y0+h/2*(fxk,y0+fxk+h,y1);Print"y(",xk+h/N,")=",y1/N;y

15、0=y1注：語句h=Input"步長(zhǎng)h="的步長(zhǎng)輸入應(yīng)該用帶小數(shù)點(diǎn)的數(shù)輸入，這樣可以加快計(jì)算速度，否則，由于Mathematica軟件本身的原因，可能會(huì)出現(xiàn)計(jì)算時(shí)間過長(zhǎng)等計(jì)算不出結(jié)果的情況。說明：本程序用改進(jìn)的Euler公式求常微方程初值問題數(shù)值解。程序執(zhí)行后，按要求通過鍵盤依次輸入輸入函數(shù)f(x,y)、初值y0、自變量區(qū)間端點(diǎn)a,b步長(zhǎng)h后，計(jì)算機(jī)則給出常微方程初值問題數(shù)值解。程序中變量說明：fx,y: 存放函數(shù)f(x,y)y0: 存放初值y0及數(shù)值解a：存放自變量區(qū)間左端點(diǎn)b: 存放自變量區(qū)間右端點(diǎn)n: 存放節(jié)點(diǎn)個(gè)數(shù)h: 存放節(jié)點(diǎn)步長(zhǎng)xk：存放節(jié)點(diǎn)xiy1: 存放數(shù)值解

16、2.5例題與實(shí)驗(yàn)例2. 用改進(jìn)的Euler方法求初值問題的數(shù)值解。取步長(zhǎng)h=0.1和0.01計(jì)算，并與準(zhǔn)確解比較。解：執(zhí)行改進(jìn)的Euler方法程序后，在輸入的窗口中按提示分別輸入-2x*y2、1、0、1、0.1，每次輸入后用鼠標(biāo)點(diǎn)擊窗口的“OK”按扭，計(jì)算機(jī)在屏幕上給出如下數(shù)值解結(jié)果：y(0.1)=0.99 y(0.2)=0.961366 y(0.3)=0.917246 y(0.4)=0.861954y(0.5)=0.800034 y(0.6)=0.735527 y(0.7)=0.671587 y(0.8)=0.610399y(0.9)=0.553289 y(1.)=0.500919再執(zhí)行一次

17、改進(jìn)的Euler方法程序后，在輸入的窗口中按提示分別輸入-2x*y2、1、0、1、0.05，每次輸入后用鼠標(biāo)點(diǎn)擊窗口的“OK”按扭，計(jì)算機(jī)在屏幕上給出如下數(shù)值解結(jié)果：y(0.05)=0.9975 y(0.1)=0.990087 y(0.15)=0.977978 y(0.2)=0.961519y(0.25)=0.941158 y(0.3)=0.917417 y(0.35)=0.890863 y(0.4)=0.862076y(0.45)=0.831624 y(0.5)=0.800043 y(0.55)=0.767819 y(0.6)=0.735383y(0.65)=0.7031 y(0.7)=0.

18、671277 y(0.75)=0.640158 y(0.8)=0.609935y(0.85)=0.580748 y(0.9)=0.552699 y(0.95)=0.52585 y(1.)=0.500236本題的準(zhǔn)確解為y(x)=(1+x2)-1它在0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0的值為y(0.1)=0.990099, y(0.2)=0.961538, y(0.3)=0.917431, y(0.4)=0.862069,y(0.5)=0.8, y(0.6)=0.735294, y(0.7)=0.671141, y(0.8)=0.609756, y(0.

19、9)= 0.552486, y(1.0)= 0.5,比較上面的計(jì)算結(jié)果，顯然，步長(zhǎng)為0.05的結(jié)果比步長(zhǎng)為0.1的結(jié)果要好。3.Runge-Kutta方法Runge-Kutta方法是顯式單步高階的求微分方程數(shù)值解的方法，其中四階Runge-Kutta方法最為常用。四階Runge-Kutta方法又稱為經(jīng)典的Runge-Kutta方法，它具有精度高、程序簡(jiǎn)單，容易編程、計(jì)算過程穩(wěn)定餓易于調(diào)節(jié)步長(zhǎng)的特點(diǎn)。3.1 Runge-Kutta方法的構(gòu)造過程： 將將微分方程y¢=f(x,y) 在xk,xk+1積分并利用積分中值定理，有：為提供增加求解的精度，把f(x,y(x)寫成一個(gè)線性組合的形式并

20、用yk代替y(xk)，就得到Runge-Kutta方法的一般形式：上式中選擇不同的m,ci,x i值就得到不同的Runge-Kutta計(jì)算公式。通常，為方便獲得Runge-Kutta計(jì)算公式，常把Runge-Kutta方法的一般形式寫為：于是，利用二元Taylor展開將公式中的f在(xk,yk)展開并適當(dāng)選擇參數(shù)就可以得到具體的Runge-Kutta計(jì)算公式了。稱由Runge-Kutta計(jì)算公式式求數(shù)值解的方法為Runge-Kutta方法。下面以m=2的二階Runge-Kutta計(jì)算公式為例來說明獲得Runge-Kutta計(jì)算公式的方法。 m=2的二階Runge-Kutta計(jì)算公式為它的增量函

21、數(shù)為j (x,y(x),h)=c1f(x,y(x)+ c2f(x+a2h,y(x)+hb21f(x,y(x)，考慮它的局部截?cái)嗾`差Tk+1 =y(xk+1)- y(xk)- hj (xk,y(xk),h）。利用f(x+a2h,y(x)+hb21f(x,y(x) 在(x,y)做二元Taylor展開使其具階為2，則有：這是有四個(gè)參數(shù)三個(gè)方程的方程組，它有無窮多組解。例如，取c1=0,可以得到c2=1，a2=0.5, b21=0.5,于是得到一個(gè)m=2的二階Runge-Kutta計(jì)算公式：它稱為中點(diǎn)公式。類似地，可以得到很多其他Runge-Kutta計(jì)算公式。經(jīng)典的Runge-Kutta計(jì)算公式是四

22、階的，形式為：于是，由初始條件y0=y(a),借助經(jīng)典的Runge-Kutta計(jì)算公式就可以依次計(jì)算出微分方程初值問題的數(shù)值解：y1 ,y2 ,yn。3.2 Runge-Kutta方法算法： 4 輸入函數(shù)f(x,y)、初值y0、自變量區(qū)間端點(diǎn)a,b步長(zhǎng)h5 計(jì)算節(jié)點(diǎn)數(shù)n和節(jié)點(diǎn)x k 6 用Runge-Kutta r公式求數(shù)值解3.3 Runge-Kutta方法程序：Clearx,y,ffx_,y_= Input"函數(shù)f(x,y)="y0=Input"初值y0 ="a=Input"左端點(diǎn)a="b=Input"右端點(diǎn)b=&quo

23、t;h=Input"步長(zhǎng)h="n=(b-a)/h;Fori=0,i<10,i+,xk=a+i*h;k1=fxk,y0;k2=fxk+h/2,y0+k1*h/2;k3=fxk+h/2,y0+k2*h/2;k4=fxk+h,y0+k3*h;y1=y0+h/6*(k1+2k2+2k3+k4);Print"y(",xk+h/N,")=",y1/N;y0=y1說明：本程序用經(jīng)典Runge-Kutta公式求常微方程初值問題數(shù)值解。程序執(zhí)行后，按要求通過鍵盤依次輸入輸入函數(shù)f(x,y)、初值y0、自變量區(qū)間端點(diǎn)a,b步長(zhǎng)h后，計(jì)算機(jī)則給出常微

24、方程初值問題數(shù)值解。程序中變量說明：fx,y: 存放函數(shù)f(x,y)y0: 存放初值y0及數(shù)值解a：存放自變量區(qū)間左端點(diǎn)b: 存放自變量區(qū)間右端點(diǎn)n: 存放節(jié)點(diǎn)個(gè)數(shù)h: 存放節(jié)點(diǎn)步長(zhǎng)xk：存放節(jié)點(diǎn)xiy1: 存放數(shù)值解注：語句h=Input"步長(zhǎng)h="的步長(zhǎng)輸入應(yīng)該用帶小數(shù)點(diǎn)的數(shù)輸入，這樣可以加快計(jì)算速度，否則，由于Mathematica軟件本身的原因，可能會(huì)出現(xiàn)計(jì)算時(shí)間過長(zhǎng)等計(jì)算不出結(jié)果的情況。1.5 例題與實(shí)驗(yàn)例3. 用經(jīng)典的Runge-Kutta方法求初值問題的數(shù)值解。取步長(zhǎng)h=0.1，并與準(zhǔn)確解比較。解：執(zhí)行經(jīng)典的Runge-Kutta方法程序后，在輸入的窗口中按提

25、示分別輸入-2x*y2、1、0、1、0.1，每次輸入后用鼠標(biāo)點(diǎn)擊窗口的“OK”按扭，計(jì)算機(jī)在屏幕上給出如下數(shù)值解結(jié)果：y(0.1)=0.990099y(0.2)=0.961538y(0.3)=0.917431y(0.4)=0.862068y(0.5)=0.799999y(0.6)=0.735294y(0.7)=0.671141y(0.8)=0.609756y(0.9)=0.552487y(1.)=0.500001本題的準(zhǔn)確解為y(x)=(1+x2)-1，它在0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0的值為y(0.1)=0.990099, y(0.2)=0.9

26、61538, y(0.3)=0.917431, y(0.4)=0.862069,y(0.5)=0.8, y(0.6)=0.735294, y(0.7)=0.671141, y(0.8)=0.609756, y(0.9)= 0.552486, y(1.0)= 0.5,比較上面的計(jì)算結(jié)果，可見，計(jì)算結(jié)果比改進(jìn)的Euler法要好得多。例4. 用經(jīng)典的Runge-Kutta方法求初值問題的數(shù)值解。分別取步長(zhǎng)h=0.1，0.025,0.01計(jì)算并與準(zhǔn)確解在xi=0.1i,i=0,1,2,10各點(diǎn)比較。解：執(zhí)行經(jīng)典的Runge-Kutta方法程序后，在輸入的窗口中按提示分別輸入-50y+50x2+2x、1

27、/3、0、1、0.1，每次輸入后用鼠標(biāo)點(diǎn)擊窗口的“OK”按扭，計(jì)算機(jī)在屏幕上給出步長(zhǎng)h=0.1的如下數(shù)值解結(jié)果：y(0.1)=4.60549 y(0.2)=63.0625 y(0.3)=864.049 y(0.4)=11843.6 y(0.5)=162355.y(0.6)=2.22561´106 y(0.7)=3.05094´107 y(0.8)=4.18232´108y(0.9)=5.73327´109 y(1.)=7.85936´1010 再執(zhí)行一次經(jīng)典的Runge-Kutta方法程序，在輸入的窗口中按提示分別輸入-50y+50x2+2x、1/3、0、1、0.025，每次輸入后用鼠標(biāo)點(diǎn)擊窗口的“OK”按扭，計(jì)算機(jī)在屏幕上給出步長(zhǎng)h=0.025的部分?jǐn)?shù)值解結(jié)果：y(0.1)=0.0130149y(0.2)=0.0400633y(0.3)=0.090037y(0.4)=0.160037y(0.5)=0.250037y(0.6)=0.360037y(0.7)=0.490037y(0.8)=0.640037y(0.9)=0.810037y(1.)=1.00004又執(zhí)行一次經(jīng)