2023-07-05 Matlab中的数值积分

在《计算方法》一书中有介绍基本的数值分析中的积分方法，我们这里重点关注Matlab是如何帮助我们快速计算数值积分。

1. integral簇函数

integral簇函数下包含integral, integral2, integral3三个函数，分别对应于一重积分，二重积分和三重积分。其基本形式为

1. %% 积分 
   

   ---

2. q = integral(fun,xmin,xmax) 
   

   ---

3. q = integral(fun,xmin,xmax,Name,Value) 
   

   ---

4.  
   

   ---

5. %% 二重积分 
   

   ---

6. q = integral2(fun,xmin,xmax,ymin,ymax) 
   

   ---

7. q = integral2(fun,xmin,xmax,ymin,ymax,Name,Value) 
   

   ---

8.  
   

   ---

9. %% 三重积分 
   

   ---

10. q = integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax) 
    

    ---

11. q = integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax,Name,Value) 
    

    ---

其中fun为目标函数，由@或者eval解析

1. function y = fsin(x) 
   

   ---

2. y = sin(5*x); 
   

   ---

3. end 
   

   ---

4. q = integral(@(x) eval('fsin(x)'), 0, 1) 
   

   ---

，xmin, xmax, ymin, ymax, zmin, zmax为积分在各个维度上的上下界，并且可由函数形式定义。Name, Value为参数对，由以下参数构成

1. AbsTol — 绝对误差容限（integral, integral2, integral3）
   绝对误差容限，指定为由 'AbsTol' 和非负实数组成的逗号分隔对组。integral 使用绝对误差容限来限制绝对误差容限估计值，即 |q – Q|，其中 q 为计算的积分值，Q 为（未知）确切值。如果减小绝对误差容限，integral 可以提供更多位小数的精度。
2. RelTol — 相对误差容限（integral, integral2, integral3）
   相对误差容限，指定为由 'RelTol' 和非负实数组成的逗号分隔对组。integral2 使用相对误差容限来限制相对误差容限估计值，即 |q – Q|/|Q|，其中 q 为计算的积分值，Q 为（未知）确切值。如果减小相对误差容限，integral2 可以为精度提供更多的有效数位。默认值为 1e-6。
3. ArrayValued — 数组值函数标志（integral）
   数组值函数标志，指定为以逗号分隔的对组，其中包含 'ArrayValued' 和数值或逻辑值 1 (true) 或 0 (false)。将该标志设为 true 或 1，以指示 fun 是接受标量输入并返回向量、矩阵或 N 维数组输出的函数。如
   1. %% x为一个数，fun为向量值函数（输入标量返回向量、矩阵等） 
      

      ---

   2. fun = @(x)sin([1 2; 3 4]*x); 
      

      ---

   3. q = integral(fun,0,1,'ArrayValued',true) 
      

      ---

   4. %% x为一个向量，fun为输入向量并返回向量的函数 
      

      ---

   5. fun = @(x)sin(5*x); 
      

      ---

   6. q = integral(fun,0,1,'ArrayValued',0) 
      

      ---

4. Waypoints — 积分路点（integral）
   似乎是指复平面上的积分路点（不慎清楚，复变函数已经早忘了）
5. Method — 积分法（integral2, integral3）
   对于大多数情况，integral2/3 使用 'tiled' 方法。当任何积分范围为无限时，它将使用 'iterated' 方法。这是默认方法。integral2/3 将积分区域变换为矩形形状，并根据需要将其划分为更小的矩形区域。积分范围必须是有限的。integral2/3 调用 integral 执行迭代积分。外积分在 xmin ≤ x ≤ xmax 范围内进行计算。内积分在 ymin(x) ≤ y ≤ ymax(x) 范围内进行计算。积分范围可以是无限的。

*1 在AbsTol和RelTol中Q是未知的，那么停止的判定应该怎么进行呢？在integral函数中，可以发现其调用integralCalc来计算积分，其中Q的设置为（errbnd = abs(err_ok) + norm(errsubs,1)），是某个误差项加上一个误差向量的1范数。论文Vectorized Adaptive Quadrature in Matlab和Weighted Quadrature by Change of Variable中似乎是用的approximate error bound，即找到一个近似误差上界。

*2 当您的函数在积分区域内不连续时，'iterated' 方法可能会更有效。但是，如果您在不连续点处拆分积分并计算多个积分的结果总和，则可实现最佳性能和精确度。

*3 在对非矩形区域积分时，如果 ymin、ymax（或两者）为函数句柄，则可实现最佳性能和精确度。避免将被积函数值设置为零来对非矩形区域积分。如果您必须执行此操作，请指定 'iterated' 方法。如果 ymin、ymax（或两者）为无边界函数，请使用 'iterated' 方法。

*4 含参函数的数值积分

1. fun = @(x,c) 1./(x.^3-2*x-c); 
   

   ---

2. q = integral(@(x) fun(x,5),0,2) 
   

   ---

2. quadgk和quad2d方法

这两个方法和上面的integral方法都可以对由具体表达式的函数进行数值积分，其区别是：

1. integral and quadgk use the same quadrature method as described in the reference paper
2. integral and quadgk use different default error tolerances for RelTol and AbsTol
3. quadgk has the MaxIntervalCount option, and uses a relatively small default value compared to integral. This means that with integral you will rarely if ever need to adjust this parameter (hence why it isn't there). This would only be with highly oscillatory integrands that are difficult to evaluate.
4. quadgk has an extra output that approxmates the absolute error. This lets you specify MaxIntervalCount, check the error in the calculation, and adjust MaxIntervalCount as needed.

另外 quadgk 和 integral 使用基本相同的积分方法。您通常应使用 integral 而不是 quadgk。不过，您可以使用 quadgk 执行以下操作：

a. 使用 errbnd 输出参数监控解的准确度。

b. 当 integral 警告达到最大区间数时，请为 MaxIntervalCount 指定较大的值。

如果奇异性不太强，quadgk 函数可对有限端点处奇异的函数求积分。例如，它可对在端点 c 处的行为类似于 log|x-c| 或 |x-c|p（其中 p >= -1/2）的函数求积分。如果函数在积分范围 [a b] 内的各点上有奇异性，则在奇异点为端点的子区间上将积分写入为积分和，通过 quadgk 计算它们，并将结果相加。

如果区间是无限区间（例如），则要存在 fun(x) 的积分，fun(x) 必须在 x 接近无限大时衰减，并且 quadgk 要求它快速衰减。

此外，对于二重积分则使用quad2d。

3. trapz 和 cumtrapz梯形数值积分

上述integral函数簇接受的参数是函数句柄，trapz接受向量或矩阵形式的点列，并根据点列用梯形法求数值积分：

1. Q = trapz(Y) 
   

   ---

2. Q = trapz(X,Y) 
   

   ---

3. Q = trapz(___,dim) 
   

   ---

具体来说

1. %% 当X，Y为向量时 
   

   ---

2. Y = [1 4 9 16 25]; 
   

   ---

3. Q = trapz(Y) % 求解[1,5]上Y的积分 
   

   ---

4. X = 1:0.25:2; 
   

   ---

5. Q = trapz(X,Y) % 求解[1,2]上Y的积分 
   

   ---

6.  
   

   ---

7. %% 当X为向量，Y为矩阵时 
   

   ---

8. % dim=1按列积分，dim=2按行积分 
   

   ---

9. X = [1 2.5 7 10]; 
   

   ---

10. Y = [5.2 7.7 9.6 13.2; 
    

    ---

11. 4.8 7.0 10.5 14.5; 
    

    ---

12. 4.9 6.5 10.2 13.8]; 
    

    ---

trapz还可以用于求解多重积分，只需要嵌套使用它

1. x = -3:.1:3;  
   

   ---

2. y = -5:.1:5;  
   

   ---

3. [X,Y] = meshgrid(x,y);%生成2D坐标网格 
   

   ---

4. F = X.^2 + Y.^2; %每个网格上的函数值 
   

   ---

5. I = trapz(y,trapz(x,F,2)) %先对x积分，再对y积分 
   

   ---

函数cumtrapz几乎和trapz一样，唯一的区别是cumtrapz会将每一步的计算结果都返回（矩阵形式）。

1. I = cumtrapz(y,trapz(x,F,2)) %先对x积分，再对y积分 
   

   ---