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

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

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

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

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

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

AbsTol

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

RelTol

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

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

RelTol

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

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 равен произведению площади треугольника на высоту, то есть 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 * 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, необходимо выполнить замену переменных.