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

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

quad

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

quadv

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

quadl

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

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

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

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

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

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

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

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

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

Ссылка: W. Gander и 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

Указывает точки, которые станут конечными точками подынтервалов в алгоритме, что может привести к значительному улучшению оценки ошибки интеграла, более быстрому вычислению или тому и другому. Это может быть полезно для указания большего количества подынтервалов вокруг области, где подынтегральная функция быстро меняется, или для обозначения мест, где имеется разрыв в первой производной функции. Например, функция signum имеет разрыв в x == 0 и путем указания контрольной точки

quadgk (@(x) sign (x), -0.5, 1, "Waypoints", [0])

предел ошибки уменьшается с 4e-7 до 1e-13.

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

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

MaxIntervalCount

quadgk первоначально делит интервал, на котором следует выполнить квадратуру, на 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])

интегрирует log (z) по квадрату, определяемому [1+1i, -1+1i, -1-1i, +1-1i].

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

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

Ссылка: Л.Ф. Шэмп, "Векторизованная адаптивная квадратура в MATLAB", Журнал вычислительной и прикладной математики, стр. 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 — это вектор из 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.

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

Ссылка: П. Гоннет, Повышение надежности адаптивной квадратуры с использованием явных интерполянтов, 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 выполнит интегрирование по прямой.

Дополнительные необязательные параметры могут быть заданы с помощью пар "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.