摘要：当前一些科学计算软件虽然能够利用数值积分得到复杂积分的近似解，但最多只能求解三重积分。为解决积分计算重数受限问题，基于累次积分的解析计算思路，结合递归算法，提出了一种任意重积分自适应递归式快速计算方法。从原理上介绍了算法思路和递归过程，给出了在MATLAB中运行的算法源程序代码。在算例和工程应用上，将该方法与现有方法进行了对比。结果表明了该方法的有效性，可使计算时间大大缩短，只需提供积分区间和被积函数即可求解，省去了将程序代码与积分重数进行匹配的编程步骤。方法操作简单，计算精度高，速度快，可满足实际工程需要。
关键词：数值积分； 多重积分； 递归算法； 工程应用；
Abstract:At present, some scientific calculation software can get approximate solution of complex integration through numerical integration. However, scientific calculation software can only solve triple integration at most. To solve the problem of limited multiplicity in integral calculation, an adaptive recursive fast calculation method for arbitrary multiple integrals was proposed. It was based on the analytic calculation idea of multiple integrals and recursive algorithm. The algorithm and recursion process of the method were introduced in principle, and the source code of the method was given in MATLAB. The method was compared with the existing methods by examples. The results indicate the effectiveness of the proposed method. This method can reduce the calculation time greatly. It can be solved only by providing integral interval and integrand function, eliminating the programming steps of matching program code with integral multiplicity. The method is simple in operation, and high in calculation accuracy and speed, which can meet the actual engineering needs.
Keyword:numerical integration; multiple integral; recursive algorithm; engineering application;
在各学科领域科学研究及工程应用中，如电力系统动态仿真[1]、准光学电磁场传播[2]、磁力驱动器磁力求解[3,4]等实际问题，积分计算是必不可少的工具和环节。传统解析积分Newton-Leibniz公式只适用于初等原函数存在情况下的定积分求解，而在实际应用中，往往会出现无初等原函数、原函数极为复杂以及无解析表达式这三种情况，如e-x2无初等原函数、x22x2+3??????√的原函数极为复杂、实验测得数据（无解析表达式）。数值积分通过收集积分区间内离散点及其对应函数值来逼近未知函数区间积分值，可以用来求解复杂定积分问题的近似值[5].
关于数值积分，中外学者还进行了各方面研究[6,7,8].张建斌等[9]设计出了基于MATLAB的数值积分求解器，实现了数值积分的可视化功能。陈付龙[10]讨论了二元数值积分的计算方法，并给出了部分C语言代码。唐珺等[11]将递归算法与自适应积分结合用于定积分计算，与辛普森积分相比减少了计算量。还有学者利用遗传算法[12]、神经网络算法[13,14,15]、粒子群算法[16]、混沌布谷鸟算法[17]计算数值积分。Angulo等[18]将数值积分应用于红细胞生长模型研究。尽管涉及数值积分的研究较多，但是对于三重以上积分数值计算的讨论相对较少。
此外，当前广泛使用的各类科学计算软件具备一定的数值积分计算功能，如著名的MATLAB、Maple、Mathematica等，其数值逼近方法大大加快了复杂积分计算速度。然而，这些科学计算软件自带函数所能计算的积分重数普遍也不高于三重。例如，MATLAB软件虽然具有强大的矩阵运算能力，能够很好地胜任数据计算这一烦琐工作，但能够进行数值积分的自带函数只有quad、dbquad、triplequad三种，分别为一重、二重、三重积分函数，也就是说MATLAB自带函数最多只能计算三重积分，远不能满足科学研究和工程应用需要。
为了提高科学计算软件的积分计算重数，根据多重积分的解析计算方式将累次积分的思想应用于数值积分中[19],现提出一种任意重积分自适应递归式快速计算方法（简称"自适应递归法"）。首先，详细介绍了任意重积分的算法思路，并将递归过程分为递推（划分积分区间）和回归（求解函数值）两部分进行阐述。其次，以MATLAB为例，该方法由积分区间、多重积分函数和被积函数表达式三部分组成。其中，多重积分函数采用递归算法和自适应辛普森公式，具有独立性和自适应性，即当积分重数改变时多重积分函数能够自动适应变化，而无需重新进行编程，避免了重复性操作。最后，结合算例，讨论了该方法与近似三重积分法[3]和辛普森法[19]三者之间的差别。
1 算法思路
递归算法在程序设计中被经常用到。程序编写过程中如果程序运行时又调用了其本身，称这种编写方式为递归[20].采用递归算法可以将一个复杂问题一层一层转化为多个较原问题更为简单的小问题求解，有效地缩短了程序编写的代码长度。同时递归需要有递归出口，当满足递归出口的约束条件时，递归结束。
该方法可用于求解已知被积函数表达式的任意重积分，得到该积分数值解。具体计算步骤如下。
设任意重积分为
I=∫b1a1∫b2a2…∫bnanf（X1,X2,…，Xn,P1,P2,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）中：第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个积分变量的积分区间划分为
Xi=[ai?ai+hi?ai+2hi?ai+bi2?bi?2hi?bi?hi?bi]???（4）
式（4）中：hi为第i个积分变量的积分步长，hi=1m（bi?ai）。
第i个积分变量的积分节点对应函数为
[fai?fai+hi?fai+2hi?fai+bi2?fbi?2hi?fbi?hi?fbi] （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）又可写成下面的形式，即
图2 各积分变量的积分节点对应函数树状图
Fig.2 Tree diagram of the corresponding functions of the integration nodes of each integration variable
函数的上角标代表赋值的不同积分节点，下角标与所求解积分变量的下角标相对应
∫b1a1∫b2a2…∫bnanf（X1,X2,…，Xn,P1,P2,P3,…）×dX1dX2…dXn=∫b1a1f1（X1,Yn-1）dX1（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）关系为
f1（X1,Yn-1）=∫b2a2f2（X2,Yn-2）dX2（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）的求解相当于求解定积分，即
式中：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所示，即前一次调用的积分变量的每个积分节点对应函数值是通过后一次调用的积分变量的所有积分节点对应函数值通过数值积分公式求解而得。
图3 积分节点对应函数的积分变量赋值树状图
Fig.3 Tree diagram of integral variable assignment of function corresponding to integral node
2.2 回归：求解函数值
整个递归过程中函数求解关系可以通过图2得知，函数中各积分变量的详细赋值如图3所示，图2中倒数第一行函数值均可通过图3中相应位置以及其祖先节点所在位置的积分节点求得。
结合图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定义如下：
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
∫0.50∫31∫0?π3∫e1∫10[e2x1x2sin（3x3）+x35/x4]dx1dx2dx3dx4dx5=???43（e?1）2+π12≈？3.6748572.
例2
∫π20∫20∫10∫40（3y2zcosx+m）dxdydxdm=??16+8π≈41.1327412.
例3 工程应用
磁力驱动器凭借极高的密封性和无接触特性，被广泛应用于制药、化工、石油、食品等行业[3].磁力驱动器的内外磁套之间的轴向磁力[3]被写为
Fa=A∫R2R1∫R4R3∫2π0∫2π0A1r1r2dr1dr2dαdβ（9）
式（9）中：r1为内磁套表面上任意一点的位置；r2为外磁套表面上任意一点的位置；α为内磁套表面上任意一点的角度；A=Br1Br24πμ0,其中β为外磁套表面上任意一点的角度，Br1为内磁套的剩余磁通密度，Br2为外磁套的剩余磁通密度，μ0为空间磁导率；A1表示为
A1=2δ0（A2+2δ02）3/2?δ0+δ2[A2+（δ0+δ2）2]3/2?δ0?δ1[A2+（δ0?δ1）2]3/2 （10）
式（10）中：δ0为内磁套和外磁套之间的相对位移，δ0=0.02 m; δ1为内磁套厚度；δ2为外磁套厚度；A2表示为
A2=（r2sinβ-r1sinα）2+（r2cosβ-r1cosα）2（11）
磁性套筒的材料参数和尺寸参考文献[3],如表1所示。
Table 1 Material parameters and dimensions of magnetic sleeve
表1 磁性套筒的材料参数和尺寸
以上三个算例的计算结果如表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时有数值解。综合以上算例，本文提出的自适应递归法的计算效率高于近似三重积分法和辛普森法，主要是由于以下两点：
（1）辛普森法中的步长h受积分区间的影响。h取值要小于各变量积分区间中的最小值，如例3中h<0.5时则无法求解。这种情况导致了积分区间划分的不合理，尤其是积分上下限差值较大时，过小的h增加了程序的运行时间。此外，辛普森法的积分区间划分方式固定，步长h恒定，增加了不必要的计算量，而自适应递归法巧妙地将各变量积分区间进行了独立划分，并单独设定误差限，实现了计算量的合理分配。自适应递归法可以在保证精度要求的同时有效地降低计算量，避免了辛普森法中因提升精度而导致计算时间呈指数形增长现象的发生。
Table 2 Calculation results of case 1, 2 and 3
表2 例1、例2和例3的计算结果
注：n为近似三重积分法[3]中沿dr方向的等分数；h为辛普森法[16]中的步长；ti为自适应递归法中各积分变量的误差容限；例3中δ0=20 mm; Inf表示经MATLAB计算后的结果为无穷大。
（2）自适应递归法的优势在于输入各变量积分上下限和被积函数表达式后，即可求解，适应性强，而近似三重积分和辛普森法缺少对积分重数的适应能力。近似三重积分仅能计算四重积分。采用辛普森法计算多重积分，积分重数与程序中的函数TNR_ITG[16]为一一对应关系，而不是自适应递归法中的多对一关系，导致函数TNR_ITG的内容需根据积分重数的变化而不断调整。积分重数与函数TNR_ITG的改动幅度为正相关。
5 结论
针对当前科学计算软件求解积分重数的有限性及传统解析积分的局限性，提出了一种任意重积分自适应递归式快速计算方法。该方法参考定积分计算方式，依次将各积分变量的积分区间进行自适应划分，当完全确定划分后各积分变量的积分节点时，再反向逐层求积，从而实现递归过程。
（1）基于该方法给出了MATLAB中的算法源程序代码，并进行了算例分析，验证了该方法在多重积分求解时的高效率。
（2）与现有方法的比较结果表明，提出的自适应递归法不仅计算时间短、结果精度高，效率高，而且凭借对积分重数的适应能力，省去了不断编写程序的烦琐过程，操作更为方便，具备通用性和普适性。
（3）该方法是直线往复运动磁力驱动器磁力求解问题的一种更简洁、快速且准确的计算方法，可应用于涉及积分运算的实际应用，具有一定的工程应用价值。
参考文献
[1] 戴汉扬，汤涌，宋新立，等。电力系统动态仿真数值积分算法研究综述[J].电网技术，2018,42（12）：3977-3984.Dai Hanyang,Tang Yong,Song Xinli,et al.Review on numerical integration algorithms for dynamic simulation of power system[J].Power System Technology,2018,42（12）：3977-3984.
[2] 梁彬，白明，金铭，等。瑞利索莫菲衍射与物理光学积分的应用比较[J].北京航空航天大学学报，2012（9）：1214-1218.Liang Bin,Bai Ming,Jin Ming,et al.Correlation between Rayleigh-Somerfield diffraction and physical optics integration[J].Journal of Beijing University of Aeronautics and Astronautics,2012（9）：1214-1218.
[3] Li J G,Tan Q C,Zhang Y Q,et al.Study on the calculation of magnetic force based on the equivalent magnetic charge method[J].Physics Procedia,2012,24（Part A）：190-197.
[4] Li J G,Tan Q C,Pei Y C,et al.Influence of non-coaxial alignment on interaction forces for the contactless magnetic driver of a reciprocating motion[J].International Journal of Applied Electromagnetic and Mechanics,2015,49（2）：169-177.
[5] 张雨浓，何良宇，金龙，等。一组一点外推的数值积分公式的提出与验证[J].数学的实践与认识，2013,43（8）：204-209.Zhang Yunong,He Liangyu,Jin Long,et al.Proposing and verification of a group of 1-node-extrapolated numerical intergration formulas[J].Mathematics in Practice and Theory,2013,43（8）：204-209.
[6] 夏炎，田波。Cowell数值积分器的变步长与自起步方法[J].科学技术与工程，2020,20（17）：6742-6747.Xia Yan,Tian Bo.The method of step size change and self-startup of Cowell integrator[J].Science Technology and Engineering,2020,20（17）：6742-6747.
[7] 张芮丹，王辉，吴昊，等。多参数积分化暂态高频电力变压器模型及其参数在线辨识方法[J].科学技术与工程，2019,19（21）：181-188.Zhang Ruidan,Wang Hui,Wu Hao,et al.Multi-parameter integral transient high-frequency power transformer model and its parameter online identification method[J].Science Technology and Engineering,2019,19（21）：181-188.
[8] 李亚东，张子军，杨凤田，等。某型电动飞机起飞爬升性能分析及飞行试验[J].科学技术与工程，2019,19（35）：364-369.Li Yadong,Zhang Zijun,Yang Fengtian,et al.Analysis of take-off and climb performance of a certain type of electric aircraft and flighting test[J].Science Technology and Engineering,2019,19（35）：364-369.
[9] 张建斌，赵静，许晓晴。基于MATLAB-GUI的数值积分界面设计[J].实验室研究与探索，2017,36（1）：13-16.Zhang Jianbin,Zhao Jing,Xu Xiaoqing.The numerical integration design based on MATLAB-GUI[J].Research and Exploration in Laboratory,2017,36（1）：13-16.
[10] 陈付龙。二元数值积分的计算方法[J].计算机工程与应用，2007,43（19）：32-34.Chen Fulong.Two-dimensional numerical value integral's computing method[J].Computer Engineering and Applications,2007,43（19）：32-34.
[11] 唐珺，杨道杰。自适应积分的递归实现研究[J].中国高新技术企业，2016（27）：13-16.Tang Jun,Yang Daojie.Research on recursive implementation of adaptive integral （in Chinese）[J].China High-Tech Enterprises,2016（27）：13-16.
[12] 赖志柱，张云艳。基于遗传算法求任意函数的数值积分[J].计算机工程与应用，2014,50（2）：54-57.Lai Zhizhu,Zhang Yunyan.Solving numerical integration based on genetic algorithms[J].Computer Engineering and Applications,2014,50（2）：54-57.
[13] 肖红，许少华，李盼池。基于数值积分的多聚合过程神经网络算法[J].信息与控制，2013,42（5）：608-621,624.Xiao Hong,Xu Shaohua,Li Panchi.Multi-aggregation process neural networks algorithm based on numerical integration[J].Information and Control,2013,42（5）：608-621,624.
[14] 罗玉雄，文卉。一种基于神经网络算法的数值积分方法[J].传感技术学报，2006,19（4）：1187-1189,1194.Luo Yuxiong,Wen Hui.A numerical integration approach using neural network algorithm[J].Chinese Journal of Sensors and Actuators,2006,19（4）：1187-1189,1194.
[15] 李盼池，施光尧。基于数值积分的离散过程神经网络算法及应用[J].系统工程理论与实践，2013,33（12）：3216-3222.Li Panchi,Shi Guangyao.Numerical integration-based discrete process neural networks algorithm and applications[J].System Engineering-Theory & Practice,2013,33（12）：3216-3222.
[16] 韦杏琼，周永权。基于粒子群算法的数值积分方法研究[J].微电子学与计算机，2009,26（7）：117-119.Wei Xingqiong,Zhou Yongquan.Numerical integral method research based on PSO[J].Microelectronics & Computer,2009,26（7）：117-119.
[17] 凌巍高，欧旭，罗德相。求解二重数值积分的混沌布谷鸟优化算法[J].微电子学与计算机，2014,31（8）：148-150,154.Ling Weigao,Ou Xu,Luo Dexiang.Two dimensional numerical integration based on chaotic cuckoo search optimization algorithm[J].Microelectronics & Computer,2014,31（8）：148-150,154.
[18] Angulo O,Crauste F,López-Marcos J C.Numerical integration of an erythropoiesis model with explicit growth factor dynamics[J].Journal of Computational and Applied Mathematics,2018,330:770-782.
[19] 朱振广。多重积分数值计算的一种递归方法[J].辽宁工学院学报，2004,24（1）：24-26,39.Zhu Zhenguang.A recurrence method of estimating multivariable integral[J].Journal of Liaoning Institute of Technology,2004,24（1）：24-26,39.
[20] 魏斌，马继辉，牛虎。基于递归算法的树型结构图的设计与实现[J].计算机应用与软件，2011,28（1）：96-98.Wei Bin,Ma Jihui,Niu Hu.Design and implementation of recursion algorithm-based tree structure diagram[J].Computer Applications and Software,2011,28（1）：96-98.
