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, nfun, 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 указывает на успешное интегрирование).

nfun указывает количество выполненных вычислений функции.

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, nfun] = 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.

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

Примечание: 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, nfun] = quadl (…)

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

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

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

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

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

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

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

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

Ссылка: 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 должна быть векторизована и возвращать вектор значений вывода при задании вектора входных значений.

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

quadgk (@(x) 1 ./ (sqrt (x) .* (x + 1)), 0, Inf)

Обратите внимание, что формулировка подынтегральной функции использует оператор поэлементного умножения ./ , и все пользовательские функции для quadgk должны делать то же самое.

Необязательный аргумент tol определяет абсолютную толерантность, используемую для остановки процедуры интегрирования. Значение по умолчанию — 1e-10 (1e-5 для одинарных).

Используемый алгоритм quadgk включает разбиение интервала интегрирования и вычисление каждого подынтервала. Если trace истинно, то после вычисления каждого из этих частичных интегралов отображаются: (1) количество подынтервалов на этом шаге, (2) текущая оценка ошибки err, (3) текущая оценка интеграла q.

Поведение алгоритма можно настроить, передав аргументы в quadgk как пары "prop", val. Допустимые свойства:

AbsTol

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

RelTol

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

MaxIntervalCount

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

WayPoints

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

quadgk (@(x) abs (1 - x.^2), 0, 2, "Waypoints", 1)

сигнализирует о точке разрыва в подынтегральной функции в x = 1.

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 — точное значение интеграла.

Ссылка: L.F. Shampine, «Векторизованная адаптивная квадратура в 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 для single), а относительная погрешность — 1e-6 (1e-4 для single).

Необязательный аргумент sing содержит список точек, в которых подынтегральное выражение имеет известные особенности или разрывы в любой из его производных внутри интервала интегрирования. Для примера выше, у которого есть разрыв в x=1, вызов quadcc будет следующим

int = quadcc (f, a, b, [], [ 1 ]);

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

err — оценка абсолютной погрешности интегрирования.

nr_points — количество точек, в которых было вычислено подынтегральное выражение.

Если адаптивное интегрирование не сошлось, значение 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, …)

Численно вычислить интеграл от f от a до b с помощью адаптивной квадратуры.

integral — это обёртка для quadcc (общие скалярные подынтегральные выражения), quadgk (интегралы с заданными путями интегрирования) и quadv (матричные подынтегральные выражения), предназначенная для обеспечения совместимости с MATLAB. Более тонкий контроль численного интегрирования может быть достигнут путём прямого вызова различных функций квадратуры.

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

a и b — нижний и верхний пределы интегрирования. Один или оба предела могут быть бесконечными или содержать слабые особенности на концах. Если один или оба предела являются комплексными, integral выполнит интеграл по прямой траектории. В качестве альтернативы, путь в комплексной области может быть задан с помощью опции "Waypoints" (см. ниже).

Дополнительные необязательные параметры можно задать с помощью пар "property", value. Допустимые свойства:

Waypoints

Указывает точки, которые будут использоваться для определения подынтервалов алгоритма квадратуры, или если a, b или waypoints являются комплексными, то квадратура вычисляется как контурный интеграл вдоль кусково-непрерывного пути. Более подробные сведения см. в quadgk.

ArrayValued

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

AbsTol

Определяет абсолютную погрешность толерантности для квадратуры. По умолчанию абсолютная погрешность равна 1e-10 (1e-5 для single).

RelTol

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

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

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

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

1. Если толерантности не указаны, а любые пределы интегрирования или контрольные точки являются типа single, то функции интегрирования Octave автоматически уменьшают значения абсолютной и относительной погрешности по умолчанию, как указано выше. Если требуются более жёсткие толерантности, они должны быть указаны. MATLAB оставляет более жёсткие толерантности, соответствующие входам double на месте независимо от класса пределов интегрирования.
2. Вследствие использования quadcc, quadgk, и quadv, определённые комбинации опций не поддерживаются. В настоящее время "ArrayValued" не может быть объединён с "RelTol" или "Waypoints".

См. также: 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.