羅麗珍 吳慶軍
摘 要 本文通過(guò)介紹四階龍格-庫(kù)塔方法,通過(guò)預(yù)報(bào)斜率和泰勒展開(kāi)式推導(dǎo)出龍格—庫(kù)塔格式。了解它的基本思想與算法步驟、MATLAB語(yǔ)言編寫(xiě)的程序。列舉一些例子,運(yùn)用四階龍格-庫(kù)塔方法的MATLAB程序在軟件中運(yùn)行,求解出常微分方程的數(shù)值解,同時(shí)將求解出的數(shù)值解與精確解進(jìn)行比較。
關(guān)鍵詞 龍格-庫(kù)塔方法 常微分方程 數(shù)值解
中圖分類(lèi)號(hào):TP337文獻(xiàn)標(biāo)識(shí)碼:A
0引言
從17世紀(jì)以來(lái)國(guó)內(nèi)外數(shù)學(xué)家對(duì)常微分方程的研究取得了很多的成果.歐拉在研究中指出常微分方程存在唯一解和無(wú)數(shù)解,他用近似值求解微分方程,發(fā)現(xiàn)用積分因子求解微積分方程的特殊算法。拉格朗日建立了一階微分方程理論,他將參數(shù)變法應(yīng)用到四階非齊次方程的求解。
我們生活中許多問(wèn)題的解決都運(yùn)用到常微分方程,常微分方程的數(shù)值解法中經(jīng)常使用的方法是四階龍格-庫(kù)塔方法。各個(gè)領(lǐng)域和工程問(wèn)題中的原理和演變規(guī)律都是用常微分方程來(lái)描述的,如在物理方面的電路中電流變化的規(guī)律、航天航空方面衛(wèi)星運(yùn)轉(zhuǎn)問(wèn)題、經(jīng)濟(jì)方面物品供給以及需求與物價(jià)的之間的關(guān)系、軍事方面研究深水炸彈在水下的運(yùn)動(dòng)等。對(duì)這些事物、現(xiàn)象變化規(guī)律的描述、認(rèn)知和分析,需要運(yùn)用常微分方程來(lái)解決。人們使用常微分方程數(shù)值解法的四階龍格-庫(kù)塔方法去研究這些問(wèn)題很實(shí)用,而且具有很重要的應(yīng)用價(jià)值。
目前,常微分方程在解決我們生活中的問(wèn)題很實(shí)用,許多問(wèn)題都運(yùn)用常微分方程來(lái)求解。中國(guó)科學(xué)技術(shù)大學(xué)學(xué)者倪興在常微分方程的研究中寫(xiě)了關(guān)于歐拉法、方法等幾種方法,他運(yùn)用常微分計(jì)算衛(wèi)星運(yùn)動(dòng)的初軌,把方法運(yùn)用到衛(wèi)星軌道改進(jìn)的例子中;揚(yáng)州大學(xué)學(xué)者馮建強(qiáng)和孫詩(shī)一研究四階方法的推導(dǎo),他寫(xiě)出了如何推導(dǎo)的過(guò)程。在高校數(shù)值分析、數(shù)值計(jì)算方法與實(shí)驗(yàn)等教材中,許多作者都出版關(guān)于常微分方程初值問(wèn)題數(shù)值解法的教材書(shū),歐拉方法、改進(jìn)歐拉法和方法等,同時(shí)在教材書(shū)中寫(xiě)入各種實(shí)際問(wèn)題的例子,運(yùn)用這些方法去解決常微分方程的初值問(wèn)題。本文主要介紹常微分方程數(shù)值解法的四階方法,使用四階方法求解常微分方程。介紹四階方法的基本思想、算法步驟和MATLAB程序,同時(shí)使用四階方法在案例中解決問(wèn)題。
1四階方法
1.1四階方法的思想與算法
對(duì)于一階常微分方程初值問(wèn)題,根據(jù)拉格朗日微分中值定理可知存在使
記為在區(qū)間上的平均斜率。
在給定的區(qū)間里預(yù)測(cè)個(gè)點(diǎn)處的斜率,再把它們的加權(quán)平均作為平均斜率的近似值。這里取在區(qū)間上若干個(gè)點(diǎn)的斜率值或者預(yù)報(bào)斜率值的加權(quán)平均值,作為平均斜率的近似值,令區(qū)間上若干個(gè)點(diǎn)的斜率值或預(yù)報(bào)斜率值為,以及權(quán)系數(shù)為,使差分格式
當(dāng)上式為階時(shí),把它叫做階龍格—庫(kù)塔格式.通過(guò)上式能夠推導(dǎo)得精確度比較高的公式來(lái)求解常微分方程。
當(dāng)時(shí),龍格-庫(kù)塔的格式為
在上式中,在點(diǎn)的斜率是在點(diǎn)的預(yù)報(bào)斜率為是,參數(shù)都是待定參數(shù)。當(dāng)待定參數(shù)使上式是二階格式時(shí),稱(chēng)上式是二階龍格-庫(kù)塔格式。
函數(shù)在點(diǎn)的泰勒展開(kāi)式為
令是準(zhǔn)確的,有,根據(jù)二元函數(shù)泰勒展開(kāi)有
要使上面的截距誤差為,把與泰勒展開(kāi)式比較可得。通過(guò)以上的推導(dǎo)可以得出待定參數(shù)的值。
由可知二階龍格-庫(kù)塔格式是一個(gè)系列的差分格式。當(dāng)取時(shí)得
這是改進(jìn)的歐拉公式。取,可得
其中
這是常用的二階龍格-庫(kù)塔方法,稱(chēng)其為中心格式。
三階龍格-庫(kù)塔格式:
三階龍格-庫(kù)塔格式是比二階格式高一階的格式,它是在二階龍格-庫(kù)塔格式的基礎(chǔ)上進(jìn)一步構(gòu)造得到的格式。當(dāng)時(shí),龍格-庫(kù)塔的格式如下
在上式里,在點(diǎn)的斜率是在點(diǎn)和點(diǎn)的預(yù)報(bào)斜率分別是、,當(dāng)待定參數(shù)和使上式為三階格式時(shí),則稱(chēng)其為三階龍格—庫(kù)塔格式。
類(lèi)似二階龍格-庫(kù)塔格式的推導(dǎo),得到三階龍格-庫(kù)塔格式為:
稱(chēng)其為三階龍格-庫(kù)塔格式。
四階龍格-庫(kù)塔格式:
同樣是根據(jù)上述的推導(dǎo),得到四階龍格-庫(kù)塔格式為:
稱(chēng)上式為四階龍格-庫(kù)塔格式。在解決許多問(wèn)題中常常使用的龍格-庫(kù)塔格式是四階龍格-庫(kù)塔格式。
四階龍格-庫(kù)塔方法的算法步驟如下:
第一步輸入?yún)^(qū)間等分?jǐn)?shù),初值;
第二步輸出在的個(gè)點(diǎn)處的近似值;
第三步令;
第四步計(jì)算
令,輸出;
第五步如果,令,轉(zhuǎn)向第四步;否則停機(jī)。
1.2四階方法的MATLAB程序
四階龍格-庫(kù)塔法MATLAB函數(shù)文件nark4.m如下:
function [x,y]=nark4(dyfun,xspan,y0,h)
x=xspan(1):h:xspan(2);
y(1)=y0;
for n=1:length(x)-1
k1=dyfun(x(n),y(n));
K2=dyfun(x(n)+h/2,y(n)+h/2*k1);
K3=dyfun(x(n)+h/2,y(n)+h/2*k2);
k4=dyfun(x(n+1),y(n)+h*k3);
y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6
end
x=x';y=y';
1.3四階方法的應(yīng)用
例2.1:用四階方法求解初值問(wèn)題的數(shù)值解,取步長(zhǎng)。其精確解為,并將解得的數(shù)值解與精確解進(jìn)行比較。