一種數(shù)值積分算法的改進(jìn)研究

吳佳惠

一種數(shù)值積分算法的改進(jìn)研究

吳佳惠

(新疆大學(xué) 電氣工程學(xué)院，新疆烏魯木齊 830047)

本文分析了Matlab軟件常用的常系數(shù)微分方程的數(shù)值積分算法和含有間斷特性的微分方程求解方法，重點(diǎn)研究了ODE45算法，并針對(duì)其求解不連續(xù)問(wèn)題有可能出現(xiàn)誤差，采用條件函數(shù)零點(diǎn)搜索法進(jìn)行了改進(jìn)，通過(guò)用ode45，ode23，ode113，ode15s，ode23s，ode23t和ode23tb求解器求解算例和本文用改進(jìn)的ODE算法求解算例比較，驗(yàn)證了本文改進(jìn)算法的有效性和適用性。

數(shù)值積分；微分方程；ODE45

0 引言

物理學(xué)、機(jī)械學(xué)、電學(xué)和熱力學(xué)的基本定律通常都是建立在經(jīng)驗(yàn)性的觀測(cè)資料之上的，這些觀測(cè)資料可以解釋物理性質(zhì)和系統(tǒng)狀態(tài)變化，定律一般不直接描述物理系統(tǒng)的狀態(tài)，而是通過(guò)時(shí)間和空間的變化來(lái)表達(dá)。這些定律定義了變化的機(jī)制，一旦與能量、質(zhì)量或動(dòng)量等連續(xù)性定律結(jié)合，就得到了常微分方程[1]。數(shù)值積分是求解這類方程非常重要而常用的數(shù)值計(jì)算方，也是在實(shí)際中得到廣泛應(yīng)用的算法[9]。

Matlab軟件中一般含有ode45(Ordinary Difference Equation, ODE)，ode23，ode113，ode15s，ode23s，ode23t和ode23tb求解器，ode45是四/五階龍格-庫(kù)塔法，適用于大多數(shù)連續(xù)或離散系統(tǒng)，但不適用于剛性（stiff）系統(tǒng)；ode23是二/三階龍格-庫(kù)塔法；ode113是一種階數(shù)可變的解法器；ode113是一種多步解法器；ode113是變階Adams法；ode15s：是一種基于數(shù)字微分公式的解法器；ode23s是一種單步解法器；ode23s是定階單步法；ode23t是梯形規(guī)則的一種自由插值實(shí)現(xiàn)；ode23tb是具有兩個(gè)階段的隱式龍格-庫(kù)塔公式。在這些公式中ode45應(yīng)用最廣泛[7][10]，但是，函數(shù)的連續(xù)性對(duì)保證計(jì)算精度是十分重要的，單步法要求在一個(gè)步長(zhǎng)內(nèi)右端函數(shù)連續(xù)，K階多步法則要求在K步內(nèi)右端函數(shù)連續(xù)，但是在實(shí)際問(wèn)題中經(jīng)常遇到不連續(xù)的情況，條件函數(shù)零點(diǎn)搜索法解決具有一般間斷特性時(shí)的求解方法[2]。若能事先確定滿足條件函數(shù)的點(diǎn)，然后就可以恰當(dāng)?shù)剡x擇步長(zhǎng)。根據(jù)條件函數(shù)的特點(diǎn)，若條件函數(shù)僅與時(shí)間有關(guān)，事先確定滿足條件函數(shù)的時(shí)間點(diǎn)值，計(jì)算推進(jìn)時(shí)，每一步均與間斷點(diǎn)時(shí)間值進(jìn)行比較，使當(dāng)前步不跨過(guò)間斷點(diǎn)[2]。

國(guó)內(nèi)很多學(xué)者針對(duì)各種情況研究了數(shù)值積分算法，取得了大量可以借鑒的成果[3-6][8]。本文針對(duì)matlab軟件中常用的ODE45算法進(jìn)行改進(jìn)研究，采用件函數(shù)零點(diǎn)搜索法與ODE45算法結(jié)合，可以有效求解含有間斷特性的常微分方程問(wèn)題。

1 RKF45算法改進(jìn)

RKF（Runge–Kutta–Fehlberg）45即4階5級(jí)龍格-庫(kù)塔-費(fèi)爾別格法，也matlab軟件中ODE45算法，具體表達(dá)式如下：

式中的各系數(shù)見表1。

主要軟件代碼為：

while (tout(i)＜=tend)

tout(i+1)=tout(i)+h;

k1=f(yout(i));

k2=f(yout(i)+h*k1*1/4);

k3=f(yout(i)+h*(k1*3/32+k2*9/32));

k4=f(yout(i)+h*(k1*1932/2197–k2*7200/21 97+k3*7296/2197));

k5=f(yout(i)+h*(k1*439/216–k2*8+k3*368 0/513–k4*845/4104));

k6=f(yout(i)+h*(–k1*8/27+k2*2–k3*3544/2 565+k4*1859/4104–k5*11/40));

yout(i+1)=yout(i)+h*(k1*25/216+k2*1408/2565 +k4*2197/4104–k5*1/5);

y_hat=yout(i)+h*(k1*16/135+k3*6656/12825+k 4*28561/56430–k5*9/50+k6*2/55);

表1 Fehlberg公式中的各系數(shù)Tab.1 The Coefficients in Fehlberg Formula

由文獻(xiàn)[2]可知，當(dāng)用經(jīng)典的數(shù)值積分法（包括單步法和多步法）進(jìn)行數(shù)值求解時(shí)，右端函數(shù)的連續(xù)性對(duì)保證仿真精度是十分重要的。單步法要求在一個(gè)步長(zhǎng)內(nèi)右端函數(shù)連續(xù)，K階多步法則要求在K步內(nèi)右端函數(shù)連續(xù)，為了使仿真結(jié)果達(dá)到一定精度，人們往往采用變步長(zhǎng)方法。然而，采用變步長(zhǎng)方法也會(huì)遇到困難。

由文獻(xiàn)[2]可知，具有一般間斷特性時(shí)的微分方程數(shù)值求解問(wèn)題。若系統(tǒng)微分方程y˙=f(y,t)中右端函數(shù)f(y,t)不連續(xù)，不失一般性，可用以下形式來(lái)描述

即右端函數(shù)的表達(dá)式依賴于條件函數(shù)的取值。對(duì)上式所描述的系統(tǒng)，若能事先確定f(y,t)滿足條件函數(shù)的點(diǎn)，即滿足φi(y,t)=0的t*，然后就可以恰當(dāng)?shù)剡x擇步長(zhǎng)。根據(jù)條件函數(shù)的特點(diǎn)，條件函數(shù)若僅與時(shí)間有關(guān)，即φi(y,t)=φ(t)=0；事先確定滿足條件函數(shù)的時(shí)間點(diǎn)值，數(shù)值求解推進(jìn)時(shí)，每一步均與間斷點(diǎn)時(shí)間值進(jìn)行比較，使當(dāng)前步不跨過(guò)間斷點(diǎn)，這就是條件函數(shù)零點(diǎn)搜索法。

在經(jīng)典RKF45算法中加入條件函數(shù)零點(diǎn)搜索法，可以使其更好地求解某些含有間斷點(diǎn)的微分方程求解問(wèn)題，其代碼為：

if ((yout(i+1)＞=yhigh)&(yout(i+1)＞yout(i)))

if (yout(i+1) – yhigh ＞ tol)

h = h *(yhigh–yout(i))/(yout(i+1) –yout(i));

continue;

end

siphon=1;

elseif(yout(i+1)＜=ylow&(yout(i+1)＜yout(i)))

if (ylow–yout(i+1) ＞ tol)

h = h *(yout(i)–ylow)/(yout(i) –yout(i+1));

continue;

end

siphon=0;

end

2 案例驗(yàn)證

本文案例來(lái)自文獻(xiàn)[1]，有一座間歇式噴泉，如圖1。水以固定的流率Qin注入圓柱形容器，水位達(dá)到y(tǒng)high時(shí)容器被注滿。這時(shí)候，水通過(guò)圓形排出管被虹吸出容器，在管子的盡頭形成噴泉。噴泉一直噴射到水位線降至ylow，此時(shí)虹吸管中裝滿了空氣，于是噴泉停止噴射。然后重復(fù)這個(gè)過(guò)程，當(dāng)水位線到達(dá)yhigh，即容器被注滿時(shí)，噴泉又開始噴射。當(dāng)噴泉噴射時(shí)，根據(jù)托里切利定律（Torricelli’s law）,流出物Qout可以由下面的公式計(jì)算，

圖1 間歇式噴泉Fig.1 Intermittent Fountain

忽略管中水的體積，計(jì)算100 s的時(shí)間內(nèi)容器水位線與時(shí)間的函數(shù)關(guān)系并繪制圖形，假設(shè)空容器的初始條件為y(0)=0，計(jì)算時(shí)使用下面的參數(shù)RT=0.05 m，r=0.007 m，ylow=0.025 m，yhigh=0.1，C=0.6，g=9.81 m/s2，Qin=50x10–6 m3/s。通過(guò)ode23，ode113，ode15 s，ode23 s，ode23 t和ode23 tb來(lái)求解結(jié)果如下圖所示。圖的結(jié)果明顯是不正確的，除了最初的注滿周期之外，水位線看起來(lái)在再次達(dá)到y(tǒng)high以前就開始下降。類似地，在排水的時(shí)候，水位線還沒(méi)有降至ylow，虹吸管就關(guān)上了。

[1]可知，當(dāng)噴泉噴射時(shí)，容器體積V(m3)的變化率由流入減去流出大的簡(jiǎn)單平衡確定

當(dāng)噴泉停止噴射時(shí)，分子的第二項(xiàng)變成0。為此，我們可以在方程中加入一個(gè)新的無(wú)量綱的變量siphon，當(dāng)噴泉停止時(shí)siphon等于0，當(dāng)噴泉工作時(shí)siphon等于1。

siphon可以被看成是控制噴泉停止和工作的開關(guān)。這種只具有兩種狀態(tài)的量被稱為布爾變量或者邏輯變量，其中，0相當(dāng)于假，1相當(dāng)于真；將siphon與應(yīng)變量y聯(lián)系起來(lái)。首先，當(dāng)水位線低于ylow時(shí)令siphon等于0。相應(yīng)地，當(dāng)水位線高于yhigh時(shí)令siphon等于1，下面的M文件函數(shù)在計(jì)算導(dǎo)數(shù)時(shí)考慮了這個(gè)邏輯關(guān)系。

Function dy=plinyode(t,y)

Global siphon

Rt=0.05; r=0.007; yhi=0.1; ylo=0.025

C=0.6; g=9.81; Qin=0.00005;

if y(1)＜=ylo; siphon=0;

elseif y(1)＞=yhi

siphon=1;

end

Qout=siphon*Csqrt(2*g*y(1))*pi*r^2;

Dy=(Qin-Qout)/(pi*Rt^2);

由于siphon的取值必須在函數(shù)調(diào)用之間被保持，所以它被申明為一個(gè)全局變量。雖然，我們不鼓勵(lì)使用全局變量（特別是對(duì)于大型的程序），但是，它在當(dāng)前的問(wèn)題中很有用。下面的腳本用內(nèi)置ode45函數(shù)積分Plinyode，并繪制解的圖形。

global siphon

siphon=0;

tspan=[0 100]; y0=0;

[tp,yp]=ode113(@plinyode, tspan, y0);

plot(tp, yp);

xlabel(‘time, (s)’) ;

ylabel(‘water level in tank, (m)’)

[tp,yp]=ode15s(@plinyode, tspan, y0);

plot(tp, yp);

xlabel(‘time, (s)’) ;

ylabel(‘water level in tank, (m)’)

[tp, yp]=ode23(@plinyode, tspan, y0);

plot(tp, yp);

xlabel(‘time, (s)’);

ylabel(‘water level in tank, (m)’)

[tp, yp]=ode23s(@plinyode, tspan, y0);

plot(tp, yp);

xlabel(‘time, (s)’);

ylabel(‘water level in tank, (m)’);

[tp, yp]=ode23t(@plinyode, tspan, y0);

plot(tp, yp) ;

xlabel(‘time,(s)’);

ylabel(‘water level in tank, (m)’);

[tp, yp]=ode23tb(@plinyode, tspan, y0);

plot(tp, yp);

xlabel(‘time, (s)’);

ylabel(‘water level in tank, (m)’);

圖的結(jié)果明顯是不正確的。除了最初的注滿周期之外，水位線看起來(lái)在再次達(dá)到y(tǒng)high以前就開始下降。類似地，在排水的時(shí)候，水位線還沒(méi)有降至ylow，虹吸管就關(guān)上了。

圖2 Matlab軟件中ODE求解算法求解的結(jié)果Fig.2 The results of ode algorithm in MATLAB software

由文獻(xiàn)[1]可知，是因?yàn)閛de在虹吸管開關(guān)時(shí)并不連續(xù)，例如，當(dāng)容器注水時(shí)，導(dǎo)數(shù)僅僅依賴于固定的流入物和取值恒為某個(gè)參數(shù)。然而，一旦水位線到達(dá)yhigh，流出開關(guān)導(dǎo)通，導(dǎo)數(shù)迅速地變?yōu)榱硪粋€(gè)值，雖然Matlab所用的自適應(yīng)步長(zhǎng)程序?qū)τ谠S多問(wèn)題微分方程求解問(wèn)題十分有效，但是，它們?cè)谔幚磉@類間斷問(wèn)題時(shí)往往會(huì)失效。我們用ODE45和改進(jìn)的ODE45算法求解的結(jié)果如下，可以看出ODE45同樣出現(xiàn)這種困難，是因?yàn)镺DE在虹吸管開關(guān)時(shí)并不連續(xù)，但是采用改進(jìn)的ODE45算法，就沒(méi)有出現(xiàn)問(wèn)題。因?yàn)椋ㄟ^(guò)改進(jìn)經(jīng)典ODE45算法，加入了條件函數(shù)零點(diǎn)搜索法，所以，使得ODE45算法可以求解某些含有間斷特性的微分方程求解問(wèn)題。

圖3 ODE45和改進(jìn)的ODE45求解結(jié)果Fig.3 Ode45 and improved ode45 solution results

3 結(jié)論

本文分析了Matlab軟件中一般含有ode45(Ordinary Difference Equation, ODE)，ode23，ode113，ode15s，ode23s，ode23t和ode23tb求解器，重點(diǎn)分析了ode45即龍格-庫(kù)塔-費(fèi)爾別格法。針對(duì)matlab軟件中常用的ODE45算法求解含有間斷特性的常微分方程會(huì)出現(xiàn)困難情況，采用件函數(shù)零點(diǎn)搜索法對(duì)ODE45算法進(jìn)行了改進(jìn)，通過(guò)算例驗(yàn)證了本文改進(jìn)算法可以有效求解某些含有間斷特性的常微分方程問(wèn)題。

參考文獻(xiàn)

[1] (美)Steven C.Chapra著; 唐玲燕, 田尊華譯.工程與科學(xué)數(shù)值方法的matlab實(shí)現(xiàn)[Z]. 北京: 清華大學(xué)出版社, 2009, 5.

[2] 肖田元, 范文慧. 計(jì)算機(jī)仿真(第2版)[Z]. 北京: 清華大學(xué)出版社, 2010, 2.

[3] 楊錄峰, 馬寧, 趙雙鎖. 一種變步長(zhǎng)和變階計(jì)算的自適應(yīng)數(shù)值積分算法[J]. 云南民族大學(xué)學(xué)報(bào)(自然科學(xué)版), 2011, 01.

[4] 張明會(huì), 高婷婷. 一個(gè)數(shù)值積分公式的推廣[J]. 四川文理學(xué)院學(xué)報(bào), 2011, 02.

[5] 許江浩, 陳志坤, 劉斌. 一個(gè)高精度數(shù)值積分公式[J]. 四川理工學(xué)院學(xué)報(bào)(自然科學(xué)版), 2011, 02.

[6] 邢誠(chéng), 王建強(qiáng), 賈志強(qiáng). 多種數(shù)值積分方法比較分析[J].城市勘測(cè), 2010, 01.

[7] 陳佩寧, 劉競(jìng). 用MATLAB求數(shù)值積分的方法[J]. 石家莊職業(yè)技術(shù)學(xué)院學(xué)報(bào), 2008, 06.

[8] 龍愛(ài)芳, 胡軍浩. 一個(gè)新的數(shù)值積分公式[J]. 數(shù)學(xué)理論與應(yīng)用, 2010, 04.

[9] 李海合, 王三福. 數(shù)值積分的一種改進(jìn)方法[J]. 甘肅聯(lián)合大學(xué)學(xué)報(bào)(自然科學(xué)版), 2009, 04.

[10] 張志涌等. 精通MATLAB6.5版(第1版)[Z]. 北京: 北京航空航天大學(xué)出版社, 2003.

Research on Improvement of a Numerical Integration Algorithm

WU Jia-hui
(College of Electronic Engineering, Xinjiang University, Urumqi 830047, China)

This paper analyzes the numerical integration algorithm of constant coefficient differential equation in MATLAB software, and solution method of differential equation with discontinuous characteristics. It focuses on the ode45 algorithm, and the error of solving the discontinuous problems, improves the zero search algorithm of conditional function, and verifies the effectiveness and applicability of the improved algorithm by comparing the results of using ode45, ode23, ode113, ode15s, ode23s, ode23t and ode23tb to solve the example with the results of using the improved ode algorithm to solve the example.

Numerical integration; Differential equation; Ode45

TP391.41

: A

10.3969/j.issn.1003-6970.2017.02.007

國(guó)家自然科學(xué)基金項(xiàng)目(51575469)

一種數(shù)值積分算法的改進(jìn)研究

本文著錄格式：吳佳惠. 一種數(shù)值積分算法的改進(jìn)研究[J]. 軟件，2017，38（2）：28-32