四階龍格-庫塔方法的程序設計與應用

羅麗珍 吳慶軍

摘 要 本文通過介紹四階龍格-庫塔方法，通過預報斜率和泰勒展開式推導出龍格—庫塔格式。了解它的基本思想與算法步驟、MATLAB語言編寫的程序。列舉一些例子，運用四階龍格-庫塔方法的MATLAB程序在軟件中運行，求解出常微分方程的數值解，同時將求解出的數值解與精確解進行比較。

關鍵詞 龍格-庫塔方法 常微分方程 數值解

中圖分類號：TP337文獻標識碼：A

0引言

從17世紀以來國內外數學家對常微分方程的研究取得了很多的成果.歐拉在研究中指出常微分方程存在唯一解和無數解，他用近似值求解微分方程，發現用積分因子求解微積分方程的特殊算法。拉格朗日建立了一階微分方程理論，他將參數變法應用到四階非齊次方程的求解。

我們生活中許多問題的解決都運用到常微分方程，常微分方程的數值解法中經常使用的方法是四階龍格-庫塔方法。各個領域和工程問題中的原理和演變規律都是用常微分方程來描述的，如在物理方面的電路中電流變化的規律、航天航空方面衛星運轉問題、經濟方面物品供給以及需求與物價的之間的關系、軍事方面研究深水炸彈在水下的運動等。對這些事物、現象變化規律的描述、認知和分析，需要運用常微分方程來解決。人們使用常微分方程數值解法的四階龍格-庫塔方法去研究這些問題很實用，而且具有很重要的應用價值。

目前，常微分方程在解決我們生活中的問題很實用，許多問題都運用常微分方程來求解。中國科學技術大學學者倪興在常微分方程的研究中寫了關于歐拉法、方法等幾種方法，他運用常微分計算衛星運動的初軌，把方法運用到衛星軌道改進的例子中;揚州大學學者馮建強和孫詩一研究四階方法的推導，他寫出了如何推導的過程。在高校數值分析、數值計算方法與實驗等教材中，許多作者都出版關于常微分方程初值問題數值解法的教材書，歐拉方法、改進歐拉法和方法等，同時在教材書中寫入各種實際問題的例子，運用這些方法去解決常微分方程的初值問題。本文主要介紹常微分方程數值解法的四階方法，使用四階方法求解常微分方程。介紹四階方法的基本思想、算法步驟和MATLAB程序，同時使用四階方法在案例中解決問題。

1四階方法

1.1四階方法的思想與算法

對于一階常微分方程初值問題，根據拉格朗日微分中值定理可知存在使

記為在區間上的平均斜率。

在給定的區間里預測個點處的斜率，再把它們的加權平均作為平均斜率的近似值。這里取在區間上若干個點的斜率值或者預報斜率值的加權平均值，作為平均斜率的近似值，令區間上若干個點的斜率值或預報斜率值為，以及權系數為，使差分格式

當上式為階時，把它叫做階龍格—庫塔格式.通過上式能夠推導得精確度比較高的公式來求解常微分方程。

當時，龍格-庫塔的格式為

在上式中，在點的斜率是在點的預報斜率為是，參數都是待定參數。當待定參數使上式是二階格式時，稱上式是二階龍格-庫塔格式。

函數在點的泰勒展開式為

令是準確的，有，根據二元函數泰勒展開有

要使上面的截距誤差為，把與泰勒展開式比較可得。通過以上的推導可以得出待定參數的值。

由可知二階龍格-庫塔格式是一個系列的差分格式。當取時得

這是改進的歐拉公式。取，可得

其中

這是常用的二階龍格-庫塔方法，稱其為中心格式。

三階龍格-庫塔格式：

三階龍格-庫塔格式是比二階格式高一階的格式，它是在二階龍格-庫塔格式的基礎上進一步構造得到的格式。當時，龍格-庫塔的格式如下

在上式里，在點的斜率是在點和點的預報斜率分別是、，當待定參數和使上式為三階格式時，則稱其為三階龍格—庫塔格式。

類似二階龍格-庫塔格式的推導，得到三階龍格-庫塔格式為：

稱其為三階龍格-庫塔格式。

四階龍格-庫塔格式：

同樣是根據上述的推導，得到四階龍格-庫塔格式為：

稱上式為四階龍格-庫塔格式。在解決許多問題中常常使用的龍格-庫塔格式是四階龍格-庫塔格式。

四階龍格-庫塔方法的算法步驟如下：

第一步輸入區間等分數，初值;

第二步輸出在的個點處的近似值;

第三步令;

第四步計算

令，輸出;

第五步如果，令，轉向第四步;否則停機。

1.2四階方法的MATLAB程序

四階龍格-庫塔法MATLAB函數文件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四階方法的應用

例2.1：用四階方法求解初值問題的數值解，取步長。其精確解為，并將解得的數值解與精確解進行比較。