四階龍格-庫(kù)塔方法的程序設(shè)計(jì)與應(yīng)用

羅麗珍 吳慶軍

摘 要 本文通過(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)行比較。