23.3 Функции от нескольких переменных

Octave включает несколько функций для вычисления интеграла функций от нескольких переменных. Этот процесс, как правило, выполняется путем создания функции, которая интегрирует f по x, а затем интегрирует полученную функцию по y. Этот процесс можно выполнить вручную, используя следующий пример, который интегрирует функцию:

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

для x и y в пределах от 0 до 1.

Используя quadgk в примере ниже, можно выполнить двойное интегрирование. (Обратите внимание, что любой из функций квадратур 1D может быть использован таким образом, за исключением 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.

Рекурсивный алгоритм квадратур, представленный выше, называется "iterated". Отдельный метод двумерного интегрирования реализован в функции quad2d. Эта функция выполняет "tiled" интегрирование, разбивая область интегрирования на прямоугольные области и выполняя отдельные интегрирования по этим областям. Области далее разбиваются в областях, требующих уточнения, чтобы достичь желаемой точности. Для некоторых функций этот метод может быть быстрее, чем двумерная итерация, используемая в других функциях выше.

: 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.

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

AbsTol

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

RelTol

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

MaxFunEvals

Максимальное число вызовов функции векторизованной функции f. Значение по умолчанию равно 5000.

Singular

Включить/отключить преобразования для ослабления особенностей на границе области интегрирования. Значение по умолчанию равно true.

Vectorized

Вариант для отключения векторизованного интегрирования, принуждая Octave использовать только скалярные входные данные при вызове подынтегральной функции. Значение по умолчанию равно false.

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

Результат — объём, который для этой подынтегральной функции с постоянным значением равен Площадь треугольника × Высота или 1/2 * Base * Width * Height.

Примечания по программированию: Если внутри области интегрирования есть особенности, лучше разбить интеграл и поместить особенности на границу.

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

Ссылка: Л.Ф. Шэмпин, MATLAB программа для квадратуры в 2D, Прикладная математика и вычисления, стр. 266–274, Том 1, 2008.

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

Наконец, функции integral2 и integral3 предоставляются в качестве общих функций для интегрирования в 2D и 3D. Они автоматически выбирают между методами интегрирования по частям и плитками и, в отличие от 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.

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

AbsTol

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

RelTol

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

Method

Указать метод двумерного интегрирования, который будет использоваться, допустимые варианты — "auto" (по умолчанию), "tiled", или "iterated". При использовании "auto", Octave выберет метод "tiled", если какой-либо из пределов интегрирования не является конечным.

Vectorized

Включить или отключить векторизованное интегрирование. Значение false принуждает Octave использовать только скалярные входные данные при вызове подынтегральной функции, что позволяет использовать подынтегральные функции f(x,y), которые не были векторизованы и принимают только x и y в качестве скаляров. Значение по умолчанию — true.

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

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

Результат — объём, который для этой подынтегральной функции с постоянным значением равен Площадь треугольника × Высота или 1/2 * Base * Width * Height.

Примечания по программированию: Если внутри области интегрирования есть особенности, лучше разбить интеграл и поместить особенности на границу.

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

Ссылка: Л.Ф. Шэмпин, MATLAB программа для квадратуры в 2D, Прикладная математика и вычисления, стр. 266–274, Том 1, 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 должна иметь вид z = f(x,y), где x — вектор, а y — скаляр. Она должна возвращать вектор той же длины и ориентации, что и x.

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

AbsTol

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

RelTol

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

Method

Указать метод двумерного интегрирования, допустимые варианты — "auto" (по умолчанию), "tiled", или "iterated". При использовании "auto", Octave выберет метод "tiled", если какие-либо пределы интегрирования не являются бесконечными.

Vectorized

Включить или выключить векторизованное интегрирование. Значение false заставляет Octave использовать только скалярные входные данные при вызове интегрируемой функции, что позволяет использовать интегрируемые функции f(x,y), которые не были векторизованы и принимают только x и y как скаляры. Значение по умолчанию — true.

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

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

err является приблизительной границей погрешности в интеграле abs (q - I), где I — точное значение интеграла.

Пример 1: интегрирование по прямоугольному объёму

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

Для этой интегрируемой функции с постоянным значением результат представляет собой объём, который равен просто 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.

Примечания по программированию: если в области интегрирования есть особенности, лучше разбить интеграл и поместить особенности на границу.

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

Ссылка: Л.Ф. Шампин, MATLAB программа для квадратур в 2D, Прикладная математика и вычисления, стр. 266–274, Том 1, 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, то требуется изменение переменных.