吳佳惠

一種數值積分算法的改進研究

吳佳惠

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

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

數值積分；微分方程；ODE45

0 引言

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

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

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

1 RKF45算法改進

RKF（Runge–Kutta–Fehlberg）45即4階5級龍格-庫塔-費爾別格法，也matlab軟件中ODE45算法，具體表達式如下：

式中的各系數見表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公式中的各系數Tab.1 The Coefficients in Fehlberg Formula

由文獻[2]可知，當用經典的數值積分法（包括單步法和多步法）進行數值求解時，右端函數的連續(xù)性對保證仿真精度是十分重要的。單步法要求在一個步長內右端函數連續(xù)，K階多步法則要求在K步內右端函數連續(xù)，為了使仿真結果達到一定精度，人們往往采用變步長方法。然而，采用變步長方法也會遇到困難。

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

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

在經典RKF45算法中加入條件函數零點搜索法，可以使其更好地求解某些含有間斷點的微分方程求解問題，其代碼為：

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 案例驗證

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

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

忽略管中水的體積，計算100 s的時間內容器水位線與時間的函數關系并繪制圖形，假設空容器的初始條件為y(0)=0，計算時使用下面的參數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。通過ode23，ode113，ode15 s，ode23 s，ode23 t和ode23 tb來求解結果如下圖所示。圖的結果明顯是不正確的，除了最初的注滿周期之外，水位線看起來在再次達到y(tǒng)high以前就開始下降。類似地，在排水的時候，水位線還沒有降至ylow，虹吸管就關上了。

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

當噴泉停止噴射時，分子的第二項變成0。為此，我們可以在方程中加入一個新的無量綱的變量siphon，當噴泉停止時siphon等于0，當噴泉工作時siphon等于1。

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

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的取值必須在函數調用之間被保持，所以它被申明為一個全局變量。雖然，我們不鼓勵使用全局變量（特別是對于大型的程序），但是，它在當前的問題中很有用。下面的腳本用內置ode45函數積分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)’);

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

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

由文獻[1]可知，是因為ode在虹吸管開關時并不連續(xù)，例如，當容器注水時，導數僅僅依賴于固定的流入物和取值恒為某個參數。然而，一旦水位線到達yhigh，流出開關導通，導數迅速地變?yōu)榱硪粋€值，雖然Matlab所用的自適應步長程序對于許多問題微分方程求解問題十分有效，但是，它們在處理這類間斷問題時往往會失效。我們用ODE45和改進的ODE45算法求解的結果如下，可以看出ODE45同樣出現這種困難，是因為ODE在虹吸管開關時并不連續(xù)，但是采用改進的ODE45算法，就沒有出現問題。因為，通過改進經典ODE45算法，加入了條件函數零點搜索法，所以，使得ODE45算法可以求解某些含有間斷特性的微分方程求解問題。

圖3 ODE45和改進的ODE45求解結果Fig.3 Ode45 and improved ode45 solution results

3 結論

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

參考文獻

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

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

[3] 楊錄峰, 馬寧, 趙雙鎖. 一種變步長和變階計算的自適應數值積分算法[J]. 云南民族大學學報(自然科學版), 2011, 01.

[4] 張明會, 高婷婷. 一個數值積分公式的推廣[J]. 四川文理學院學報, 2011, 02.

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

[6] 邢誠, 王建強, 賈志強. 多種數值積分方法比較分析[J].城市勘測, 2010, 01.

[7] 陳佩寧, 劉競. 用MATLAB求數值積分的方法[J]. 石家莊職業(yè)技術學院學報, 2008, 06.

[8] 龍愛芳, 胡軍浩. 一個新的數值積分公式[J]. 數學理論與應用, 2010, 04.

[9] 李海合, 王三福. 數值積分的一種改進方法[J]. 甘肅聯合大學學報(自然科學版), 2009, 04.

[10] 張志涌等. 精通MATLAB6.5版(第1版)[Z]. 北京: 北京航空航天大學出版社, 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

國家自然科學基金項目(51575469)

一種數值積分算法的改進研究

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