郭志萍

（浙江工業(yè)大學(xué)化材學(xué)院，浙江 杭州 310014；山西水利職業(yè)技術(shù)學(xué)院，山西 運(yùn)城 044004）

Runge-Kutta法在微分方程初值問題中的應(yīng)用

郭志萍

（浙江工業(yè)大學(xué)化材學(xué)院，浙江 杭州 310014；山西水利職業(yè)技術(shù)學(xué)院，山西 運(yùn)城 044004）

本文利用Taylor級數(shù)推導(dǎo)了二階Runge-Kutta算法，并得到了經(jīng)典的四階Runge-Kutta法。結(jié)合Matlab語言，在微分方程的初值問題中應(yīng)用Runge-Kutta法數(shù)值計(jì)算了一階微分方程和高階微分方程組。Runge-Kutta法能夠較高精度地?cái)?shù)值模擬數(shù)學(xué)建模所建立的微分方程。

Runge-Kutta法；Matlab；微分方程；數(shù)值計(jì)算

一、引言

微分方程是自然科學(xué)和社會科學(xué)研究中非常重要且十分實(shí)用的一種數(shù)學(xué)工具。Runge-Kutta法是一類具有相當(dāng)實(shí)用價值的方法，也是最常用的一種數(shù)值方法。在數(shù)值求解微分方程初值問題中應(yīng)用廣泛,張鵬利對兩空泡的共振頻率,共振振幅及空化噪聲聲壓利用Runge-Kutta法進(jìn)行了數(shù)值求解；Runge-Kutta法在振蕩電路求解中的應(yīng)用；淹沒磨料射流的運(yùn)動分析，等等。應(yīng)用Runge-Kutta法來數(shù)值求解，能夠較高精度地?cái)?shù)值模擬數(shù)學(xué)建模所建立的微分方程，其計(jì)算與實(shí)驗(yàn)結(jié)果吻合較好。這種方法能夠解決的問題越來越多，越來越復(fù)雜，在微分方程數(shù)值求解方法中扮演著非常重要的角色。

Matlab的基本單位是矩陣，它的表達(dá)式與數(shù)學(xué)、工程計(jì)算中常用的形式十分相似，極大地方便了人們學(xué)習(xí)和使用。利用Matlab提供各種數(shù)學(xué)工具，可以避免做繁瑣的數(shù)學(xué)推導(dǎo)和計(jì)算，方便地解決了很多數(shù)學(xué)問題。

本文從理論上對Runge-Kutta法進(jìn)行了推導(dǎo)，運(yùn)用Matlab語言，在具體的微分方程初值問題算例中詳細(xì)地給以求解。

二、初值問題數(shù)值解

除常系數(shù)線性微分方程可用特征根求解，少數(shù)特殊方程可用初等積分法求解外，大部分微分方程無顯式解，應(yīng)用中主要靠數(shù)值解?？紤]一階常微分方程組初值問題：

其中 y=(y1,y2,L,ym)T，f=(f1,f2,L,fm)T，y0=(y10,y20,L,ym0)T，這里表示轉(zhuǎn)置。所謂數(shù)值解，就是尋求y(t)在一系列離散節(jié)點(diǎn)t0＜t1＜L＜tn〈tf上的近似值 yk(k=0,1,L,n)。稱 hk=tk+1-tk為步長，通常取為常量h。

高階的常微分方程初值問題可以化為一階常微分方程組，已給一個n階方程：y(n)=f(t,y,y',L,y(n-1))T

設(shè)y1-y,y2=y',y1-y,L,yn=y(n-1),式(9)化為一階方程組：

三、Runge-kutta法推理過程

Runge-kutta方法是一種單步法，是用函數(shù)f在tk和tk-1之間值的有限差分近似代替高階導(dǎo)數(shù)。下面以二階為例來推導(dǎo)：Taylor級數(shù)

中多保留一項(xiàng)，可得到二階Taylor級數(shù)法

通過鏈?zhǔn)椒▽η髮?dǎo)，即

其中的下標(biāo)表示的是給定變量求偏導(dǎo)。

由此得到

用上式近似二階導(dǎo)數(shù)，得到Runge-Kutta法

其中 k1=f(tk,yk)，k2=f(tk+hk,yk+hkk1)

Runge-kutta方法優(yōu)點(diǎn)：計(jì)算tk+1時刻的值時，并不需要tk時刻以前的值。

最著名的是四階經(jīng)典的Runge-kutta格式

四、Runge-Kutta方法在Matlab語言中的應(yīng)用

下面以常用的ode45為例講解：

其調(diào)用格式為：[t,y]=ode45(odefun,tspan,y0)其中，odefun為用以表示的函數(shù)句，為標(biāo)量，y為標(biāo)量或向量；tspan如果是二維向量[t0,tf]，表示自變量初值和終值；如果是高維向量[t0,t1,…,tn]，則表示輸出節(jié)點(diǎn)列向量；y0表示節(jié)點(diǎn)列向量t=(t0,t1,L,tn)T；y表示數(shù)值解矩陣，每一列對應(yīng)y的一個分量；若無輸出參數(shù)，則作出圖形。

其中，options為計(jì)算參數(shù)，默認(rèn)可用空矩陣[]表示；p1,p2,…為附加傳遞參數(shù)，這時odefun的表示為f(t,y,pl,p2,L)。

運(yùn)行程序輸出和的關(guān)系。

Runge-kutta法有利于解決數(shù)學(xué)中的許多實(shí)際微分初值問題，方便、易操作，而且圖像處理清晰明了。Runge-kutta法是聯(lián)系數(shù)學(xué)與實(shí)際問題的橋梁，能夠較高精度地?cái)?shù)值模擬數(shù)學(xué)建模所建立的微分方程，其計(jì)算與實(shí)驗(yàn)結(jié)果吻合較好。這種方法能夠解決的問題越來越多，越來越復(fù)雜，在微分方程數(shù)值求解方法中扮演著非常重要的角色。

五、Runge-kutta法的意義

［1］戴嘉尊，邱建賢.微分方程數(shù)值解法［M］.南京:南京出版社，2002.

［2］何強(qiáng)，何英.MATLAB擴(kuò)展編程［M］.北京：清華大學(xué)出版社，2002.

［3］趙靜，但琦.數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn)（第2版）［M］.北京：高等教育出版社，2003.

［4］江世宏.MATLAB語言與數(shù)學(xué)實(shí)驗(yàn)［M］.北京：科學(xué)出版社，2006.

［5］董振海.精通MATLAB 7編程與數(shù)據(jù)庫應(yīng)用［M］.北京：電子工業(yè)出版社，2007.

O1

A

1673-0046（2011）01-0189-02