《龍格庫(kù)塔方法及其matlab實(shí)現(xiàn)(共5頁(yè))》由會(huì)員分享,可在線閱讀,更多相關(guān)《龍格庫(kù)塔方法及其matlab實(shí)現(xiàn)(共5頁(yè))(5頁(yè)珍藏版)》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。
1、精選優(yōu)質(zhì)文檔-----傾情為你奉上
龍格-庫(kù)塔方法及其matlab實(shí)現(xiàn)
摘要:本文的目的數(shù)值求解微分方程精確解,通過(guò)龍格-庫(kù)塔法,加以利用matlab為工具達(dá)到求解目的。龍格-庫(kù)塔(Runge-Kutta)方法是一種在工程上應(yīng)用廣泛的高精度單步算法,用于數(shù)值求解微分方程。MatLab軟件是由美國(guó)Mathworks公司推出的用于數(shù)值計(jì)算和圖形處理的科學(xué)計(jì)算系統(tǒng)環(huán)境。MatLab是英文MATrix LABoratory(矩陣實(shí)驗(yàn)室)的縮寫。在MratLab環(huán)境下,用戶可以集成地進(jìn)行程序設(shè)計(jì)、數(shù)值計(jì)算、圖形繪制、輸入輸出、文件管理等各項(xiàng)操作。
關(guān)鍵詞:龍格-庫(kù)塔 matlab 微分方程
1.
2、 前言
1.1:知識(shí)背景
龍格-庫(kù)塔法(Runge-Kutta)是用于非線性常微分方程的解的重要的一類隱式或顯式迭代法。這些技術(shù)由數(shù)學(xué)家卡爾·龍格和馬丁·威爾海姆·庫(kù)塔于1900年左右發(fā)明。通常所說(shuō)的龍格庫(kù)塔方法是相對(duì)四階龍格庫(kù)塔而言的,成為經(jīng)典四階龍格庫(kù)塔法。該方法具有精度高,收斂,穩(wěn)定,計(jì)算過(guò)程中可以改變步長(zhǎng)不需要計(jì)算高階導(dǎo)數(shù)等優(yōu)點(diǎn),但是仍需計(jì)算在一些點(diǎn)上的值,比如四階龍格-庫(kù)塔法沒(méi)計(jì)算一步需要計(jì)算四步,在實(shí)際運(yùn)用中是有一定復(fù)雜性的。
Matlab是在20世紀(jì)七十年代后期的事:時(shí)任美國(guó)新墨西哥大學(xué)計(jì)算機(jī)科學(xué)系主任的Cleve Moler教授出于減輕學(xué)生編程負(fù)擔(dān)的動(dòng)機(jī),為學(xué)生設(shè)計(jì)了一組
3、調(diào)用LINPACK和EISPACK庫(kù)程序的“通俗易用”的接口,此即用FORTRAN編寫的萌芽狀態(tài)的MATLAB。
經(jīng)幾年的校際流傳,在Little的推動(dòng)下,由Little、Moler、Steve Bangert合作,于1984年成立了MathWorks公司,并把MATLAB正式推向市場(chǎng)。從這時(shí)起,MATLAB的內(nèi)核采用C語(yǔ)言編寫,而且除原有的數(shù)值計(jì)算能力外,還新增了數(shù)據(jù)圖視功能。
MATLAB以商品形式出現(xiàn)后,僅短短幾年,就以其良好的開放性和運(yùn)行的可靠性,使原先控制領(lǐng)域里的封閉式軟件包(如英國(guó)的UMIST,瑞典的LUND和SIMNON,德國(guó)的KEDDC)紛紛淘汰
4、,而改以MATLAB為平臺(tái)加以重建。在時(shí)間進(jìn)入20世紀(jì)九十年代的時(shí)候,MATLAB已經(jīng)成為國(guó)際控制界公認(rèn)的標(biāo)準(zhǔn)計(jì)算軟件。
到九十年代初期,在國(guó)際上30幾個(gè)數(shù)學(xué)類科技應(yīng)用軟件中,MATLAB在數(shù)值計(jì)算方面,而Mathematica和Maple則分居符號(hào)計(jì)算軟件的前兩名。Mathcad因其提供計(jì)算、圖形、文字處理的統(tǒng)一環(huán)境而深受中學(xué)生歡迎。
1.2研究的意義
精確求解數(shù)值微分方程,對(duì)龍格庫(kù)塔的深入了解與正確運(yùn)用,主要是在已知方程導(dǎo)數(shù)和初值信息,利用計(jì)算機(jī)仿真時(shí)應(yīng)用,省去求解微分方程的復(fù)雜過(guò)程。利用matlab強(qiáng)大的數(shù)值計(jì)算功能,省去認(rèn)為計(jì)算的過(guò)程,達(dá)到快速精確求解數(shù)值微分方程。在實(shí)際生活
5、中可以利用龍格庫(kù)塔方法和matlab的完美配合解決問(wèn)題。
1.3研究的方法
對(duì)實(shí)例的研究對(duì)比,實(shí)現(xiàn)精度的要求,龍格庫(kù)塔是并不是一個(gè)固定的公式,所以只是對(duì)典型進(jìn)行分析
2. 龍格-庫(kù)塔方法
2.1龍格-庫(kù)塔公式
在一階精度的的拉格朗日中值定理有:
對(duì)于函數(shù)y=f(x,y)
y'=f(x,y)
y(n+1)=y(n)+h*K1
K1=f(xn, yn)
這就是一階龍格-庫(kù)塔方法
形如 y(n+1)=y(n)+h*i=1rciki
k1 =f(xn,yn)
ki=f(xn+hai ,yn+h*j=1i-1bijki)
6、 i=2…r
故二階龍格-庫(kù)塔公式
y(n+1)=y(n)+h(c1k1+c2k2)
k1= f(xn,yn) (2)
k2= f(xn+ha2 ,yn+ha2 k1)
將y(x)在xn處展成冪級(jí)數(shù)
y(xn+1)=y(xn)+hy'(xn)+h22y'’ (xn)+o(h3)
y'(x)= f(x,y(x))
y'’x= fx‘(x,y(x))+ fy‘(x,y(x))·f(x,y(x))
y(xn+1)=y(xn)+hf+
7、h22(fx‘+fy‘f)+ o(h3) (3)
將(2)式中的k2在(xn,yn)點(diǎn)展成冪級(jí)數(shù)
k2= f(xn+ha2 ,yn+ha2 k1)
=f+ha2fx‘+ ha2fy‘f+ o(h2)
將k1,k2代入(2)式,得
yn+1=yn +h(c1+c2)f
+ha2c2(fx‘+fy‘f)+ o(h3) (4)
對(duì)比(3)(4),當(dāng)y(xn)= yn時(shí)
只有c1+c2=1,a2c2=12 (5)
形如(2)存在常數(shù)滿足(5)式,局部
8、截?cái)嗾`差為o(h3)的求解方法稱為二階龍格-庫(kù)塔法。
滿足(5)式,若取c1=12,則得到c2=12,a2=1,則公式則恰為預(yù)估-校正法公式
若取c1=0,則c2=1,a2=12,
yn+1=yn+hk2
k1= f(xn,yn) (6)
k2= f(xn+h2,yn+h2k)
9、 n=0,1…N-1
由(5)式,可知龍格-庫(kù)塔法不是唯的
三階龍格-庫(kù)塔法
yn+1=yn+h(c1k1+c2k2+ c3k3)
k1= f(xn,yn)
k2= f(xn+ha2,yn+ha2k1) (7)
k3= f(xn+ha3,yn+hb31k1+hb32k2)
若c1,c2, c3,a2,a3,b31, b32且滿足b31+ b
10、32=a3,,并使得局部截?cái)嗾`差為o(h4)。類似二階龍格-庫(kù)塔法推導(dǎo)的
c1+c2+ c3=1
a2c2+a3c3=12
a2b32c3=16 (8)
a22c2+a32c3=13
形如(7),常數(shù)滿足(8),局部截?cái)嗾`差為o(h4)的求解方法稱為三階龍格-庫(kù)塔法
在(8)式
11、中若取c1=16,c3=16,則得c2=23,a2=12,a3=1,b31=-1,b32=2
代入(7)中得三階龍格-庫(kù)塔法公式
yn+1=yn+h6(k1+4k2+ k3)
k1= f(xn,yn)
k2= f(xn+h2,yn+h2k1) (9)
k3= f(xn+ha3,yn-hk1+2hk2)
四階龍格庫(kù)塔法的推導(dǎo)類似于三階龍格-庫(kù)塔法,但相對(duì)復(fù)雜
12、這里不再進(jìn)行推導(dǎo),公式如下
yn+1=yn+h6(k1+2k2+2 k3+k4)
k1= f(xn,yn)
k2= f(xn+h2,yn+h2k1) (10)
k3= f(xn+h2,yn+h2k2)
k4= f(xn+h,yn+h k3)
n=0,1…N-1
這就是標(biāo)準(zhǔn)
13、四階龍格庫(kù)塔公式
2.1 對(duì)實(shí)例的研究
利用龍格-庫(kù)塔法求解方程
y'=8-3yy0=2的數(shù)值,其中h=0.2,計(jì)算y(0.4)的近似值。至少保留四位小數(shù)。
解:f?(x,y)=8-3y,利用四階龍格-庫(kù)塔公式有
yn+1=yn+h6(k1+2k2+2 k3+k4)
k1= f(xn,yn)=8-3yn
k2= f(xn+h2,yn+h2k1)=5.6-2.1yn
14、
k3= f(xn+h2,yn+h2k2)=6.32-2.37yn
k4= f(xn+h,yn+h k3)=4.208-1.578yn
n=0,1…N-1
yn+1=yn1.2016+0.5494yn
當(dāng)x0=0,y0=2,
y(0.2)≈y1=1.2016+0.5494y0=1.2016+0.5494×2=2.3004
y(0.4)≈y2=1.2016+0.5494y1=1.2016+0.5494×2.3004=2.4654
專心---專注---專業(yè)