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

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

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

для x и y между 0 и 1.

Используя quadgk в примере ниже, можно выполнить двойное интегрирование. (Обратите внимание, что любая из функций квадратуры 1-D может быть использована таким образом, за исключением 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 для интегралов по двум переменным. 3-мерный эквивалент этого процесса реализован в 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". Отдельный метод интегрирования 2-D реализован в функции quad2d. Эта функция выполняет "tiled" интегрирование, разбивая область интегрирования на прямоугольные области и выполняя отдельные интегрирования по этим областям. Области далее разбиваются в областях, требующих уточнения, чтобы достичь желаемой численной точности. Для некоторых функций этот метод может быть быстрее, чем итерация 2-D, используемая в других функциях выше.

: 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 и y будут либо массивами одинакового размера, либо один будет скаляром.) Подлежащие интеграторы будут вводить массивы точек интегрирования в f и/или использовать внутренние векторные расширения для ускорения вычислений, что может привести к непредсказуемым результатам, если f не ограничена элементарными операциями. Для интегралов, где это неизбежно, ("Vectorized") option described below may produce более надёжные результаты.

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

AbsTol

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

RelTol

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

MaxFunEvals

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

Singular

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

Vectorized

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

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

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

Пример 3 : интегрирование не-векторизованной функции по прямоугольной области

f = @(x,y) sinc (x) * sinc (y));
q = quad2d (f, -1, 1, -1, 1)
  ⇒ q =  12.328  (incorrect)
q = quad2d (f, -1, 1, -1, 1, "Vectorized", false)
  ⇒ q =  1.390 (correct)
f = @(x,y) sinc (x) .* sinc (y);
q = quad2d (f, -1, 1, -1, 1)
  ⇒ q =  1.390  (correct)

Первый результат неверен, так как неэлементная операция между функциями sinc в f создает непреднамеренные матричные умножения между внутренними массивами интегрирования, используемыми quad2d. Во втором результате установка "Vectorized" в false заставляет quad2d выполнять скалярные внутренние операции для вычисления интеграла, что приводит к правильному численному результату за счёт примерно 20-кратного увеличения времени вычислений. В третьем результате векторизация интегрируемой функции f с использованием элементарного оператора умножения даёт правильный результат без увеличения времени вычислений.

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

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

Ссылка: L.F. Shampine, MATLAB программа для квадратуры в 2D, Applied Mathematics and Computation, стр. 266–274, том 1, 2008.

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

Наконец, функции integral2 и integral3 предоставляются как общие функции 2-D и 3-D интегрирования. Они автоматически выбирают между методами интегрирования по итерациям и по плиткам и, в отличие от 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 и y будут либо массивами одного размера, либо один будет скаляром.) Базовые интегрирующие модули будут вводить массивы интеграционных точек в f и/или использовать внутренние векторные расширения для ускорения вычислений, что может привести к непредсказуемым результатам, если f не ограничена элементарными операциями. Для интегральных слагаемых, где это неизбежно, ("Vectorized") option described below may produce более надёжные результаты.

Дополнительные необязательные параметры можно указать с помощью "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. Обратите внимание, что это достигается путём обертывания f(x,y) функцией arrayfun, что может значительно снизить скорость вычислений.

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

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

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

Пример 3: интегрирование не векторизованной функции по прямоугольной области

f = @(x,y) sinc (x) * sinc (y));
q = integral2 (f, -1, 1, -1, 1)
  ⇒ q =  12.328  (incorrect)
q = integral2 (f, -1, 1, -1, 1, "Vectorized", false)
  ⇒ q =  1.390 (correct)
f = @(x,y) sinc (x) .* sinc (y);
q = integral2 (f, -1, 1, -1, 1)
  ⇒ q =  1.390  (correct)

Первый результат некорректен, так как неэлементный оператор между функциями sinc в f создаёт нежелательные матричные умножения между внутренними массивами интегрирования, используемыми integral2. Во втором результате установка "Vectorized" в false заставляет integral2 выполнять скалярные внутренние операции для вычисления интеграла, что приводит к правильному численному результату ценой примерно в 20-кратного увеличения времени вычислений. В третьем результате векторизация интегрального слагаемого f с использованием элементарного оператора умножения даёт правильный результат без увеличения времени вычислений.

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

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

Ссылка: L.F. Shampine, MATLAB программа для квадратуры в 2D, Applied Mathematics and Computation, стр. 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,z), и все операции должны быть векторизованными, так что x, y и z принимают массивы в качестве входных данных и возвращают массивы выходных данных того же размера. (Можно предположить, что x, y и z будут либо массивами одинакового размера, либо скалярами.) Базовые интегрирующие модули будут вводить массивы точек интегрирования в f и/или использовать внутренние векторные расширения для ускорения вычислений, что может привести к непредсказуемым результатам, если f не ограничена элементарными операциями. Для тех интеграндов, где это неизбежно, описанные ниже ("Vectorized") option могут дать более надёжные результаты.

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

AbsTol

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

RelTol

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

Method

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

Vectorized

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

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

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.

Пример 3: интегрирование не-векторизованной функции по кубическому объёму

f = @(x,y) sinc (x) * sinc (y), * sinc (z);
q = integral3 (f, -1, 1, -1, 1, -1, 1)
  ⇒ q =  14.535  (incorrect)
q = integral3 (f, -1, 1, -1, 1, -1, 1, "Vectorized", false)
  ⇒ q =  1.6388 (correct)
f = @(x,y,z) sinc (x) .* sinc (y), .* sinc (z);
q = integral3 (f, -1, 1, -1, 1, -1, 1)
  ⇒ q =  1.6388  (correct)

Первый результат неверен, так как неэлементная операция между функциями sinc в f создаёт нежелательные матричные умножения между внутренними массивами интегрирования, используемыми integral3. Во втором результате установка "Vectorized" в false заставляет integral3 выполнять скалярные внутренние операции для вычисления интеграла, что приводит к правильному численному результату ценой примерно в 30-кратного увеличения времени вычисления. В третьем результате векторизация интегранда f с использованием элементарного оператора умножения даёт правильный результат без увеличения времени вычисления.

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

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

Ссылка: L.F. Shampine, MATLAB программа для квадратур в 2D, Applied Mathematics and Computation, стр. 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, то потребуется замена переменных.