23.3 多变量函数 ¶

Octave 包含几个用于计算多变量函数积分的函数。通常，这一过程可以通过先创建一个对 x 进行积分的函数 f，再将该函数对 y 进行积分来实现。下面的例子演示了如何手动完成这一过程，被积函数为：

f(x, y) = sin(pi*x*y) * sqrt(x*y)

其中 x 和 y 的取值范围均为 0 到 1。

下面的示例使用 quadgk 来执行二重积分。（注意，除 quad 外，任何一维求积函数都可以这样使用，因为 quad 是用 Fortran 编写的，不能被递归调用。）

function q = g(y)
  q = ones (size (y));
  for i = 1:length (y)
    f = @(x) sin (pi*x.*y(i)) .* sqrt (x.*y(i));
    q(i) = quadgk (f, 0, 1);
  endfor
endfunction

I = quadgk ("g", 0, 1)
      ⇒  0.30022

上述算法已在 dblquad 函数中实现，用于两个变量的积分。该过程的三维等效版本在 triplequad 中实现，用于三个变量的积分。例如，上面的结果可以通过调用 dblquad 来复现，如下所示。

I = dblquad (@(x, y) sin (pi*x.*y) .* sqrt (x.*y), 0, 1, 0, 1)
      ⇒  0.30022

q = dblquad (f, xa, xb, ya, yb) ¶

q = dblquad (f, xa, xb, ya, yb, tol) ¶

q = dblquad (f, xa, xb, ya, yb, tol, quadf) ¶

q = dblquad (f, xa, xb, ya, yb, tol, quadf, …) ¶

数值计算 f 的二重积分。

f 是一个函数句柄、内联函数或包含要评估的函数名称的字符串。函数 f 必须具有形式 z = f(x,y)，其中 x 是一个向量，y 是一个标量。它应返回一个与 x 长度和方向相同的向量。

xa、ya 和 xb、yb 分别是 x 和 y 的积分下限和上限。底层积分器决定是否接受无穷边界。

可选参数 tol 定义用于积分每个子积分的绝对容差。默认值为 1e-6。

可选参数 quadf 指定要使用的底层积分器函数。除 quad 外任何函数都可用，默认值为 quadcc。

额外的参数将直接传递给 f。若要使用 tol 或 quadf 的默认值，可以传入 ':' 或空矩阵 ([])。

另请参阅： integral2、integral3、triplequad、quad、quadv、quadl、quadgk、quadcc、trapz。

q = triplequad (f, xa, xb, ya, yb, za, zb) ¶

q = triplequad (f, xa, xb, ya, yb, za, zb, tol) ¶

q = triplequad (f, xa, xb, ya, yb, za, zb, tol, quadf) ¶

q = triplequad (f, xa, xb, ya, yb, za, zb, tol, quadf, …) ¶

数值计算 f 的三重积分。

f 是一个函数句柄、内联函数或包含要评估的函数名称的字符串。函数 f 必须具有形式 w = f(x,y,z)，其中 x 或 y 中有一个是向量，其余输入为标量。它应返回一个与 x 或 y 长度和方向相同的向量。

xa、ya、za 和 xb、yb、zb 分别是 x、y 和 z 的积分下限和上限。底层积分器决定是否接受无穷边界。

可选参数 tol 定义用于积分每个子积分的绝对容差。默认值为 1e-6。

可选参数 quadf 指定要使用的底层积分器函数。除 quad 外任何函数都可用，默认值为 quadcc。

额外的参数将直接传递给 f。若要使用 tol 或 quadf 的默认值，可以传入 ':' 或空矩阵 ([])。

另请参阅： integral3、integral2、dblquad、quad、quadv、quadl、quadgk、quadcc、trapz。

上述用于求积的递归算法被称为"迭代法"。函数 quad2d 中实现了一种不同的二维积分方法。此函数执行"分块法"积分，将积分域细分为矩形块并在这些块上分别进行积分。域在需要达到所需数值精度的区域被进一步细分。对于某些函数，这种方法可能比上述其他函数中使用的二维迭代法更快。

q = quad2d (f, xa, xb, ya, yb) ¶

q = quad2d (f, xa, xb, ya, yb, prop, val, …) ¶

[q, err, iter] = quad2d (…) ¶

数值计算 f 的二维积分，在由 xa、xb、ya、yb 定义的二维域上使用自适应求积进行分块积分。此外，ya 和 yb 可以是 x 的标量函数，从而允许在非矩形区域上进行积分。

f 是一个函数句柄、内联函数或包含要评估的函数名称的字符串。函数 f 必须具有形式 z = f(x,y)，并且所有运算必须向量化，使得 x 和 y 能接受相同大小的数组输入并返回数组输出。（可以假设 x 和 y 要么是相同大小的数组，要么是标量。）底层积分器会将积分点数组输入到 f 和/或使用内部向量展开来加速计算，但如果 f 不限于逐元素运算，则可能出错。对于被积函数无法避免这种情况时，下面描述的 ("Vectorized") 选项可能会产生更可靠的结果。

可以使用 property, value 对指定额外的可选参数。有效属性包括：

AbsTol

定义求积的绝对误差容限。默认值为 1e-10（single 类型为 1e-5）。

RelTol

定义求积的相对误差容限。默认值为 1e-6（single 类型为 1e-4）。

MaxFunEvals

对向量化函数 f 的最大调用次数。默认值为 5000。

Singular

启用/禁用在积分域边缘削弱奇点的变换。默认值为 true。

Vectorized

启用或禁用向量化积分。设置为 false 将强制 Octave 在调用被积函数时仅使用标量输入，这使得可以对尚未向量化或仅接受 x 或 y 标量值的 f(x,y) 进行积分。默认值为 true。请注意，这是通过使用 arrayfun 函数包装 f(x,y) 来实现的，这可能会显著降低计算速度。

FailurePlot

如果 quad2d 在达到 MaxFunEvals 之前未能收敛到所需的误差容限，则创建显示仍需要细化的区域的图。默认值为 false。

自适应求积用于最小化误差估计，直到满足以下条件：

        error <= max (AbsTol, RelTol*|q|)

可选输出 err 是积分误差 abs (q - I) 的近似上限，其中 I 是积分的精确值。可选输出 iter 是对函数 f 的向量化函数调用次数。

示例 1：在 x-y 平面上积分矩形区域

f = @(x,y) 2*ones (size (x));
q = quad2d (f, 0, 1, 0, 1)
  ⇒   q =  2

结果是一个体积，对于此常值被积函数，即为 Length * Width * Height。

示例 2：在 x-y 平面上积分三角形区域

f = @(x,y) 2*ones (size (x));
ymax = @(x) 1 - x;
q = quad2d (f, 0, 1, 0, ymax)
  ⇒   q =  1

结果是一个体积，对于此常值被积函数 f = 2，即为三角形面积乘以高度，即 1/2 * Base * Width * Height。

示例 3：在正方形区域上积分非向量化函数

f = @(x,y) sinc (x) * sinc (y));
q = quad2d (f, -1, 1, -1, 1)
  ⇒   q =  12.328  (不正确)
q = quad2d (f, -1, 1, -1, 1, "Vectorized", false)
  ⇒   q =  1.390 (正确)
f = @(x,y) sinc (x) .* sinc (y);
q = quad2d (f, -1, 1, -1, 1)
  ⇒   q =  1.390  (正确)

第一个结果不正确，因为函数 f 中 sinc 函数之间的非逐元素运算符在 quad2d 使用的内部积分数组之间产生了意外的矩阵乘法。在第二个结果中，将 "Vectorized" 设置为 false 强制 quad2d 执行标量内部运算来计算积分，以增加约 20 倍的计算时间为代价得到了正确的数值结果。在第三个结果中，使用逐元素乘法运算符对被积函数 f 进行向量化，在不增加计算时间的情况下获得了正确的结果。

编程注意事项：如果积分区域内有奇点，最好将积分拆分，并将奇点置于边界上。

已知 MATLAB 不兼容性：如果未指定容差，且任何积分限为 single 类型，则 Octave 的积分函数会自动降低上述默认的绝对和相对误差容限。如果需要更严格的容差，则必须明确指定。MATLAB 则无论积分限的类型如何，都会保留适用于 double 输入的更严格容差。

参考文献：L.F. Shampine, "MATLAB program for quadrature in 2D", Applied Mathematics and Computation, Vol. 1, pp. 266–274, 2008.

另请参阅： integral2、dblquad、integral、quad、quadgk、quadv、quadl、quadcc、trapz、integral3、triplequad。

最后，函数 integral2 和 integral3 作为通用的二维和三重积分函数提供。它们会在迭代法和分块法之间自动选择，并在非矩形积分域上使用 dblquad 和 triplequad。

q = integral2 (f, xa, xb, ya, yb) ¶

q = integral2 (f, xa, xb, ya, yb, prop, val, …) ¶

[q, err] = integral2 (…) ¶

数值计算 f 的二维积分，在由 xa、xb、ya、yb 定义的二维域上使用自适应求积（标量可以是有限值或无穷大）。此外，ya 和 yb 可以是 x 的标量函数，从而允许在非矩形区域上进行积分。

f 是一个函数句柄、内联函数或包含要评估的函数名称的字符串。函数 f 必须具有形式 z = f(x,y)，并且所有运算必须向量化，使得 x 和 y 能接受相同大小的数组输入并返回数组输出。（可以假设 x 和 y 要么是相同大小的数组，要么是标量。）底层积分器会将积分点数组输入到 f 和/或使用内部向量展开来加速计算，但如果 f 不限于逐元素运算，则可能出错。对于被积函数无法避免这种情况时，下面描述的 ("Vectorized") 选项可能会产生更可靠的结果。

可以使用 "property", value 对指定额外的可选参数。有效属性包括：

AbsTol

定义求积的绝对误差容限。默认值为 1e-10（single 类型为 1e-5）。

RelTol

定义求积的相对误差容限。默认值为 1e-6（single 类型为 1e-4）。

Method

指定要使用的二维积分方法，有效值为 "auto"（默认）、"tiled" 或 "iterated"。使用 "auto" 时，Octave 将选择 "tiled" 方法，除非某个积分限为无穷大。

Vectorized

启用或禁用向量化积分。设置为 false 将强制 Octave 在调用被积函数时仅使用标量输入，这使得可以对尚未向量化或仅接受 x 或 y 标量值的 f(x,y) 进行积分。默认值为 true。请注意，这是通过使用 arrayfun 函数包装 f(x,y) 来实现的，这可能会显著降低计算速度。

自适应求积用于最小化误差估计，直到满足以下条件：

        error <= max (AbsTol, RelTol*|q|)

err 是积分误差 abs (q - I) 的近似上限，其中 I 是积分的精确值。

示例 1：在 x-y 平面上积分矩形区域

f = @(x,y) 2*ones (size (x));
q = integral2 (f, 0, 1, 0, 1)
  ⇒   q =  2

结果是一个体积，对于此常值被积函数，即为 Length * Width * Height。

示例 2：在 x-y 平面上积分三角形区域

f = @(x,y) 2*ones (size (x));
ymax = @(x) 1 - x;
q = integral2 (f, 0, 1, 0, ymax)
  ⇒   q =  1

结果是一个体积，对于此常值被积函数 f = 2，即为三角形面积乘以高度，即 1/2 * Base * Width * Height。

示例 3：在正方形区域上积分非向量化函数

f = @(x,y) sinc (x) * sinc (y));
q = integral2 (f, -1, 1, -1, 1)
  ⇒   q =  12.328  (不正确)
q = integral2 (f, -1, 1, -1, 1, "Vectorized", false)
  ⇒   q =  1.390 (正确)
f = @(x,y) sinc (x) .* sinc (y);
q = integral2 (f, -1, 1, -1, 1)
  ⇒   q =  1.390  (正确)

第一个结果不正确，因为函数 f 中 sinc 函数之间的非逐元素运算符在 integral2 使用的内部积分数组之间产生了意外的矩阵乘法。在第二个结果中，将 "Vectorized" 设置为 false 强制 integral2 执行标量内部运算来计算积分，以增加约 20 倍的计算时间为代价得到了正确的数值结果。在第三个结果中，使用逐元素乘法运算符对被积函数 f 进行向量化，在不增加计算时间的情况下获得了正确的结果。

编程注意事项：如果积分区域内有奇点，最好将积分拆分，并将奇点置于边界上。

已知 MATLAB 不兼容性：如果未指定容差，且任何积分限为 single 类型，则 Octave 的积分函数会自动降低上述默认的绝对和相对误差容限。如果需要更严格的容差，则必须明确指定。MATLAB 则无论积分限的类型如何，都会保留适用于 double 输入的更严格容差。

参考文献：L.F. Shampine, "MATLAB program for quadrature in 2D", Applied Mathematics and Computation, Vol. 1, pp. 266–274, 2008.

另请参阅： quad2d、dblquad、integral、quad、quadgk、quadv、quadl、quadcc、trapz、integral3、triplequad。

q = integral3 (f, xa, xb, ya, yb, za, zb) ¶

q = integral3 (f, xa, xb, ya, yb, za, zb, prop, val, …) ¶

数值计算 f 的三重积分，在由 xa、xb、ya、yb、za、zb 定义的三维域上使用自适应求积（标量可以是有限值或无穷大）。此外，ya 和 yb 可以是 x 的标量函数，za 和 zb 可以是 x 和 y 的标量函数，从而允许在非矩形区域上进行积分。

f 是一个函数句柄、内联函数或包含要评估的函数名称的字符串。函数 f 必须具有形式 w = f(x,y,z)，并且所有运算必须向量化，使得 x、y 和 z 能接受相同大小的数组输入并返回数组输出。（可以假设 x、y 和 z 要么是相同大小的数组，要么是标量。）底层积分器会将积分点数组输入到 f 和/或使用内部向量展开来加速计算，但如果 f 不限于逐元素运算，则可能出错。对于被积函数无法避免这种情况时，下面描述的 ("Vectorized") 选项可能会产生更可靠的结果。

可以使用 "property", value 对指定额外的可选参数。有效属性包括：

AbsTol

定义求积的绝对误差容限。默认值为 1e-10（single 类型为 1e-5）。

RelTol

定义求积的相对误差容限。默认值为 1e-6（single 类型为 1e-4）。

Method

指定要使用的二维积分方法，有效值为 "auto"（默认）、"tiled" 或 "iterated"。使用 "auto" 时，Octave 将选择 "tiled" 方法，除非某个积分限为无穷大。

Vectorized

启用或禁用向量化积分。设置为 false 将强制 Octave 在调用被积函数时仅使用标量输入，这使得可以对尚未向量化或仅接受 x、y 或 z 标量值的 f(x,y,z) 进行积分。默认值为 true。请注意，这是通过使用 arrayfun 函数包装 f(x,y,z) 来实现的，这可能会显著降低计算速度。

自适应求积用于最小化误差估计，直到满足以下条件：

        error <= max (AbsTol, RelTol*|q|)

err 是积分误差 abs (q - I) 的近似上限，其中 I 是积分的精确值。

示例 1：在 x-y 平面上积分矩形区域

f = @(x,y) 2*ones (size (x));
q = integral3 (f, 0, 1, 0, 1, 0, 1)
  ⇒  q =  2

结果是一个体积，对于此常值被积函数，即为 Length * Width * Height。

示例 2：在球形体积上积分

f = @(x,y) ones (size (x));
ymax = @(x) sqrt (1 - x.^2);
zmax = @(x,y) sqrt (1 - x.^2 - y.^2);
q = integral3 (f, 0, 1, 0, ymax, 0, zmax)
  ⇒   q =  0.52360

对于此常值被积函数，结果即为体积，它是单位球体的 1/8，即 1/8 * 4/3 * pi。

示例 3：在立方体积上积分非向量化函数

f = @(x,y) sinc (x) * sinc (y), * sinc (z);
q = integral3 (f, -1, 1, -1, 1, -1, 1)
  ⇒  q =  14.535  (不正确)
q = integral3 (f, -1, 1, -1, 1, -1, 1, "Vectorized", false)
  ⇒  q =  1.6388 (正确)
f = @(x,y,z) sinc (x) .* sinc (y), .* sinc (z);
q = integral3 (f, -1, 1, -1, 1, -1, 1)
  ⇒  q =  1.6388  (正确)

第一个结果不正确，因为函数 f 中 sinc 函数之间的非逐元素运算符在 integral3 使用的内部积分数组之间产生了意外的矩阵乘法。在第二个结果中，将 "Vectorized" 设置为 false 强制 integral3 执行标量内部运算来计算积分，以增加约 30 倍的计算时间为代价得到了正确的数值结果。在第三个结果中，使用逐元素乘法运算符对被积函数 f 进行向量化，在不增加计算时间的情况下获得了正确的结果。

编程注意事项：如果积分区域内有奇点，最好将积分拆分，并将奇点置于边界上。

已知 MATLAB 不兼容性：如果未指定容差，且任何积分限为 single 类型，则 Octave 的积分函数会自动降低上述默认的绝对和相对误差容限。如果需要更严格的容差，则必须明确指定。MATLAB 则无论积分限的类型如何，都会保留适用于 double 输入的更严格容差。

参考文献：L.F. Shampine, "MATLAB program for quadrature in 2D", Applied Mathematics and Computation, Vol. 1, pp. 266–274, 2008.

另请参阅： triplequad、integral、quad、quadgk、quadv、quadl、quadcc、trapz、integral2、quad2d、dblquad。

上述积分可能相当慢，而且这一问题随积分维度的增加呈指数级增长。二维积分的另一种可能的解决方案是使用上一节中描述的正交配置（参见正交配置）。函数 f(x,y) 在 x 和 y 介于 0 和 1 之间的积分，可以使用 n 个点近似为 i=1:n 和 j=1:n 上 q(i)*q(j)*f(r(i),r(j)) 的总和，其中 q 和 r 由 colloc (n) 返回。对两个以上变量的推广是直接的。以下代码使用 n=8 个点计算了所讨论的积分。

f = @(x,y) sin (pi*x*y') .* sqrt (x*y');
n = 8;
[t, ~, ~, q] = colloc (n);
I = q'*f(t,t)*q;
      ⇒  0.30022

应该注意的是，点的数量决定了近似的质量。如果积分需要在 a 和 b 之间（而非 0 和 1 之间）进行，则需要进行变量替换。