羅麗珍 吳慶軍

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

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

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

0引言

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

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

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

1四階方法

1.1四階方法的思想與算法

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

記為在區(qū)間上的平均斜率。

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

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

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

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

函數(shù)在點的泰勒展開式為

令是準確的，有，根據(jù)二元函數(shù)泰勒展開有

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

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

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

其中

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

三階龍格-庫塔格式：

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

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

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

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

四階龍格-庫塔格式：

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

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

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

第一步輸入區(qū)間等分數(shù)，初值;

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

第三步令;

第四步計算

令，輸出;

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

1.2四階方法的MATLAB程序

四階龍格-庫塔法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四階方法的應用

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