23.1 Функции одной переменной

Octave поддерживает пять различных адаптивных алгоритмов квадратуры для вычисления интеграла функции f на интервале от a до b. Это:

quad

Численное интегрирование на основе квадратуры Гаусса.

quadv

Численное интегрирование с использованием адаптивного векторизованного правила Симпсона.

quadl

Численное интегрирование с использованием адаптивного правила Лобатто.

quadgk

Численное интегрирование с использованием адаптивного правила Гаусса — Кронрода.

quadcc

Численное интегрирование с использованием адаптивных правил Кленшоу — Кёртиса.

Кроме того, доступны следующие функции:

integral

Функция-обёртка для совместимости, которая выбирает между quadv и quadgk в зависимости от подынтегральной функции и выбранных параметров.

trapz, cumtrapz

Численное интегрирование данных методом трапеций.

Выбор наиболее подходящего алгоритма квадратуры зависит от подынтегральной функции. Если у вас есть эмпирические данные, а не функция, выбирайте trapz или cumtrapz. Если характеристики подынтегральной функции неизвестны, наиболее надёжным вариантом будет quadcc, поскольку он может обрабатывать разрывы, сингулярности, осциллирующие функции и бесконечные интервалы. Для гладкой подынтегральной функции самым быстрым алгоритмом может быть quadgk.

Ниже приведён пример использования 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 с использованием процедур Fortran из QUADPACK.

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.

: quad_options ()

: val = quad_options (opt)

: quad_options (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 с использованием адаптивного правила Симпсона.

f — это дескриптор функции, встроенная функция или строка, содержащая имя вычисляемой функции. quadv — векторизованная версия quad, и функция, заданная через f, должна принимать на вход скаляр или вектор и возвращать скаляр, вектор или массив.

a и b — нижний и верхний пределы интегрирования. Оба предела должны быть конечными.

Необязательный аргумент tol задаёт абсолютную погрешность, при достижении которой процедура адаптации останавливается. Значение по умолчанию — 1e-6.

Используемый алгоритм quadv рекурсивно делит интервал интегрирования на подынтервалы и применяет правило Симпсона к каждому из них. Если trace имеет значение true, после вычисления каждого из частичных интегралов отображаются: (1) общее число вычислений функции, (2) левый конец подынтервала, (3) длина подынтервала, (4) приближённое значение интеграла на подынтервале.

Дополнительные аргументы p1 и т. д. передаются непосредственно функции f. Чтобы использовать значения tol и trace по умолчанию, можно передать пустые матрицы ([]).

Результат интегрирования возвращается в q.

Необязательный выходной аргумент nfev указывает общее количество выполненных вычислений функции.

Примечание: quadv написана на языке сценариев Octave и, в отличие от функции quad, может рекурсивно использоваться в dblquad и triplequad.

См. также: 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 с использованием адаптивного правила Лобатто.

f — это дескриптор функции, встроенная функция или строка, содержащая имя вычисляемой функции. Функция f должна быть векторизованной и возвращать вектор выходных значений при получении вектора входных значений.

a и b — нижний и верхний пределы интегрирования. Оба предела должны быть конечными.

Необязательный аргумент tol задаёт абсолютную погрешность, с которой выполняется интегрирование. Значение по умолчанию — 1e-6.

Используемый алгоритм quadl рекурсивно делит интервал интегрирования. Если задан trace, для каждого подынтервала отображаются: (1) общее число вычислений функции, (2) левый конец подынтервала, (3) длина подынтервала, (4) приближённое значение интеграла на подынтервале.

Дополнительные аргументы p1 и т. д. передаются непосредственно функции f. Чтобы использовать значения tol и trace по умолчанию, можно передать пустые матрицы ([]).

Результат интегрирования возвращается в q.

Необязательный выходной аргумент nfev указывает общее количество выполненных вычислений функции.

Источник: W. Gander и W. Gautschi, Адаптивная квадратура — повторный обзор, BIT, т. 40, № 1, март 2000 г., с. 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.

Если какая-либо из промежуточных точек комплексная, выполняется интегрирование по контуру, как описано ниже.

MaxIntervalCount

quadgk изначально делит интервал, на котором выполняется квадратура, на 10 интервалов либо — если заданы промежуточные точки — по каждой промежуточной точке. Подынтервалы с недопустимо большой погрешностью делятся и вычисляются повторно. Если в какой-либо момент число подынтервалов превышает 650, сообщается о плохой сходимости и возвращается текущая оценка интеграла. Свойство "MaxIntervalCount" позволяет изменить максимальное число подынтервалов, после достижения которого вычисления прекращаются.

Trace

Если имеет логическое значение true, 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])

вычисляет интеграл log (z) вдоль квадрата, заданного через [1+1i, -1+1i, -1-1i, +1-1i].

Результат интегрирования возвращается в q.

err — это приближённая граница погрешности интеграла abs (q - I), где I — точное значение интеграла. Если адаптивное интегрирование не сошлось, значение err будет больше заданной погрешности. Если запрошен только один выходной аргумент, при недостижении требуемой точности будет выдано предупреждение. Если запрошен второй выходной аргумент err, предупреждение не выдаётся; программист должен самостоятельно проверить, удовлетворительны ли результаты.

Источник: L.F. Shampine, «Векторизованная адаптивная квадратура в MATLAB», Journal of Computational and Applied Mathematics, с. 131–140, т. 211, вып. 2, февраль 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 с использованием двухступенчатой квадратуры Кленшоу — Кёртиса.

f — это дескриптор функции, встроенная функция или строка, содержащая имя вычисляемой функции. Функция f должна быть векторизованной и возвращать вектор выходных значений при получении вектора входных значений. Например:

f = @(x) x .* sin (1./x) .* sqrt (abs (1 - x));

в этом примере для всех операторов используется поэлементная «точечная» форма.

a и b — нижний и верхний пределы интегрирования. Один или оба предела могут быть бесконечными. quadcc обрабатывает бесконечный предел, заменяя переменную интегрирования на x = tan (pi/2*u).

Необязательный аргумент tol — это вектор из одного или двух элементов, задающий требуемую точность результата. Первый элемент вектора — требуемая абсолютная погрешность, второй — требуемая относительная погрешность. Чтобы проверять только относительную погрешность, задайте абсолютную погрешность равной нулю. Чтобы проверять только абсолютную погрешность, задайте относительную погрешность равной нулю. Абсолютная погрешность по умолчанию равна 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.

Алгоритм использует квадратурные правила Кленшоу — Кёртиса возрастающей степени на каждом интервале и делит интервал пополам, если функция не выглядит гладкой или достигнута максимальная степень правила. Оценка погрешности вычисляется по L2-норме разности между двумя последовательными интерполяциями подынтегральной функции в узлах соответствующих квадратурных правил.

Источник: P. Gonnet, Повышение надёжности адаптивной квадратуры с помощью явных интерполянтов, ACM Transactions on Mathematical Software, т. 37, вып. 3, статья № 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.

ArrayValued

integral предполагает, что f возвращает скалярное значение, если параметру arrayvalued не присвоено значение true. Этот параметр заставляет integral выполнять интегрирование по всему массиву и возвращать q с теми же размерами, что и результат функции f. Подробнее см. quadgk.

AbsTol

Задает абсолютный допуск ошибки для квадратуры. Абсолютный допуск по умолчанию равен 1e-10 (1e-5 для одинарной точности).

RelTol

Задает относительный допуск ошибки для квадратуры. Относительный допуск по умолчанию равен 1e-6 (1e-4 для одинарной точности).

Необязательный выходной параметр err содержит оценку абсолютной ошибки, использованную вызванным интегратором.

Адаптивная квадратура используется для минимизации оценки ошибки, пока не будет выполнено следующее условие:

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

Известные несовместимости с MATLAB:

1. Если допуски не заданы и какой-либо предел интегрирования или промежуточная точка имеет тип single, функции интегрирования Octave автоматически уменьшают абсолютный и относительный допуски по умолчанию, как указано выше. Если требуются более строгие допуски, их необходимо задать явно. MATLAB сохраняет более строгие допуски, подходящие для входных данных типа double, независимо от класса пределов интегрирования.

См. также: integral2, integral3, quad, quadgk, quadv, quadl, quadcc, trapz, dblquad, triplequad.

Иногда функция недоступна, и имеются только исходные точки (x, y), по которым нужно выполнить интегрирование. Такое может произойти при сборе данных в эксперименте. Функция trapz может интегрировать эти значения, как показано в следующем примере, где данные "data" собраны для функции косинуса на интервале [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).

См. также: trapz, cumsum.