Octave 支持五种不同的自适应求积算法,用于计算函数 f 在区间 a 到 b 上的积分。这些算法是:
quad基于高斯求积的数值积分。
quadv使用自适应向量化 Simpson 规则的数值积分。
quadl使用自适应 Lobatto 规则的数值积分。
quadgk使用自适应高斯-科恩罗德规则的数值积分。
quadcc使用自适应克伦肖-柯蒂斯规则的数值积分。
此外,还提供了以下函数:
integral一个兼容性包装函数,将根据被积函数和所选选项在 quadv 和 quadgk 之间进行选择。
trapz, cumtrapz使用梯形方法对数据进行数值积分。
最佳求积算法的选择取决于被积函数。如果你有经验数据而非函数,则应选择 trapz 或 cumtrapz。如果你不确定被积函数的特性,quadcc 是最稳健的选择,因为它能处理不连续性、奇点、振荡函数和无限区间。当被积函数平滑时,quadgk 可能是最快的算法。
| 函数 | 特性 | |
|---|---|---|
| quad | 非平滑被积函数精度低 | |
| quadv | 平滑被积函数中等精度 | |
| quadl | 平滑被积函数中等精度。比 quadgk 慢。 | |
| quadgk | 平滑被积函数中等精度(1e-6 – 1e-9)。 | |
| 处理振荡函数和无限区间 | ||
| quadcc | 非平滑/平滑被积函数低到高精度 | |
| 处理振荡函数、奇点和无限区间 |
以下是使用 quad 来积分函数的示例:
f(x) = x * sin (1/x) * sqrt (abs (1 - x))
从 x = 0 到 x = 3。
这是一个相当困难的积分(绘制积分范围内的函数以了解原因)。
第一步是定义函数:
function y = f (x) y = x .* sin (1./x) .* sqrt (abs (1 - x)); endfunction
注意运算符的"点"形式的使用。这对于 quad 积分器来说不是必需的,但其他积分器需要。在任何情况下,它都可以更容易地生成一组用于绘制的点,因为可以使用向量参数调用函数来生成向量结果。
第二步是使用积分限调用 quad:
[q, ier, nfun, err] = quad ("f", 0, 3)
⇒ 1.9819
⇒ 1
⇒ 5061
⇒ 1.1522e-07
虽然 quad 返回了 ier 的非零值,但结果相当准确(要了解原因,请检查如果将下限移动到 0.1、0.01、0.001 等时结果会发生什么变化)。
函数 "f" 可以是函数的字符串名称或函数句柄。这些选项使得无需在 m 文件中完全定义函数即可轻松进行积分。例如:
# Verify gamma function = (n-1)! for n = 4
f = @(x) x.^3 .* exp (-x);
quadcc (f, 0, Inf)
⇒ 6.0000
q = quad (f, a, b) ¶q = quad (f, a, b, tol) ¶q = quad (f, a, b, tol, sing) ¶[q, ier, nfev, err] = quad (…) ¶数值计算函数 f 从 a 到 b 的积分,使用 QUADPACK 的 Fortran 例程。
f 是一个函数句柄、内联函数或包含要计算的函数名称的字符串。函数必须具有形式 y = f (x),其中 y 和 x 是标量。
a 和 b 是积分的下限和上限。两者都可以是无穷大。
可选参数 tol 是一个向量,用于指定结果的期望精度。向量的第一个元素是所需的绝对容差,第二个元素是所需的相对容差。要仅选择相对测试,请将绝对容差设置为零。要仅选择绝对测试,请将相对容差设置为零。两个容差默认值为 sqrt (eps),约为 1.5e-8。
可选参数 sing 是一个向量,包含被积函数已知存在奇点的值。
积分的结果在 q 中返回。
ier 包含一个整数错误代码(0 表示积分成功)。
nfev 表示已执行的函数评估次数。
err 包含对解的误差估计。
函数 quad_options 可以为 quad 设置其他可选参数。
注意:由于 quad 是用 Fortran 编写的,它不能递归调用。这阻止了它在通过 dblquad 和 triplequad 例程对多个变量进行积分时使用。
另请参阅: quad_options, quadv, quadl, quadgk, quadcc, trapz, dblquad, triplequad.
() ¶val = quad_options (opt) ¶(opt, val) ¶查询或设置函数 quad 的选项。
当在没有参数的情况下调用时,将显示所有可用选项的名称及其当前值。
给定一个参数,返回选项 opt 的值。
当用两个参数调用时,quad_options 将选项 opt 设置为值 val。
选项包括
"absolute tolerance"绝对容差;对于纯相对误差测试可以为零。
"relative tolerance"非负相对容差。如果绝对容差为零,则相对容差必须大于或等于 max (50*eps, 0.5e-28)。
"single precision absolute tolerance"单精度的绝对容差;对于纯相对误差测试可以为零。
"single precision relative tolerance"单精度的非负相对容差。如果绝对容差为零,则相对容差必须大于或等于 max (50*eps, 0.5e-28)。
q = quadv (f, a, b) ¶q = quadv (f, a, b, tol) ¶q = quadv (f, a, b, tol, trace) ¶q = quadv (f, a, b, tol, trace, p1, p2, …) ¶[q, nfev] = quadv (…) ¶数值计算函数 f 从 a 到 b 的积分,使用自适应 Simpson 规则。
f 是一个函数句柄、内联函数或字符串,其中包含要计算的函数的名称。quadv 是 quad 的向量化版本,由 f 定义的函数必须接受标量或向量作为输入,并返回标量、向量或数组作为输出。
a 和 b 是积分的下限和上限。两个极限都必须是有限的。
可选参数 tol 定义了用于停止自适应过程的绝对容差。默认值为 1e-6。
quadv 使用的算法涉及递归细分积分区间,并在每个子区间上应用 Simpson 规则。如果 trace 为 true,则在计算这些部分积分中的每一个之后显示:(1) 函数评估的总数,(2) 子区间的左端,(3) 子区间长度,(4) 子区间上积分的近似值。
其他参数 p1 等直接传递给函数 f。要使用 tol 和 trace 的默认值,可以传递空矩阵 ([])。
积分的结果在 q 中返回。
可选输出 nfev 表示执行的函数评估的总数。
注意:quadv 是用 Octave 的脚本语言编写的,可以在 dblquad 和 triplequad 中递归使用,这与 quad 函数不同。
另请参阅: quad, quadl, quadgk, quadcc, trapz, dblquad, triplequad, integral, integral2, integral3.
q = quadl (f, a, b) ¶q = quadl (f, a, b, tol) ¶q = quadl (f, a, b, tol, trace) ¶q = quadl (f, a, b, tol, trace, p1, p2, …) ¶[q, nfev] = quadl (…) ¶数值计算函数 f 从 a 到 b 的积分,使用自适应 Lobatto 规则。
f 是一个函数句柄、内联函数或字符串,其中包含要计算的函数的名称。函数 f 必须向量化,并在给定输入值向量时返回输出值向量。
a 和 b 是积分的下限和上限。两个极限都必须是有限的。
可选参数 tol 定义执行积分的绝对容差。默认值为 1e-6。
quadl 使用的算法涉及递归细分积分区间。如果定义了 trace,则对每个子区间显示:(1) 函数评估的总数,(2) 子区间的左端,(3) 子区间长度,(4) 子区间上积分的近似值。
其他参数 p1 等直接传递给函数 f。要使用 tol 和 trace 的默认值,可以传递空矩阵 ([])。
积分的结果在 q 中返回。
可选输出 nfev 表示执行的函数评估的总数。
参考文献:W. Gander and W. Gautschi, Adaptive Quadrature - Revisited, BIT Vol. 40, No. 1, March 2000, pp. 84–101. https://www.inf.ethz.ch/personal/gander/
另请参阅: quad, quadv, quadgk, quadcc, trapz, dblquad, triplequad, integral, integral2, integral3.
q = quadgk (f, a, b) ¶q = quadgk (f, a, b, abstol) ¶q = quadgk (f, a, b, abstol, trace) ¶q = quadgk (f, a, b, "prop", val, …) ¶[q, err] = quadgk (…) ¶数值计算函数 f 从 a 到 b 的积分,使用自适应高斯-科恩罗德求积。
f 是一个函数句柄、内联函数或字符串,其中包含要计算的函数的名称。函数 f 必须向量化,并在给定输入值向量时返回输出值向量(有关此规则的例外,请参见属性 "ArrayValued")。
a 和 b 是积分的下限和上限。其中一个或两个极限可以是无穷大或包含弱端点奇点。变量变换将用于处理任何无限区间并削弱奇点。例如:
quadgk (@(x) 1 ./ (sqrt (x) .* (x + 1)), 0, Inf)
注意,被积函数的公式使用逐元素运算符 ./,所有传递给 quadgk 的用户函数也应如此。
可选参数 abstol 定义用于停止积分过程的绝对容差。默认值为 1e-10(单精度为 1e-5)。
quadgk 使用的算法涉及细分积分区间并评估每个子区间。如果 trace 为 true,则在计算每个部分积分后显示:(1) 此步骤的子区间数,(2) 误差 err 的当前估计,(3) 积分 q 的当前估计。
算法的行为可以通过向 quadgk 传递参数对 "prop", val 来配置。有效属性为:
AbsTol定义求积的绝对误差容限。默认绝对容差为 1e-10(单精度为 1e-5)。
RelTol定义求积的相对误差容限。默认相对容差为 1e-6(单精度为 1e-4)。
ArrayValued当设置为 true 时,函数 f 为标量输入生成数组输出。默认值为 false,这要求 f 生成与输入大小相同的输出。例如:
quadgk (@(x) x .^ (1:5), 0, 2, "ArrayValued", 1)
将在一个函数调用中积分 [x.^1, x.^2, x.^3, x.^4, x.^5],而不必重复定义单个匿名函数并正常调用 quadgk。
WayPoints指定将成为算法中子区间端点的点,这可以显著改进积分误差的估计、加快计算速度,或两者兼而有之。在被积函数快速变化的区域周围指定更多的子区间,或者在函数的一阶导数中存在不连续性的标志位置,这可能很有用。例如,符号函数在 x == 0 处有不连续性,通过指定一个路径点:
quadgk (@(x) sign (x), -0.5, 1, "Waypoints", [0])
误差范围从 4e-7 减小到 1e-13。
如果函数在积分区域内具有奇点,则不应通过路径点来处理。相反,应将整体积分分解为几个较小积分的总和,使得奇点出现在调用 quadgk 时的一个积分边界上。
如果任何路径点是复数,则执行围道积分,如下所述。
MaxIntervalCountquadgk 最初将执行求积的区间细分为 10 个区间,或者如果给出了路径点,则在每个路径点处细分。具有不可接受误差的子区间被细分并重新评估。如果子区间的数量在任何时候超过 650 个子区间,则发出收敛性差的信号,并返回积分的当前估计。属性 "MaxIntervalCount" 可用于更改退出前允许存在的子区间数。
Trace如果为逻辑真,quadgk 在每次迭代时打印关于求积收敛性的信息。
如果 a、b 或 waypoints 中的任何一个是复数,则求积被视为沿由 [a, waypoints(1), waypoints(2), …, b] 定义的分段线性路径的围道积分。在这种情况下,假设积分没有边缘奇点。例如:
quadgk (@(z) log (z), 1+1i, 1+1i, "WayPoints",
[-1+1i, -1-1i, +1-1i])
沿由 [1+1i, -1+1i, -1-1i, +1-1i] 定义的正方形积分 log (z)。
积分的结果在 q 中返回。
err 是积分误差的近似界 abs (q - I),其中 I 是积分的精确值。如果自适应积分没有收敛,则 err 的值将大于所请求的容差。如果只请求一个输出,则在未满足请求的容差时会发出警告。如果请求了第二个输出 err,则不发出警告,由程序员负责检查并确定结果是否令人满意。
参考文献:L.F. Shampine, "Vectorized adaptive quadrature in MATLAB", Journal of Computational and Applied Mathematics, Vol. 211, Issue 2, pp. 131–140, Feb 2008.
另请参阅: quad, quadv, quadl, quadcc, trapz, dblquad, triplequad, integral, integral2, integral3.
q = quadcc (f, a, b) ¶q = quadcc (f, a, b, tol) ¶q = quadcc (f, a, b, tol, sing) ¶[q, err, nr_points] = quadcc (…) ¶数值计算函数 f 从 a 到 b 的积分,使用双自适应 Clenshaw-Curtis 求积。
f 是一个函数句柄、内联函数或字符串,其中包含要计算的函数的名称。函数 f 必须向量化,并且在给定输入值向量时必须返回输出值向量。例如:
f = @(x) x .* sin (1./x) .* sqrt (abs (1 - x));
它对所有运算符使用逐元素的"点"形式。
a 和 b 是积分的下限和上限。其中一个或两个极限都可以是无穷大。quadcc 通过将积分变量替换为 x = tan (pi/2*u) 来处理无限极限。
可选参数 tol 是一个 1 元或 2 元向量,用于指定所需的结果精度。向量的第一个元素是所需的绝对容差,第二个元素是所需的相对容差。要仅选择相对测试,请将绝对容差设置为零。要仅选择绝对测试,请将相对容差设置为零。默认绝对容差为 1e-10(单精度为 1e-5),默认相对容差为 1e-6(单精度为 1e-4)。
可选参数 sing 包含被积函数在积分区间内具有已知奇点或其任何导数中存在不连续性的点的列表。对于上面在 x=1 处有不连续性的示例,调用 quadcc 的方式如下:
int = quadcc (f, a, b, [], [ 1 ]);
积分的结果在 q 中返回。
err 是绝对积分误差的估计值。
nr_points 是评估被积函数的点数。
如果自适应积分没有收敛,则 err 的值将大于所请求的容差。如果只请求一个输出,则在未满足所请求的容差时会发出警告。如果请求了第二个输出 err,则不发出警告,由程序员负责检查并确定结果是否令人满意。
quadcc 能够处理被积函数的非数值,例如 NaN 或 Inf。如果积分发散,并且 quadcc 检测到这种情况,则会发出警告,并返回 Inf 或 -Inf。
注意:quadcc 是一种通用求积算法,因此对于平滑或其他表现良好的被积函数,它可能不如其他方法(如 quadgk)高效。
该算法在每个区间中使用次数递增的 Clenshaw-Curtis 求积规则,如果函数看起来不光滑或已达到最大次数规则,则将区间一分为二。误差估计是根据被积函数在相应求积规则节点上的两个连续插值之间的差值的 L2 范数计算的。
参考文献:P. Gonnet, "Increasing the Reliability of Adaptive Quadrature Using Explicit Interpolants", ACM Transactions on Mathematical Software, Vol. 37, Issue 3, Article No. 3, 2010.
另请参阅: quad, quadv, quadl, quadgk, trapz, dblquad, triplequad.
q = integral (f, a, b) ¶q = integral (f, a, b, prop, val, …) ¶[q, err] = integral (…) ¶数值计算函数 f 从 a 到 b 的积分,使用自适应求积。
integral 是 quadcc(一般实值标量被积函数和极限)和 quadgk(具有指定积分路径和数组值被积函数的积分)的包装器,旨在提供 MATLAB 兼容性。通过直接调用各种求积函数,可以对数值积分进行更多控制。
f 是一个函数句柄、内联函数或字符串,其中包含要计算的函数的名称。函数 f 必须向量化,并在给定输入值向量时返回输出值向量。
a 和 b 是积分的下限和上限。其中一个或两个极限可以是无穷大或包含弱端点奇点。如果其中一个或两个极限是复数,integral 将执行直线路径积分。或者,可以使用 "Waypoints" 选项指定复数域路径(见下文)。
可以使用 "property", value 对指定其他可选参数。有效属性为:
Waypoints指定用于定义求积算法子区间的点,或者如果 a、b 或 waypoints 是复数,则求积将作为沿分段连续路径的围道积分进行计算。更多详细信息,请参见 quadgk。
ArrayValuedintegral 期望 f 返回标量值,除非 arrayvalued 被指定为 true。此选项将导致 integral 在整个数组上执行积分,并返回与 f 返回的维度相同的 q。更多详细信息,请参见 quadgk。
AbsTol定义求积的绝对误差容限。默认绝对容差为 1e-10(单精度为 1e-5)。
RelTol定义求积的相对误差容限。默认相对容差为 1e-6(单精度为 1e-4)。
可选输出 err 包含被调用的积分器使用的绝对误差估计。
使用自适应求积来最小化误差估计,直到满足以下条件:
error <= max (AbsTol, RelTol*|q|).
已知的 MATLAB 不兼容性:
single 类型,则 Octave 的积分函数会自动降低上面指定的默认绝对和相对误差容限。如果需要更严格的容差,则必须指定它们。MATLAB 无论积分极限的类别如何,都会保持适用于 double 输入的更严格容差。另请参阅: integral2, integral3, quad, quadgk, quadv, quadl, quadcc, trapz, dblquad, triplequad.
有时我们没有函数,而只有原始的 (x, y) 点来进行积分。这在实验中收集数据时可能会出现。trapz 函数可以对这些值进行积分,如下面的示例所示,其中"数据"是在范围 [0, pi/2) 上的余弦函数上收集的。
x = 0:0.1:pi/2; # Uniformly spaced points
y = cos (x);
trapz (x, y)
⇒ 0.99666
答案相当接近精确值 1。普通求积对被积函数的特性很敏感。经验积分不仅取决于被积函数,还取决于为表示函数而选择的特定点。在范围 [0, pi/2) 上使用正弦函数重复上面的示例会产生差得多的结果。
x = 0:0.1:pi/2; # Uniformly spaced points
y = sin (x);
trapz (x, y)
⇒ 0.92849
然而,数据点的选择略有不同会显著改变结果。相同的积分,相同的点数,但间隔不同,可以得到更精确的答案。
x = linspace (0, pi/2, 16); # Uniformly spaced, but including endpoint
y = sin (x);
trapz (x, y)
⇒ 0.99909
一般来说,可能无法预先知道点的最佳分布。或者这些点可能来自一个实验,其中没有选择最佳分布的自由。在任何情况下,使用 trapz 时都必须意识到这个问题。
q = trapz (y) ¶q = trapz (x, y) ¶q = trapz (…, dim) ¶数值计算点 y 的积分,使用梯形方法。
trapz (y) 计算 y 沿第一个非奇异维度的积分。当省略参数 x 时,假定为等间距的 x 向量,间距为 1。trapz (x, y) 根据 x 的间距和 y 的值计算积分。如果 y 中的点采样不均匀,这很有用。
如果给出了可选参数 dim,则沿此维度操作。
应用说明:如果未指定 x,则将使用单位间距。要将积分缩放到正确的值,必须乘以实际间距值 (deltaX)。例如,x^3 在范围 [0, 1] 上的积分是 x^4/4,即 0.25。以下代码使用 trapz 以三种不同的方式计算该积分。
x = 0:0.1:1; y = x.^3; ## No scaling q = trapz (y) ⇒ q = 2.5250 ## Approximation to integral by scaling q * 0.1 ⇒ 0.25250 ## Same result by specifying x trapz (x, y) ⇒ 0.25250
另请参阅: cumtrapz.
q = cumtrapz (y) ¶q = cumtrapz (x, y) ¶q = cumtrapz (…, dim) ¶点 y 的累积数值积分,使用梯形方法。
cumtrapz (y) 计算 y 沿第一个非奇异维度的累积积分。其中 trapz 仅报告整体积分和,而 cumtrapz 报告 y 在每个点处的当前部分和。
当省略参数 x 时,假定为等间距的 x 向量,间距为 1。cumtrapz (x, y) 根据 x 的间距和 y 的值计算积分。如果 y 中的点采样不均匀,这很有用。
如果给出了可选参数 dim,则沿此维度操作。
应用说明:如果未指定 x,则将使用单位间距。要将积分缩放到正确的值,必须乘以实际间距值 (deltaX)。
版权所有 © 2024-2026 Octave中文网
ICP备案/许可证号:黑ICP备2024030411号-4