龍格庫(kù)塔方法及其matlab實(shí)現(xiàn)

《龍格庫(kù)塔方法及其matlab實(shí)現(xiàn)》由會(huì)員分享，可在線閱讀，更多相關(guān)《龍格庫(kù)塔方法及其matlab實(shí)現(xiàn)（7頁珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

1、 龍格-庫(kù)塔方法及其matlab實(shí)現(xiàn) 摘要：本文的目的數(shù)值求解微分方程精確解，通過龍格-庫(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. 前言 1.1：知識(shí)背景 龍格

2、－庫(kù)塔法（Runge-Kutta）是用于非線性常微分方程的解的重要的一類隱式或顯式迭代法。這些技術(shù)由數(shù)學(xué)家卡爾龍格和馬丁威爾海姆庫(kù)塔于1900年左右發(fā)明。通常所說的龍格庫(kù)塔方法是相對(duì)四階龍格庫(kù)塔而言的，成為經(jīng)典四階龍格庫(kù)塔法。該方法具有精度高，收斂，穩(wěn)定，計(jì)算過程中可以改變步長(zhǎng)不需要計(jì)算高階導(dǎo)數(shù)等優(yōu)點(diǎn)，但是仍需計(jì)算在一些點(diǎn)上的值，比如四階龍格-庫(kù)塔法沒計(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ì)了一組調(diào)用LINPACK和EISPACK庫(kù)程序

3、的“通俗易用”的接口，此即用FORTRAN編寫的萌芽狀態(tài)的MATLAB。 經(jīng)幾年的校際流傳，在Little的推動(dòng)下，由Little、Moler、Steve Bangert合作，于1984年成立了MathWorks公司，并把MATLAB正式推向市場(chǎng)。從這時(shí)起，MATLAB的內(nèi)核采用C語言編寫，而且除原有的數(shù)值計(jì)算能力外，還新增了數(shù)據(jù)圖視功能。 MATLAB以商品形式出現(xiàn)后，僅短短幾年，就以其良好的開放性和運(yùn)行的可靠性，使原先控制領(lǐng)域里的封閉式軟件包（如英國(guó)的UMIST，瑞典的LUND和SIMNON，德國(guó)的KEDDC）紛紛淘汰，而改以MATLAB為平臺(tái)加以重建。在時(shí)

4、間進(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ì)算方面獨(dú)占鰲頭，而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ù)雜過程。利用matlab強(qiáng)大的數(shù)值計(jì)算功能，省去認(rèn)為計(jì)算的過程，達(dá)到快速精確求解數(shù)值微分方程。在實(shí)際生活中可以利用龍格庫(kù)塔方法和matl

5、ab的完美配合解決問題。 1.3研究的方法 對(duì)實(shí)例的研究對(duì)比，實(shí)現(xiàn)精度的要求，龍格庫(kù)塔是并不是一個(gè)固定的公式，所以只是對(duì)典型進(jìn)行分析 請(qǐng)預(yù)覽后下載！ 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+h22(fx‘+fy‘

7、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） 請(qǐng)預(yù)覽后下載！ 對(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ù)塔法公式 請(qǐng)預(yù)覽后下載！ 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)類似于三階龍

12、格-庫(kù)塔法，但相對(duì)復(fù)雜這里不再進(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

13、…N-1 這就是標(biāo)準(zhǔn)四階龍格庫(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， 請(qǐng)預(yù)覽后下載！ y(0.2)≈y1＝1.2016＋0.5494y0＝1.2016＋0.54942＝2.3004 y(0.4)≈y2＝1.2016＋0.5494y1＝1.2016＋0.54942.3004＝2.4654 （注：可編輯下載，若有不當(dāng)之處，請(qǐng)指正，謝謝!） 請(qǐng)預(yù)覽后下載！