一種任意重積分自適應(yīng)遞歸式快速計算方法

裴永臣， 關(guān)景晗， 王佳煒

(吉林大學(xué)機(jī)械與航空航天工程學(xué)院， 長春 130025)

關(guān)于數(shù)值積分，中外學(xué)者還進(jìn)行了各方面研究[6-8]。張建斌等[9]設(shè)計出了基于MATLAB的數(shù)值積分求解器，實現(xiàn)了數(shù)值積分的可視化功能。陳付龍[10]討論了二元數(shù)值積分的計算方法，并給出了部分C語言代碼。唐珺等[11]將遞歸算法與自適應(yīng)積分結(jié)合用于定積分計算，與辛普森積分相比減少了計算量。還有學(xué)者利用遺傳算法[12]、神經(jīng)網(wǎng)絡(luò)算法[13-15]、粒子群算法[16]、混沌布谷鳥算法[17]計算數(shù)值積分。Angulo等[18]將數(shù)值積分應(yīng)用于紅細(xì)胞生長模型研究。盡管涉及數(shù)值積分的研究較多，但是對于三重以上積分?jǐn)?shù)值計算的討論相對較少。

此外，當(dāng)前廣泛使用的各類科學(xué)計算軟件具備一定的數(shù)值積分計算功能，如著名的MATLAB、Maple、Mathematica等，其數(shù)值逼近方法大大加快了復(fù)雜積分計算速度。然而，這些科學(xué)計算軟件自帶函數(shù)所能計算的積分重數(shù)普遍也不高于三重。例如，MATLAB軟件雖然具有強(qiáng)大的矩陣運(yùn)算能力，能夠很好地勝任數(shù)據(jù)計算這一煩瑣工作，但能夠進(jìn)行數(shù)值積分的自帶函數(shù)只有quad、dbquad、triplequad三種，分別為一重、二重、三重積分函數(shù)，也就是說MATLAB自帶函數(shù)最多只能計算三重積分，遠(yuǎn)不能滿足科學(xué)研究和工程應(yīng)用需要。

為了提高科學(xué)計算軟件的積分計算重數(shù)，根據(jù)多重積分的解析計算方式將累次積分的思想應(yīng)用于數(shù)值積分中[19]，現(xiàn)提出一種任意重積分自適應(yīng)遞歸式快速計算方法(簡稱“自適應(yīng)遞歸法”)。首先，詳細(xì)介紹了任意重積分的算法思路，并將遞歸過程分為遞推(劃分積分區(qū)間)和回歸(求解函數(shù)值)兩部分進(jìn)行闡述。其次，以MATLAB為例，該方法由積分區(qū)間、多重積分函數(shù)和被積函數(shù)表達(dá)式三部分組成。其中，多重積分函數(shù)采用遞歸算法和自適應(yīng)辛普森公式，具有獨(dú)立性和自適應(yīng)性，即當(dāng)積分重數(shù)改變時多重積分函數(shù)能夠自動適應(yīng)變化，而無需重新進(jìn)行編程，避免了重復(fù)性操作。最后，結(jié)合算例，討論了該方法與近似三重積分法[3]和辛普森法[19]三者之間的差別。

1 算法思路

遞歸算法在程序設(shè)計中被經(jīng)常用到。程序編寫過程中如果程序運(yùn)行時又調(diào)用了其本身，稱這種編寫方式為遞歸[20]。采用遞歸算法可以將一個復(fù)雜問題一層一層轉(zhuǎn)化為多個較原問題更為簡單的小問題求解，有效地縮短了程序編寫的代碼長度。同時遞歸需要有遞歸出口，當(dāng)滿足遞歸出口的約束條件時，遞歸結(jié)束。

該方法可用于求解已知被積函數(shù)表達(dá)式的任意重積分，得到該積分?jǐn)?shù)值解。具體計算步驟如下。

設(shè)任意重積分為

P3,…)dX1dX2…dXn

(1)

式(1)中：a1,a2,…,an為積分下限；b1,b2,…,bn為積分上限；X1,X2,…,Xn為積分變量；P1,P2,P3,…為可變參數(shù)。

(1)根據(jù)實際生產(chǎn)或研究需求，輸入被積函數(shù)為

f(X1,X2,…,Xn,P1,P2,P3,…)

(2)

(2)輸入問題求解每個積分變量的積分參數(shù)(積分下限、積分上限和誤差容限)，并以積分參數(shù)矩陣的形式給出。積分參數(shù)矩陣A為n行3列矩陣，即

(3)

式(3)中：第i行元素ai、bi、ti為第i個積分變量的積分下限、積分上限、誤差容限。

(3)設(shè)置遞歸出口，若積分參數(shù)矩陣A行數(shù)為0，則遞歸終止。

(4)讀取積分參數(shù)矩陣A中第1個積分變量的積分參數(shù)。

(5)刪除積分參數(shù)矩陣A中第1個積分變量的積分參數(shù)，即刪去第一行。

(6)根據(jù)步驟(4)中讀取的第1個積分變量的積分下限、積分上限，將積分區(qū)間劃分成m份，即積分節(jié)點(diǎn)為m+1個。由步驟(4)、步驟(5)知，如若重復(fù)執(zhí)行步驟(4)時，相當(dāng)于第i次執(zhí)行時讀取的是積分參數(shù)矩陣A中第i個積分變量的積分參數(shù)。第i個積分變量的積分區(qū)間劃分為

bi-2hibi-hibi]

(4)

第i個積分變量的積分節(jié)點(diǎn)對應(yīng)函數(shù)為

(5)

(7)再次執(zhí)行步驟(4)～步驟(6)，直到滿足步驟(3)中約束條件，遞歸終止，此時利用數(shù)值積分公式回歸求解。

以上步驟的執(zhí)行流程如圖1所示。

圖1 積分流程示意圖Fig.1 Schematic diagram of integration process

2 遞歸過程

整個遞歸過程分為兩步：第一步，遞推，依次劃分各積分變量的積分區(qū)間，從積分變量X1到積分變量Xn；第二步，回歸，反向求解各積分節(jié)點(diǎn)處函數(shù)值，從函數(shù)fn到函數(shù)f1。

2.1 遞推：劃分積分區(qū)間

求解任意重積分與定積分?jǐn)?shù)值計算方法的具體實施過程相似。定積分?jǐn)?shù)值計算方法先確定被劃分后的積分節(jié)點(diǎn)及其所對應(yīng)的函數(shù)值，然后利用數(shù)值積分公式得到數(shù)值解。此外，該方法與解析方法雖然在求解形式上有所不同，但兩者的求解思路存在相同之處，即遞歸方法，求解n重積分時，先求解n-1重積分；而求解n-1重積分時，可先求n-2重積分；依次遞推直到最后求解的是定積分。于是任意重積分問題可轉(zhuǎn)化為定積分問題來求解。將任意重積分看作求X1的積分，將其他未賦值積分變量暫時看作已知量，用Yn-1表示，此時的被積函數(shù)用f1(X1,Yn-1)表示。于是式(1)又可寫成下面的形式，即

總是會有不同程度路面路基問題出現(xiàn)在道路橋梁在交付使用后，受到外界環(huán)境的影響，會降低道路橋梁使用年限，影響橋梁的質(zhì)量并加重之前質(zhì)量問題。所以，要求有關(guān)技術(shù)人員定期檢查破損的路基路面情況，維修排水系統(tǒng)，重視好后期的保養(yǎng)及維護(hù)，及時上報并做好相應(yīng)記錄，便于市政工程團(tuán)隊人員的修復(fù)?？偠灾?，需全面提升路基路面施工質(zhì)量，加強(qiáng)后期維護(hù)及保養(yǎng)。

(6)

此時，任意重積分轉(zhuǎn)化成定積分問題，利用定積分?jǐn)?shù)值計算方法即可得到該問題數(shù)值解。然而，被劃分后的積分節(jié)點(diǎn)對應(yīng)的函數(shù)值中含有未賦值積分變量X2,X3,…,Xn，無法確定f1(X1,Yn-1)的值。因此，只有當(dāng)f1(X1,Yn-1)的值確定后，才能通過數(shù)值積分公式得到該問題數(shù)值解。由于f1(X1,Yn-1)中含有未賦值積分變量，故再次執(zhí)行步驟(4)～步驟(6)。f1(X1,Yn-1)中X1已知，所以求解f1(X1,Yn-1)相當(dāng)于求解n-1重積分，可以看作是X1被賦值后求X2的積分，其他未賦值積分變量X3,X4,…,Xn暫時看作已知量，用Yn-2表示。此時的被積函數(shù)表達(dá)式用f2(X2,Yn-2)表示。f1(X1,Yn-1)與f2(X2,Yn-2)關(guān)系為

(7)

此時，n-1重積分轉(zhuǎn)化成定積分問題，利用定積分?jǐn)?shù)值計算方法即得到該問題數(shù)值解。然而，被劃分后的積分節(jié)點(diǎn)對應(yīng)的函數(shù)值中含有未賦值積分變量X3,X4,…,Xn，無法確定f2(X2,Yn-2)的值。因此，只有當(dāng)f2(X2,Yn-2)的值確定后，才能通過數(shù)值積分公式得到該問題數(shù)值解。由于f2(X2,Yn-2)中含有未賦值的積分變量，故再次執(zhí)行步驟(4)～步驟(6)。

以此類推，fi(Xi,Yn-i)的求解相當(dāng)于求解n-i重積分，……，fn-1(Xn-1,Y1)的求解相當(dāng)于求解定積分，即

(8)

式中：Yn-i為Xi+1,Xi+2,…,Xn的集合。

X1的各積分節(jié)點(diǎn)對應(yīng)函數(shù)中X1的賦值均為該積分節(jié)點(diǎn)，故可將f1(X1,Yn-1)分成m+1組，然后再根據(jù)X2的積分節(jié)點(diǎn)將之前每組函數(shù)再次分成了m+1組。以此類推，直到Xn的積分節(jié)點(diǎn)將Xn-1每組函數(shù)再次分成了m+1組。根據(jù)積分變量數(shù)和各積分變量的積分節(jié)點(diǎn)數(shù)，最終求解的函數(shù)共有(m+1)n個。前后兩次調(diào)用的函數(shù)求解關(guān)系如圖2所示，即前一次調(diào)用的積分變量的每個積分節(jié)點(diǎn)對應(yīng)函數(shù)值是通過后一次調(diào)用的積分變量的所有積分節(jié)點(diǎn)對應(yīng)函數(shù)值通過數(shù)值積分公式求解而得。

函數(shù)的上角標(biāo)代表賦值的不同積分節(jié)點(diǎn)，下角標(biāo)與所求解積分變量的下角標(biāo)相對應(yīng)圖2 各積分變量的積分節(jié)點(diǎn)對應(yīng)函數(shù)樹狀圖Fig.2 Tree diagram of the corresponding functions of the integration nodes of each integration variable

2.2 回歸：求解函數(shù)值

整個遞歸過程中函數(shù)求解關(guān)系可以通過圖2得知，函數(shù)中各積分變量的詳細(xì)賦值如圖3所示，圖2中倒數(shù)第一行函數(shù)值均可通過圖3中相應(yīng)位置以及其祖先節(jié)點(diǎn)所在位置的積分節(jié)點(diǎn)求得。

圖3 積分節(jié)點(diǎn)對應(yīng)函數(shù)的積分變量賦值樹狀圖Fig.3 Tree diagram of integral variable assignment of function corresponding to integral node

結(jié)合圖2、圖3，圖2中只有函數(shù)fn中無未賦值積分變量。也就是說，僅當(dāng)函數(shù)fn的值確定后，才能通過數(shù)值積分公式進(jìn)一步求解函數(shù)fn-1的值。函數(shù)fn-1的值確定后，可根據(jù)積分變量Xn-1的積分節(jié)點(diǎn)，求得函數(shù)fn-2的值。以此類推，最終求得該任意重積分?jǐn)?shù)值解。

3 MATLAB算法實現(xiàn)

在MATLAB中，本文方法所用的數(shù)值積分基層算法是quad函數(shù)自帶的自適應(yīng)辛普森公式。程序共有三個腳本，文件名分別為mm.m、multiquad.m、myfun.m，其中mm.m為主程序文件，multiquad.m為多重積分函數(shù)multiquad定義文件(在計算中無需改動)，myfun.m為多變量被積函數(shù)myfun定義文件，運(yùn)行時只需運(yùn)行主程序即可。

3.1 主程序說明

該方法主程序如下：

clear;clc

P1=* ; P2=* ; … % P1, P2,…是多變量被積函數(shù)myfun必要的可變參數(shù)

A=[a1,b1,t1;a2,b2,t2;…;ai,bi,ti;…;an,bn,tn] %積分參數(shù)矩陣

Ip=multiquad([],[],A,@myfun,P1,P2,…) %調(diào)用多重積分函數(shù)multiquad

3.2 多重積分函數(shù)multiquad定義

主程序中調(diào)用的多重積分函數(shù)multiquad，定義為I=multiquad([],[],A,@myfun,P1,P2,…)。前兩個參數(shù)用于程序運(yùn)行過程中的數(shù)據(jù)傳遞，主程序中初始設(shè)定為空矩陣，myfun表示為多變量被積函數(shù)。多重積分函數(shù)multiquad以MATLAB自帶積分函數(shù)quad為基層算法求解核心，因此具有自適應(yīng)性，通過多次調(diào)用達(dá)到求解任意重積分目的。多重積分函數(shù)multiquad定義如下：

function [I]=multiquad(x,X0,A,int_fun,varargin)

a=A(1,1); %積分下限

b=A(1,2); %積分上限

t=A(1,3); %誤差容限

A(1,:)=[]; %刪除積分參數(shù)矩陣A第一行

if size(A,1)==0

I=zeros(size(x));

for i=1:length(x)

I(i)=quad(@f_multiquad,a,b,t,[],[X0;x(i)],int_fun,varargin{:}); %使用函數(shù)quad進(jìn)行計算

end

elseif size(A,1)>0

if length(x)==0

I=quad(@multiquad, a,b,t,[],X0,A,int_fun,varargin{:}); %使用函數(shù)quad進(jìn)行計算

else

I=zeros(size(x));

for i=1:length(x)

I(i)=quad(@multiquad,a,b,t,[],[X0;x(i)],A,int_fun,varargin{:}); %使用函數(shù)quad進(jìn)行計算

end

end

end

函數(shù)f_multiquad為多重積分函數(shù)multiquad的內(nèi)部調(diào)用函數(shù)，用于計算積分區(qū)間劃分后的各積分節(jié)點(diǎn)對應(yīng)函數(shù)值。

function [y]=f_multiquad(x,X0,int_fun,varargin)

if size(X0,1)>0

X=[X0*ones(size(x));x];

else

X=[x];

end

for i=1:size(X,2)

y(i)=feval(int_fun,X(:,i),varargin{:}); %計算指定函數(shù)在某點(diǎn)的函數(shù)值

end

3.3 多變量被積函數(shù)myfun定義

多變量被積函數(shù)myfun括號中X為積分變量X1,X2,…,Xn的集合。多變量被積函數(shù)myfun接受X,P1,P2,…并返回被積函數(shù)值y。主函數(shù)中的多變量被積函數(shù)myfun定義如下：

function [y]=myfun(X,P1,P2,…)

X1=X(1,:); %積分變量X1

X2=X(2,:); %積分變量X2

?

Xn=X(n,:); %積分變量Xn

y=f(X1,X2,…,Xn,P1,P2,P3,…); %被積函數(shù)f(X1,X2,…,Xn,P1,P2,P3,…)

4 算例分析

選取三種具有代表性的函數(shù)進(jìn)行程序驗證，并將自適應(yīng)遞歸法與解析法，近似三重積分法[3]和辛普森法[19]在計算時間和精度上進(jìn)行了對比。

例1

例2

16+8π≈41.132 741 2。

例3 工程應(yīng)用

磁力驅(qū)動器憑借極高的密封性和無接觸特性，被廣泛應(yīng)用于制藥、化工、石油、食品等行業(yè)[3]。磁力驅(qū)動器的內(nèi)外磁套之間的軸向磁力[3]被寫為

(9)

(10)

式(10)中：δ0為內(nèi)磁套和外磁套之間的相對位移，δ0=0.02 m；δ1為內(nèi)磁套厚度；δ2為外磁套厚度；A2表示為

A2=(r2sinβ-r1sinα)2+(r2cosβ-r1cosα)2

(11)

磁性套筒的材料參數(shù)和尺寸參考文獻(xiàn)[3]，如表1所示。

表1 磁性套筒的材料參數(shù)和尺寸Table 1 Material parameters and dimensions of magnetic sleeve

以上三個算例的計算結(jié)果如表2所示。

通過表2得知，例1和例2中，辛普森法(h=0.05)與自適應(yīng)遞歸法(ti=10-6)精度相同，辛普森法(h=0.01)與自適應(yīng)遞歸法(ti=10-7)精度相同，但計算時間差距明顯。例1中，自適應(yīng)遞歸法的計算時間與辛普森法相比，分別縮短了69.6%和99.8%。例2中，自適應(yīng)遞歸法的計算時間與辛普森法相比，分別縮短了87.8%和99.9%。例3中，近似三重積分法計算磁力驅(qū)動器的軸向磁力耗時嚴(yán)重，而辛普森法和自適應(yīng)遞歸法能夠?qū)⒂嬎銜r間減低到20 s以內(nèi)。值得注意的是，例3中辛普森法的h只在0.01時有數(shù)值解。綜合以上算例，本文提出的自適應(yīng)遞歸法的計算效率高于近似三重積分法和辛普森法，主要是由于以下兩點(diǎn)：

表2 例1、例2和例3的計算結(jié)果Table 2 Calculation results of case 1, 2 and 3

(1)辛普森法中的步長h受積分區(qū)間的影響。h取值要小于各變量積分區(qū)間中的最小值，如例3中h<0.5時則無法求解。這種情況導(dǎo)致了積分區(qū)間劃分的不合理，尤其是積分上下限差值較大時，過小的h增加了程序的運(yùn)行時間。此外，辛普森法的積分區(qū)間劃分方式固定，步長h恒定，增加了不必要的計算量，而自適應(yīng)遞歸法巧妙地將各變量積分區(qū)間進(jìn)行了獨(dú)立劃分，并單獨(dú)設(shè)定誤差限，實現(xiàn)了計算量的合理分配。自適應(yīng)遞歸法可以在保證精度要求的同時有效地降低計算量，避免了辛普森法中因提升精度而導(dǎo)致計算時間呈指數(shù)形增長現(xiàn)象的發(fā)生。

(2)自適應(yīng)遞歸法的優(yōu)勢在于輸入各變量積分上下限和被積函數(shù)表達(dá)式后，即可求解，適應(yīng)性強(qiáng)，而近似三重積分和辛普森法缺少對積分重數(shù)的適應(yīng)能力。近似三重積分僅能計算四重積分。采用辛普森法計算多重積分，積分重數(shù)與程序中的函數(shù)TNR_ITG[16]為一一對應(yīng)關(guān)系，而不是自適應(yīng)遞歸法中的多對一關(guān)系，導(dǎo)致函數(shù)TNR_ITG的內(nèi)容需根據(jù)積分重數(shù)的變化而不斷調(diào)整。積分重數(shù)與函數(shù)TNR_ITG的改動幅度為正相關(guān)。

5 結(jié)論

針對當(dāng)前科學(xué)計算軟件求解積分重數(shù)的有限性及傳統(tǒng)解析積分的局限性，提出了一種任意重積分自適應(yīng)遞歸式快速計算方法。該方法參考定積分計算方式，依次將各積分變量的積分區(qū)間進(jìn)行自適應(yīng)劃分，當(dāng)完全確定劃分后各積分變量的積分節(jié)點(diǎn)時，再反向逐層求積，從而實現(xiàn)遞歸過程。

(1)基于該方法給出了MATLAB中的算法源程序代碼，并進(jìn)行了算例分析，驗證了該方法在多重積分求解時的高效率。

(2)與現(xiàn)有方法的比較結(jié)果表明，提出的自適應(yīng)遞歸法不僅計算時間短、結(jié)果精度高，效率高，而且憑借對積分重數(shù)的適應(yīng)能力，省去了不斷編寫程序的煩瑣過程，操作更為方便，具備通用性和普適性。

(3)該方法是直線往復(fù)運(yùn)動磁力驅(qū)動器磁力求解問題的一種更簡潔、快速且準(zhǔn)確的計算方法，可應(yīng)用于涉及積分運(yùn)算的實際應(yīng)用，具有一定的工程應(yīng)用價值。