各种数值积分方法总结（龙贝格积分、自适应积分、高斯积分等）
1 (一重积分)常用的数值积分方法
1.1 牛顿-科茨（Newton-Cotes）积分公式
1.1.1 梯形法则（2点积分）
1.1.2-3 辛普森法则（3点积分和4点积分）
1.1.2 辛普森1/3法则（3点积分）
1.1.3 辛普森3/8法则（4点积分）
1.1.4 布尔法则（5点积分）
1.1.5牛顿-科茨积分公式总结
（穿插）积分改进思路1
1.2 复合积分（复合梯形积分、复合辛普森积分等）
1.2.1 复合梯形积分
1.2.2 复合辛普森积分
（穿插）积分改进思路2
1.3常用的高精度积分
1.3.1 龙贝格积分
（穿插）积分改进思路3
1.3.2 自适应积分
1.3.3 高斯积分
高斯-勒让德公式
高斯积分误差：
高斯-切比雪夫公式
广义积分：高斯-拉盖尔公式
广义积分：高斯-埃尔米特公式
2 多重积分方法
二重积分的辛普森公式
二重积分的蒙特卡洛方法
多重（高维）积分的蒙特卡洛方法

本文整理了各种类型的数值积分，从简单的梯形积分、辛普森积分到高精度的自适应积分、龙贝格积分和高斯积分；对于多重积分（高维积分），本文先简单介绍一下蒙特卡洛方法（适用于积分维数大于4的积分），后续会发帖专门针对多重积分方法做介绍。

（如有疏漏，欢迎指正，谢谢~）
1 (一重积分)常用的数值积分方法
常用积分方法包括梯形积分、辛普森积分、自适应积分、龙贝格积分和高斯积分等；其中自适应积分、龙贝格积分和高斯积分属于高精度积分。
1.1 牛顿-科茨（Newton-Cotes）积分公式
牛顿-科茨（Newton-Cotes）积分公式：



包含了梯形法则、辛普森法则、布尔法则等。
1.1.1 梯形法则（2点积分）
梯形法则积分计算公式：
其中，



于是得：



梯形公式要计算两个点的函数值f(a)和f(b)。所以说是两点积分



梯形法则的误差为：


1.1.2-3 辛普森法则（3点积分和4点积分）
辛普森法则分为辛普森1/3法则和辛普森3/8法则，就是用二次/三次插值曲线来代替直线进行积分。


1.1.2 辛普森1/3法则（3点积分）
辛普森1/3法则法则积分计算公式：


1.1.3 辛普森3/8法则（4点积分）
辛普森3/8法则法则积分计算公式：


1.1.4 布尔法则（5点积分）
布尔法则与辛普森法则类似，只是更高阶的多项式积分，其公式为：


1.1.5牛顿-科茨积分公式总结


注意：1.对于3点积分与4点积分具有相同的误差阶；5点积分与6点积分具有相同的误差阶。对于点数更多的积分也一样。因此优先使用奇数个点的积分。

2.对于点数大于等于9的牛顿-科茨积分，求积公式的稳定性得不到保证（具体原因可以找数值积分的数看看，本文源自《现代数值积分》第5章 数值积分与数值微分），积分不能保证收敛，因此实际计算中一般不采用高阶牛顿-科茨积分公式。（处于效率考虑，一般也很少用超过5个点的牛顿-科茨积分公式）。
（穿插）积分改进思路1
牛顿-科茨积分公式思想：根据积分点构造近似的插值多项式，进而利用多项式积分表示原积分。

对于高次的多项式插值，会出现龙格现象（对于高次的多项式插值，插值多项式会出现不收敛的现象，成为龙格现象，想深入了解可以去查询相关资料。并且次数越高，越容易出现不收敛现象）。
针对龙格现象，更好的选择切比雪夫节点来进行插值（想深入了解可以去查询相关资料，本文部分源自《现代数值积分》第3章 多项式插值与样条插值）。


由于切比雪夫节点的特殊性，对于高阶插值的可以收到较好的效果，因此，对于“牛顿-科茨积分公式总结”中提到的点数大于等于9的牛顿-科茨积分，可以采用切比雪夫计算节点，然后构造插值多项式，再计算积分（此方法实现和推广比较麻烦，此处不做扩展，感兴趣的可以自己下去推导实现）。
1.2 复合积分（复合梯形积分、复合辛普森积分等）
上述主要介绍的是单个区间的积分。为了计算更准确，我们可以将区间划分为若干个小的区间，然后对每个小区间进行积分，然后对各个小区间的积分求和。由于对于单个小区间的积分可分为梯形法则，辛普森法则，和布尔法则，因此复合积分的每个小区间的积分可以是这些方法中的任意一种。

下面主要讲一下 复合梯形积分和复合辛普森积分。
1.2.1 复合梯形积分
小区间的划分又分为等长区间和不等长区间两种划分方法。

不等长区间，复合梯形积分计算公式：



区间不等长情况应用场景很有限。
等长区间，复合梯形积分计算公式：



等长区间，复合梯形积分计算误差：


1.2.2 复合辛普森积分
同样的，复合辛普森积分区间的划分也可以分为等长区间和不等长区间两种划分方法。

由于不等长区间应用很少，因此直接介绍等长区间，复合辛普森1/3积分计算公式：


等长区间，复合辛普森1/3积分计算误差：


由于3点积分（辛普森1/3积分）与4点积分（辛普森3/8积分）具有相同的误差阶，因此一般优先使用复合辛普森1/3积分，复合辛普森3/8积分应用极少，因此本文不做介绍。感兴趣的同学可以下去自己推一下复合辛普森3/8积分计算公式，应该也很简单。
（穿插）积分改进思路2
可以看出，将区间划分为更小的区间可以有效的提高积分计算的精度。那么我们是不是可以将积分区间不断地细分，从而来得到更精确的积分结果？理论上是可以的！！

方法1：（区间折半法）

a.计算2个区间的积分，用复合辛普森1/3法则（当然，也可以用梯形法则或者布尔法则等都可以，但是推荐复合辛普森1/3法则）;

b.将两个区间等分为4个区间，用复合辛普森1/3法则计算四个子区间的积分和；

c.比较a和b两步计算的结果相对误差是否接近一个很小的数:(Ia - Ib)/Ib<e。如果相对误差满足要求，则结束，否则继续细分区间，用复合辛普森1/3法则计算积分，直到满足要求为止。
方法2：改进方法1
	a.对每个子区间一分为二后，比较两个区间的积分值与一个区间的积分值是否相近;
	b.如果相近则返回该值；否则对每个子区间执行a。（递归实现）

上面说，理论上是可以的！！理论上是可以的！！但是实际不行！！！

因为计算机浮点数运算是存在误差的。当区间数量很大的时候，浮点数计算的舍入误差（不理解的可以自行去学习了解一下）变得很大，会限制积分的提高。下面就是一个很好的例子。


1.3常用的高精度积分
上述区间二分的思想是很有用的。在龙贝格积分和自适应积分中都会找到它的影子。
1.3.1 龙贝格积分
龙贝格积分的思想，及推导过程：











类推：联合两个误差为O(h^4)的积分，改进后可以得到一个误差为O(h ^ 6)的积分：



由此产生龙贝格积分：




（穿插）积分改进思路3
龙贝格积分是用梯形公式推导得到的，那么我们是不是可以用辛普森公式做类似的推导，也能到达一套新的公式呢？答案是：理论上是可以的！！
由辛普森1/3公式推导出来的加入误差修正的积分计算式：



那么问题就来了。

我们的计算机表达的double值是有范围的，而对于辛普森外，采用递归，递归8次就是16^8 = 2 ^ 32。计算机递归计算的结果很容易就产生溢出了。除非采用多精度浮点数计算才行，而采用多精度浮点数可能存在一些计算效率的问题。。。更深入的问题我就不继续深入探讨了（有兴趣的可以试试）。。。好了，总之来说，用辛普森推导出来的方法不是很实用。还是梯形龙贝格积分经典好使。
1.3.2 自适应积分
龙贝格积分中有两点：

a。每次都会将所有的步长折半；（那么我们能不能只对计算误差较大的那个区间进行折半成两个子区间再积分呢？）

b。辛普森公式推导的方法，由于计算机double表示范围有限没用上。
下面都会在自适应积分中派上用场。自适应积分思想：（递归方法）

step1.辛普森1/3法则计算区间[a,b]的积分；

step2.区间步长折半，复合辛普森1/3法则计算区间[a,b]的积分；

如果两次误差小于给定的误差限，停止执行，返回计算结果

；

如果两次误差大于给定的误差限，对两个子区间分别执行step1和step2。
伪代码：


1.3.3 高斯积分
高斯积分应用比较广泛的是高斯-勒让德公式，此外还有高斯-切比雪夫公式；以及区间【0，正无穷】的广义积分 高斯-拉盖尔公式；以及区间【负无穷，正无穷】的广义积分 高斯-埃尔米特公式。
高斯-勒让德公式
介绍高斯积分的文档就很多了，本文不细讲推导的过程，给出计算方法：



多点高斯-勒让德公式：





高斯勒让德公式的区间为【-1，1】.对于任意区间为【a，b】；需要做一下简单的变量代换，


高斯积分误差：
高斯勒让德公式的区间为【-1，1】，误差：



对于任意区间为【a，b】，误差：





由上式继续往下推导，可以得出，当区间长度【b-a】> 2(n+1)时，误差会发散。因此对于高斯点为n的高斯积分，积分区间长度不要大于2(n+1)，否则积分结果很可能会不准确。

此外，高斯积分的误差还与被积函数的光顺性有关，被积函数光顺性越差，误差越大。详细可参考《现代数值积分》第5章 数值积分与数值微分，P139页）。
高斯积分的优点：高斯积分时给定节点数下代数精度最高的求积公式。
高斯积分的不足：高斯积分每次改变积分点的个数，所以的积分点的函数值都需要重新计算。
高斯-切比雪夫公式
高斯点为切比雪夫节点。
（未完待续。。。）
广义积分：高斯-拉盖尔公式
广义积分：高斯-埃尔米特公式
2 多重积分方法
二重积分的辛普森公式
二重积分的蒙特卡洛方法
多重（高维）积分的蒙特卡洛方法

数值积分的基本思想
由积分中值定理可知，在积分区间 
 内存在一点 
 ，成立
 式的几何意义即为：底为 
 而高为 
 的矩形的面积恰等于所求曲边梯形的面积 
 。因此，要想求出 
 式左端积分，我们只需要知道三个值： 
 即可，这里 
 是显然的，问题在于 
 的具体位置一般是不清楚的，从而 
 未知。我们暂且将 
 称为区间 
 上的平均高度，我们的目标就是寻求一种求出平均高度 
 的算法，这样我们的问题也就解决了。
首先我们想到的是用区间两端点的“高度” 
 与 
 的算术平均值作为平均高度 
 的近似值，从而导出下述求积公式
 式即为我们我们熟知的
梯形公式。写出MATLAB函数为
function I = Trapezoid(f, a, b, n)
% Trapezoid Formula
% (b-a)/2 * (f(a)+f(b))
h = (b - a)/n;
x = a:h:b;
I = sum((x(2:end) - x(1:end-1))/2.*(f(x(2:end)) + f(x(1:end-1))));
end
而若用区间中点 
 的“高度” 
 近似取代平均高度 
 ，则可导出
中矩形公式（简称矩形公式）
其MATLAB函数为
function I = Rectangle(f, a, b, n)
% Rectangle Formula
% (b-a) * f((a+b)/2)
h = (b - a)/n;
x = a:h:b;
I = sum((x(2:end) - x(1:end-1)).*(f((x(2:end) + x(1:end-1))./2)));
end
总结：更一般地，我们可以在区间 
 上适当选取某些节点 
 ，然后用 
 的加权平均得到平均高度 
 的近似值，这样构造出来的求积公式具有下列通式
其中 
 称为
求积节点， 
 称为
求积系数，也称作伴随节点 
 的
权。
这种数值积分方法通常称为机械求积，其优势在于将积分求值问题归结为被积函数值的计算，很适合在计算机上使用。
代数精度
由于数值积分方法是近似方法，为保证精度，自然希望我们的求积公式能够对“尽可能多”的函数准确成立，这便提出了代数精度的概念。
如果某个求积公式对于次数不超过 
 的多项式均能够准确成立，但对于 
 次多项式就不准确成立，则称该求积公式具有
 次代数精度
（或代数精确度）。
前面提到的梯形公式 
 和矩形公式 
 均具有一次代数精度。
求积公式的收敛性与稳定性
在求积公式 
 中，若 
 其中 
 ，则称求积公式 
 是
收敛的。
对任给 
 ，若 
 ，只要 
 就有 
 成立，其中 
 为计算 
 时产生误差 
 后实际得到的值，即 
 ，则称求积公式 
 是
稳定的。
定理 
若求积公式 
中系数 
，则此求积公式是稳定的。
牛顿-柯特斯(Newton-Cotes)公式
将积分区间 
 划分为 
 等份，此时步长为 
 ，选取等距节点 
 构造出的插值型求积公式
称为牛顿-柯特斯(Newton-Cotes)公式，其中 
 称为
柯特斯系数。可以通过下式确定
当 
 时， 
 这时的求积公式便是之前的梯形公式 
 。
当 
 时， 
 此时对应的求积公式便是
辛普森(Simpson)公式： 
 写成MATLAB函数为
function I = Simpson(f, a, b, n)
% Simpson Formula
% (b-a)/6 * (f(a) + 4*f((a+b)/2) + f(b))
h = (b - a)/n;
x = a:h:b;
I = sum((x(2:end) - x(1:end-1))/6.*(f(x(1:end-1)) + 4*f((x(2:end)+x(1:end-1))./2) + f(x(2:end))));
end
当 
 时的牛顿-柯特斯公式则特别称为
柯特斯公式，其形式为 
 这里 
 。
其MATLAB代码为
function I = Cotes(f, a, b, n)
% Cotes Formula
% (b-a)/90 * (7*f(x0) + 32*f(x1) + 12*f(x2) + 32*f(x3) + 7*f(x4))
h = (b - a)/n;
x = a:h:b;
I = sum((x(2:end) - x(1:end-1))/90.*(7*f(x(1:end-1)) + 32*f((x(2:end)+3*x(1:end-1))./4) + 12*f(2*(x(2:end)+x(1:end-1))./4) + 32*f((3*x(2:end)+x(1:end-1))./4) + 7*f(x(2:end))));
end
复合求积公式
由于牛顿-柯特斯公式在 
 时不具有稳定性，故不可能再通过提高阶的方法来提高求积精度。
复合求积法便是通过把积分区间分成若干个子区间（通常是等分），再在每个子区间上使用低阶求积公式的方法来提高精度的。
其实细心的同学可以发现，前面我所展示的MATLAB代码使用的便是对每个求积公式使用复合法完成的。常用的有复合梯形公式和复合辛普森公式。
龙贝格(Romberg)求积公式
前面介绍的复合求积方法可提高求积精度，如若精度仍不够，则可通过将步长逐次减半的方式来提高精度。如对复合梯形公式可导出其递推公式
其中 
 表示在 
 基础上步长 
 减半后的积分值。
定理 
 设 
 ，则有
  其中
 
 ，
系数 
 与 
 无关。
从 
 式可以看出， 
 是 
 阶，若用 
 代替 
 ，有
再用 
 乘 
 式再减去 
 式后再除以 
 得
这里 
 是与 
 无关的系数。
从 
 式可看出，用近似积分值 
 ，其误差阶为 
 ，显然误差阶是提高了。类似这种将计算 
 的近似值的误差阶由 
 提高到 
 的方法称为
外推算法，也称为理查德森(Richardson)外推算法。
这是“数值分析”中一个重要的技巧，只要真值与近似值的误差能够表示成 
 的幂级数，如 
 式所示，都可以使用外推算法，提高精度。
龙贝格(Romberg)算法
 次外推加速为
余项为
此方法常称为理查德森外推加速方法。
设用 
 表示二分 
 次后求得的梯形值，且以 
 表示序列 
 的 
 次加速值，则由递推公式 
 可得
上次则称为龙贝格求积算法，计算过程如下：
取 
 ，求 
 。令 
 ( 
 记区间 
 的二分次数)。
求梯形值 
 ，即按递推公式 
 计算 
 。
求加速值。
若 
 (预先给定的精度)，则终止计算，并且 
 ；否则令 
 转 
 继续计算。
MATLAB代码为
function [T,n,I] = Romberg(f,a,b,eps)
% Romberg Integral Algorithm
%--输入项：
% f: 被积函数
% a,b: 积分上下限
% eps: 积分精度
%--输出项：
% T: 逐次积分表
% n: 迭代次数
% R: 最后积分值

if nargin < 4
    eps = 1e-8; % 默认精度
end

h = b-a; T(1,1) = h*(f(a)+f(b))/2;
n = 1; k = 0; err = 1;
while err > eps
    k = k+1; S = 0;
    for ii = 1:n
        x = a+h*(2*ii-1)/2;
        S = S+f(x);
    end
    T(k+1,1) = (T(k,1)+h*S)/2; % 递推公式 n=128
    for m = 1:k
        T(k+1,m+1) = (4^m*T(k+1,m)-T(k,m)) / (4^m-1); % Richardson Extrapolation
    end
    err = abs(T(k+1,k+1)-T(k+1,k)); % 积分精度
    n=2*n; h = h/2;
end
I=T(k+1,k+1);
end
自适应积分方法
复合求积方法通常适用于被积函数变化不太大的积分，如果在求积区间中被积函数变化很大，有的部分函数值变化剧烈，另一部分却变化平缓。这时如果统一将区间等分再用复合求积公式计算积分将会导致计算量很大，我们想实现的是在满足误差要求的前提下，对变化剧烈部分将区间细分，而平缓部分则可使用大步长，也即针对被积函数在区间上不同情形采用不同的步长，使得再满足精度前提下积分计算工作量尽可能小。其算法技巧是在不同区间上预测被积函数变化的剧烈程度确定相应步长。这就是自适应积分方法的基本思想。
下面是一个基于复合辛普森法的自适应积分算法的MATLAB代码：
function I = Self_adaptiveSimpson(a, b, epsilon, S)
% Self-adaptive Simpson Integral Method
% a,b: 区间上下限
% epsilon: 误差精度
% S: 初始 Simpson 积分值

h = (b-a)/2;
S1 = h/6*(f(a)+4*f(a+h/2)+f(a+h));
S2 = h/6*(f(a+h)+4*f(a+3*h/2)+f(b));
if abs(S-S1-S2) <= epsilon
    I = (S1+S2)+(S1+S2-S)/15;
else
    I = Self_adaptiveSimpson(a,(a+b)/2,epsilon/2,S1)+Self_adaptiveSimpson((a+b)/2,b,epsilon/2,S2);
end

%定义被积函数
    function y = f(x)
        y = 1./(x.^2);
    end
end
下面我们计算积分 
 ，可输入下列语句调用上述MATLAB函数进行计算
clear
f1 = @(x) 1./(x.^2);
S = ComplexSimpson(f1,0.2,1,1); % 初始 Simpson 积分值
I = Self_adaptiveSimpson(0.2, 1, 0.02, S)
结果为
I =

   4.000059571596276
高斯求积公式
下面研究带权积分 
 （ 
 为权函数）的求积公式
下面看定义
如果求积公式 
 具有 
 次代数精度，则称其节点 
 为
高斯点，相应公式 
 称为
高斯型求积公式。
这里主要包括四种高斯型求积公式，即
高斯-勒让德求积公式
高斯-切比雪夫求积公式
高斯-拉盖尔求积公式
高斯-埃尔米特求积公式
多重积分
这里主要说一下多重积分计算的思路，归结为一句话，就是：多重积分化累次积分，再由里到外使用数值积分公式进行求积计算。

	

导读：介绍数值计算的基本思想及常用方法。
基本思想数值计算就把原来空间及时间坐标上连续的物理场(速度场、温度场、压力场等)，用一系列有限个离散点(节点)上的值的集合来代替，通过离散方程建立这些离散点上变量值之间的关系，求解这些离散方程，最终获得所求解变量的近似值。具体流程如下图所示。数值方法1.有限差分法(FDM，finite difference method)
将求解区域用与坐标轴平行的一系列网格线条的交点所组成的点的集合来代替。
每个节点上，将控制方程中每一个导数用相应的差分表达式来代替，从而在每个节点上形成一个代数方程。
每个方程中都包括了本节点及其附近一些节点上的未知值，求解这些代数方程就获得了所需的数值解。
缺陷：数值解的守恒性无法保证、复杂几何的适应性。2.有限容积法(FVM，finite volume method)
计算区域划分为一系列控制容积，每个控制容积都有一个节点来代表。
通过将守恒型的控制方程对控制容积做积分来导出离散方程
导出过程中，需要对界面上被求函数本身及其一阶导数作出假定，这种构成方式就是有限容积法中的离散格式。
优势：保证了守恒性，并且离散方程具有明确的物理意义。3.有限元法(FEM，finite element method)
将计算域分成一系列元体，每个元体取数个点作为节点，然后通过对控制方程做积分来获得离散方程。
与有限容积法不同，FEM需要选定形状函数，通过节点上被求变量来表示形状函数，积分前将其带入控制方程。
控制方程积分前需要乘上一个权函数。
对不规则区域适应性好，但同时也加大了计算量。4.有限分析法(FAM，finite analytic method)
与有限差分法类似，用一系列网格线将区域离散；
不同的是，每个节点与相邻的4个网格(二维)组成计算单元，也就是说一个计算单元由一个中心节点与8个邻点组成。
在计算单元中将控制方程的非线性项局部线性化，并对单元上未知函数的变化型线做假设，把所选定的线型表达式中的系数和常数项用单元边界节点上未知变量值表示，这样计算单元中的被求问题就可以转化为第一类边界条件下的一个定解问题，可以得到解析解。
利用这个解析解，得出该单元中点及边界上的8个邻点上未知值间的代数方程，此即为单元中点的离散方程。
有限分析法的系数不像有限容积法中那样有明确的物理意义，对不规则区域的适应性也较差。

二重积分的数值方法
《数值分析课程设计》
报  告
专业：
学号：
学生姓名：
指导教师：
一、题目
数值积分中二重积分探究。
二、理论
数值积分就是用数值方法近似计算定积分。其原理很简单，就是将积分核用插值多项式替代，用多项式的结果近似定积分的值。
一般常用的方法是，将积分区间等分为个子区间，即取步长
，
子区间端点为(k=0,1,…，n)，在每个子区间上套用插值积分公式，再将个区间的结果累加起来。
比较常用的有梯形公式，是在每个子区间上用1阶多项式(即直线段)近似并积分的结果：
另外一种在实际应用中很受欢迎的方法是，在每个子区间上用2阶多项式(即抛物线)近似并积分，得到著名的辛普森(Simpson)公式：
其中。
三、方法、算法与程序设计
Ⅰ.辛普森公式求二重积分
考虑二重积分，它是曲面与平面区域R围成的体积，对于矩形区域，可将它写成累次积分
。
若用复合辛普森公式，可分别将,分成N,M等份，步长
,，
先对积分，应用复合辛普森公式，令，，则
从而得
。对每个积分再分别用复合辛普森公式即可求得积分值。MATLAB程序见附录1，MATLAB中自带自适应辛普森公式dblquad()，对于变量区域同样适用。
对于变量区域，写成累次积分的形式：
进行数值计算的表达式为：
上面的表达式中、表示权重，取决于一维积分方法。我们常用复合辛普森公式，先对内积分进行计算，在计算外积分，与矩形区域情况基本一致。
Ⅱ高斯求积公式求二重积分
在高斯求积公式中，若取权函数，区间为，则得公式
。
勒让德多项式是区间上的正交多项式，因此，勒让德多多项式的零点就是求积公式的高斯点。
若取的零点做节点构造求积公式
；
若取的零点构造求积公式
；
当时，求积公式为
同样先用高斯求积公式求内积分，再求外积分，可得二重积分值。
四、算例、应用实例
算例：
计算二重积分。
(1)若区域，试分别用复合辛普森公式(取n=4)及高斯求积公式(取n=4)求积分。
(2)若区域用复合辛普森公式(取n=4)求积分。
解：
(1)
=
对各个积分应用复合辛普森公式。
也可应用MATLAB中的函数进行计算，程序见附录2。
先将区域变换为区域，其中
，等价于，有
。
对于取时的高斯求积公式节点及系数，即
，，，，
，，，
用的高斯积分公式计算积分I，
(2)
，
等分为4等份，对应值为的值，用
节点应用辛普森公式对内积分求积，再用复合辛普森公式对外积分求积，
也可用MATLAB中的函数实现，结果和程序如下(附录3)。
五、参考文献
【1】  数值分析   李庆扬，王朝能，易大义    清华大学出版社
【2】  数值分析课程设计   陈越，童若锋    浙江大学出版社
【3】  MATLAB教程    张志涌       北京航空航天大学出版社
六、附录
附录1：
function q=DblSimpson(f,a,A,b,B,m,n)if(m==1 && n==1)??%辛普森公式q=((B-b)*(A-a)/9)*(subs(sym(f),findsym(sym(f)),{a,b})+...subs(sym(f),findsym(sym(f)),{a,B})+...subs(sym(f),findsym(sym(f)),{A,b})+...subs(sym(f),findsym(sym(f)),{A,B})+...4*subs(sym(f),findsym(sym(f)),{(A-a)/2,b})+...4*subs(sym(f),findsym(sym(f)),{(A-a)/2,B})+...4*subs(sym(f),findsym(sym(f)),{a,(B-b)/2})+...   4*subs(sym(f),findsym(sym(f)),{A,(B-b)/2})+...16*subs(sym(f),findsym(sym(f)),{(A-a)/2,(B-b)/2}));else?? %复合辛普森公式q=0;for i=0:n-1for j=0:m-1x=a+2*i*(A-a)/2/n;y=b+2*j*(B-b)/2/m;?????????