一種任意重積分自適應遞歸式快速計算方法

裴永臣， 關景晗， 王佳煒

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

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

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

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

1 算法思路

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

該方法可用于求解已知被積函數表達式的任意重積分，得到該積分數值解。具體計算步驟如下。

設任意重積分為

P3,…)dX1dX2…dXn

(1)

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

(1)根據實際生產或研究需求，輸入被積函數為

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

(2)

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

(3)

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

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

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

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

(6)根據步驟(4)中讀取的第1個積分變量的積分下限、積分上限，將積分區間劃分成m份，即積分節點為m+1個。由步驟(4)、步驟(5)知，如若重復執行步驟(4)時，相當于第i次執行時讀取的是積分參數矩陣A中第i個積分變量的積分參數。第i個積分變量的積分區間劃分為

bi-2hibi-hibi]

(4)

第i個積分變量的積分節點對應函數為

(5)

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

以上步驟的執行流程如圖1所示。

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

2 遞歸過程

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

2.1 遞推：劃分積分區間

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

總是會有不同程度路面路基問題出現在道路橋梁在交付使用后，受到外界環境的影響，會降低道路橋梁使用年限，影響橋梁的質量并加重之前質量問題。所以，要求有關技術人員定期檢查破損的路基路面情況，維修排水系統，重視好后期的保養及維護，及時上報并做好相應記錄，便于市政工程團隊人員的修復。總而言之，需全面提升路基路面施工質量，加強后期維護及保養。

(6)

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

(7)

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

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

(8)

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

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

函數的上角標代表賦值的不同積分節點，下角標與所求解積分變量的下角標相對應圖2 各積分變量的積分節點對應函數樹狀圖Fig.2 Tree diagram of the corresponding functions of the integration nodes of each integration variable

2.2 回歸：求解函數值

整個遞歸過程中函數求解關系可以通過圖2得知，函數中各積分變量的詳細賦值如圖3所示，圖2中倒數第一行函數值均可通過圖3中相應位置以及其祖先節點所在位置的積分節點求得。

圖3 積分節點對應函數的積分變量賦值樹狀圖Fig.3 Tree diagram of integral variable assignment of function corresponding to integral node

結合圖2、圖3，圖2中只有函數fn中無未賦值積分變量。也就是說，僅當函數fn的值確定后，才能通過數值積分公式進一步求解函數fn-1的值。函數fn-1的值確定后，可根據積分變量Xn-1的積分節點，求得函數fn-2的值。以此類推，最終求得該任意重積分數值解。

3 MATLAB算法實現

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

3.1 主程序說明

該方法主程序如下：

clear;clc

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

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

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

3.2 多重積分函數multiquad定義

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

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

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

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

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

A(1,:)=[]; %刪除積分參數矩陣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{:}); %使用函數quad進行計算

end

elseif size(A,1)>0

if length(x)==0

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

else

I=zeros(size(x));

for i=1:length(x)

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

end

end

end

函數f_multiquad為多重積分函數multiquad的內部調用函數，用于計算積分區間劃分后的各積分節點對應函數值。

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{:}); %計算指定函數在某點的函數值

end

3.3 多變量被積函數myfun定義

多變量被積函數myfun括號中X為積分變量X1,X2,…,Xn的集合。多變量被積函數myfun接受X,P1,P2,…并返回被積函數值y。主函數中的多變量被積函數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,…); %被積函數f(X1,X2,…,Xn,P1,P2,P3,…)

4 算例分析

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

例1

例2

16+8π≈41.132 741 2。

例3 工程應用

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

(9)

(10)

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

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

(11)

磁性套筒的材料參數和尺寸參考文獻[3]，如表1所示。

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

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

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

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

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

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

5 結論

針對當前科學計算軟件求解積分重數的有限性及傳統解析積分的局限性，提出了一種任意重積分自適應遞歸式快速計算方法。該方法參考定積分計算方式，依次將各積分變量的積分區間進行自適應劃分，當完全確定劃分后各積分變量的積分節點時，再反向逐層求積，從而實現遞歸過程。

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

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

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