Журнал вычислительной математики и математической физики, 2019, T. 59, № 7, стр. 1108-1124

1. ВВЕДЕНИЕ

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

При решении сложных задач численного интегрирования с заданной точностью может иметь значение увеличение размерности решаемой задачи даже на единицу. Это может быть реализовано в случае, когда интеграл от подынтегральной функции легко вычисляется для области интегрирования, являющейся произвольным координатным параллелепипедом (т.е. имеющим вид (1.3)). В некоторых случаях, описанных ниже, для вычисления таких интегралов потребуется не более 2n арифметических операций. Достаточно общим случаем, когда такие интегралы вычисляются сравнительно легко, является тот, когда подынтегральная функция является “псевдомногочленом”, т.е. имеет вид (1.14). В этом случае можно применить метод многоэтапного дробления исходного координатного куба (в общем случае – параллелепипеда; см. замечание 3) $\Pi $, описанного около области интегрирования X (путем последовательного деления ребер куба $\Pi $ пополам), на координатные кубы меньшего размера и использовать так называемые условия “отсечения–погружения”, идейно близкие к методам “ветвей и границ”, используемым в задачах дискретной оптимизации (см., например, [1]), и некоторых непрерывных задачах глобальной оптимизации (см., например, [2]–[5]). Такой подход позволяет сузить область “сложного” интегрирования до границы интегрируемой области (см. рис. 8 из [4], или рис. 1 из [6]), т.е. снизить размерность задачи на 1. Впервые эти условия были описаны автором в случае многомерного интегрирования в работах [4], [6], [7].

Кроме того, метод многоэтапного дробления позволяет в случае, когда область интегрирования $X$ имеет вид (1.2), где ${{g}_{i}}(x)$ – достаточно гладкие функции, производить линеаризацию ограничений ${{g}_{i}}(x) \leqslant 0$, $i = 1,2, \ldots ,r$, на кубах последнего этапа дробления, т.е. заменять их на линейные ограничения. Отметим, что с увеличением количества этапов дробления количество “активных” ограничений уменьшается и наибольший вклад дают кубы с количеством “активных” ограничений 0, затем – с количеством “активных” ограничений 1, затем с количеством “активных” ограничений 2 и т.д. Под “активным” или существенным ограничением для данного куба будем понимать ограничение такое, что не все точки этого куба удовлетворяют этому ограничению. При этом нередко удается определить с высокой точностью нижнюю и верхнюю оценки для значения интеграла на кубах последнего этапа дробления с небольшим количеством (1–2) “активных” ограничений. В этом случае используются рекурсивные алгоритмы точного вычисления объема кубов с одним или двумя линейными ограничениями (см. алгоритмы 3, 4), и указанные интегралы вычисляются с достаточно малой погрешностью. Объем объединения остальных кубов весьма мал, и сумму интегралов по этим кубам (вычисляемых уже без учета ограничений) можно отнести к погрешности метода. Все это дает возможность в некоторых случаях организовать вычисление значения интеграла 2-го порядка точности, т.е. с точностью порядка $O({{\tau }^{2}})$ при достаточно гладких функциях ${{g}_{i}}(x)$ и подынтегральной функцией, являющейся достаточно гладким “псевдомногочленом”, и 3-го порядка точности в некоторых случаях, когда подынтегральная функция является достаточно гладким “псевдомногочленом” и областью интегрирования является многогранник, т.е. с точностью порядка $O({{\tau }^{3}})$, где $\tau $ – полудлина ребер кубов последнего этапа дробления, что и показывается в настоящей работе на примере, когда подынтегральной функцией является плотность нормального распределения.

Решение многих практических задач приводит к вычислению значений вероятностных критериев, которые относятся к сложным статистическим критериям (см. [8]) качества систем управления при наличии случайных факторов. Наиболее распространенными вероятностными критериями качества технических систем при случайных воздействиях являются функционалы квантили (см. [9]) и вероятности. Например, квантильный критерий может характеризовать точность системы управления при случайных возмущениях, гарантируемую с заданной вероятностью (см. [9]), гарантируемую площадь взлетно-посадочной полосы, на которую летательный аппарат садится с заданной вероятностью в условиях ветровых помех (см. [10]) и т.д. Функционал вероятности характеризует обычно вероятность выполнения каких-либо условий рассматриваемой задачи.

Во многих случаях, как отмечено в [9], решение прямой вероятностной задачи анализа сводится к вычислению многомерного интеграла от плотности вероятности некоторого случайного вектора по заданному множеству, что представляет собой сложную вычислительную проблему. В инженерной практике в таких случаях часто используется метод статистических испытаний (Монте-Карло), обладающий простотой и универсальностью. Однако при использовании этого метода требуются большие вычислительные затраты, и при этом заданная точность вычисления интеграла не гарантируется. В связи с этим имеет смысл разработка специальных методов, применимых к некоторым достаточно узким классам подынтегральных функций, но тем не менее достаточно важным для практических приложений.

Предлагаемый метод интегрирования применялся ранее в работах [4], [6], [7]. В работе [6] он применялся для случая, когда областью интегрирования является эллипсоид, а подынтегральной функцией – плотность нормального распределения. Для приведенного в работе [6] численного метода показано, что погрешность метода является величиной порядка $O({{\tau }^{2}})$, где $\tau $ – полудлина ребер кубов последнего этапа дробления, т.е. при увеличении количества этапов дробления на 1 погрешность уменьшается примерно в 4 раза (см. табл. 1, 2 из [6]). Этот результат не зависит от размерности задачи n – числа переменных. В настоящей работе аналогичный результат описан для случая, когда областью интегрирования является многогранник (см. разд. 4). Более того, в разд. 5 для случая с многогранником приводится численный метод третьего порядка точности, т.е. с погрешностью, являющейся величиной порядка $O({{\tau }^{3}})$ независимо от числа переменных n. В разд. 2 подробно излагается алгоритм “дробления–отсечения–погружения”, его модификации, приводятся двухсторонние оценки вычисляемого многомерного интеграла. В разд. 3 рассматривается случай, когда подынтегральная функция является гауссовской плотностью вероятности. Для этого случая (или аналогичного ему, когда подынтегральная функция имеет вид $f(x) = \prod\nolimits_{i = 1}^n {{{f}_{0}}({{x}_{i}})} $, где ${{f}_{0}}(t)$ – непрерывная на соответствующем отрезке функция (см. разд. 3)) можно заранее сформировать массив значений соответствующих одномерных интегралов, и тогда вычисление значений интегралов для кубов произвольного этапа дробления будет сводиться к извлечению из массива заранее вычисленных значений и произведению над ними не более 2n арифметических операций.

Следует однако отметить, что предлагаемые методы являются эффективными лишь для задач малой размерности n ($n \leqslant 5$), поскольку при решении задач большей размерности вычислительные затраты неприемлемо велики и в этом случае следует искать более эффективные методы.

Аналогичные задачи рассматривались и в других работах, в частности, в работе [13], но при этом использовались иные алгоритмы, не позволяющие описывать их для произвольного n. Например, случай с эллипсоидом описан в [13] для n = 2, а с многогранником для n = 3. Между тем, в [6] приводятся результаты численного эксперимента для эллипсоида при n = 3 (см. табл. 1, 2 ).

2. ОБЩАЯ СХЕМА МЕТОДА

Рассматривается задача приближенного вычисления интеграла

где $f,{{g}_{1}}, \ldots ,{{g}_{r}}$ – непрерывные на $\Pi $ функции, $\Pi $ является n-мерным кубом (см. замечание 3). Пусть для простоты $f(x) \geqslant 0$  $\forall x \in \Pi $. В дальнейшем это предположение будет снято.

Часто оказывается (см. замечание 4), что для произвольного куба $Q \subseteq \Pi $ вида (1.3) сравнительно легко вычисляется значение

В этом случае для приближенного вычисления величины $I$ можно воспользоваться следующим алгоритмом.

Алгоритм 1

Шаг 1. Положим $\tilde {I} = 0$, $\delta = 0$. Выберем k – количество этапов дробления куба $\Pi $. Затем дробим $\Pi $ на ${{2}^{n}}$ кубов (делением ребер куба $\Pi $ пополам), которые назовем кубами первого этапа дробления (соответственно кубы, являющиеся элементами аналогичного дробления кубов первого этапа дробления, назовем кубами второго этапа дробления и т.д.). Положим $j = 1$ ($j$ – номер текущего этапа дробления, $1 \leqslant j \leqslant k$). Пусть ${{V}_{1}}$ – множество кубов первого этапа дробления.

Шаг 2. Если ${{V}_{j}} \ne \emptyset $ (${{V}_{j}}$ – текущее множество кубов $j$-го этапа дробления, содержащее не более ${{2}^{n}}$ кубов), то перейдем к шагу 3. В противном случае присвоим $j: = j - 1$. Если $j = 0$, то закончим работу алгоритма; в противном случае повторим выполнение шага 2.

Шаг 3. Выбрав некоторый куб $Q$ из ${{V}_{j}},$ проверим выполнение условия “погружения”

Если (1.5) выполняется, то положим исключим куб $Q$ из ${{V}_{j}}$ и перейдем к шагу 2. В противном случае проверим выполнение условия “отсечения” Если (1.7) выполняется, то исключим куб $Q$ из ${{V}_{j}}$ и перейдем к шагу 2. Если для куба $Q$ не выполняется ни условие погружения (1.5), ни условие отсечения (1.7), то перейдем к шагу 4.

Шаг 4. В случае $j = k$ положим

исключим куб $Q$ из ${{V}_{k}}$ и перейдем к шагу 2. В противном случае дробим куб $Q$ (аналогично дроблению куба $\Pi $ на ${{2}^{n}}$ кубов первого этапа дробления), в результате чего получим множество ${{V}_{{j + 1}}}$, состоящее из ${{2}^{n}}$ кубов $(j + 1)$-го этапа дробления. Присвоим $j: = j + 1$ и перейдем к шагу 3.

Замечание 1. В общем случае проверка условия (1.5) является очень сложной операцией. Однако, если область интегрирования является объединением небольшого числа выпуклых попарно не пересекающихся областей, то эта проверка осуществляется достаточно легко. Здесь важно, чтобы проверка принадлежности некоторой точки любой из этих областей была достаточно несложной.

Замечание 2. Нетрудно показать (см. утверждение 1 из [7]), что одновременное выполнение условий погружения (1.5) и отсечения (1.7) для некоторого куба $Q \subseteq \Pi $ вида (1.3) равносильно тому, что $\partial X \cap \operatorname{int} Q \ne \emptyset $, где $\partial X$ – граница множества $X.$

Оценим погрешность метода. После выполнения алгоритма 1 справедливо

Здесь и далее суммирование осуществляется по всем кубам $Q$   $k$-го (т.е. последнего) этапа дробления, для которых $\partial X \cap \operatorname{int} Q \ne \emptyset $ (т.е. для которых не выполняется ни условие отсечения, ни условие погружения; см. замечание 2).

Нетрудно показать (см. утверждение 7 из [4]), что в случае $\mu (\{ x \in \Pi \,|\,g(x) = 0\} ) = 0$, где $g(x) = \max \{ {{g}_{1}}(x),\; \ldots ,\;{{g}_{r}}(x)\} $, $\mu (Y)$ – мера Лебега множества $Y \subset {{{\text{R}}}^{n}}$, выполняется

Замечание 3. Алгоритм 1 остается в силе, если $\Pi $ является n-мерным параллелепипедом (т.е. в случае, когда для некоторых  $i \in \{ 2, \ldots ,n\} $   ${{b}_{i}} - {{a}_{i}} \ne {{b}_{1}} - {{a}_{1}}$). В этом случае элементами многоэтапного дробления $\Pi $ будут не кубы, а параллелепипеды. Однако наиболее удобные программные реализации алгоритма 1 описываются для случая, когда $\Pi $ является n-мерным кубом и соответственно элементами многоэтапного дробления $\Pi $ будут кубы. Приведем некоторые рассуждения, показывающие, что условие, когда $\Pi $ – параллелепипед является мало ограничительным, и в этом случае небольшая модификация алгоритма 1 приводит к тому, что элементами многоэтапного дробления некоторого куба $\Pi {\text{'}}$, включающего в себя $\Pi $ будут являться кубы. Действительно, в случае, когда $\Pi $ – параллелепипед, рассмотрим куб $\Pi {\text{'}} = \{ x \in {{{\text{R}}}^{n}}\,|\,{{a}_{i}} \leqslant {{x}_{i}} \leqslant b_{i}^{'},$ $i = 1,2, \ldots ,n\} $, где $b_{i}^{'} = {{a}_{i}} + \mathop {\max }\limits_{1 \leqslant i \leqslant n} \{ {{b}_{i}} - {{a}_{i}}\} $. Тогда $b_{i}^{'} - {{a}_{i}} = \mathop {\max }\limits_{1 \leqslant i \leqslant n} \{ {{b}_{i}} - {{a}_{i}}\} $, $i = 1,2, \ldots ,n$, $\Pi \subseteq \Pi {\text{'}}$, а следовательно, $X = \{ x \in \Pi {\text{'}}\,|\,x \in \Pi ,{{g}_{1}}(x) \leqslant 0, \ldots ,{{g}_{r}}(x) \leqslant 0\} $, и для нахождения приближенного значения интеграла $I = \int_X {f(x)dx} $ можно использовать многоэтапное дробление куба $\Pi {\text{'}}$. Заметим, однако, что мера куба $\Pi {\text{'}}$ может многократно превосходить меру параллелепипеда $\Pi $. В связи с этим модифицируем шаг 3 алгоритма 1. Будем теперь перед проверкой выполнения условий (1.5), (1.7) дополнительно проверять выполнение условия

Если условие (1.12) выполняется, то справедливо условие (1.7), куб $Q$ исключается на шаге 3 из ${{V}_{j}}$ и тогда перейдем к шагу 2. Такая модификация шага 3 позволит за незначительное время отсечь “лишние” области в кубе $\Pi {\text{'}}$ (по сравнению с $\Pi $). Между тем, условие (1.12) проверяется тривиальным образом за $O(n)$ простых операций. Действительно, пусть $Q = \{ x \in {{{\text{R}}}^{n}}\,|\,{{\alpha }_{i}} \leqslant {{x}_{i}} \leqslant {{\beta }_{i}},$ $i = 1,2, \ldots ,n\} $. Тогда условие (1.12) выполняется тогда и только тогда, когда

Замечание 4. Сделанное выше предположение о наличии достаточно простой процедуры вычисления интеграла (1.4), где $Q$ – произвольный n-мерный куб (или параллелепипед) выполняется в случае, когда справедливо разложение

где ${{f}_{{ij}}}(t)$ – непрерывные на $[{{a}_{j}},{{b}_{j}}]$ функции. В частности, разложению (1.14) удовлетворяет произвольный многочлен, с переменными ${{x}_{1}}, \ldots ,{{x}_{n}}$. Заметим в связи с этим, что в некоторых случаях непрерывная функция может быть приближена многочленом с небольшим числом слагаемых. Другим важным для приложений является случай, когда подынтегральная функция является плотностью нормального распределения, т.е. $f(x) = {{(2\pi )}^{{ - n/2}}}{{e}^{{ - |x{{|}^{2}}/2}}}$, где ${{\left| x \right|}^{2}} = x_{1}^{2} + \ldots + x_{n}^{2}$.

Алгоритм 1 (а также описанную в замечании 3 модификацию этого алгоритма) можно дополнительно модифицировать, с тем чтобы уменьшить погрешность $\delta $ при том же количестве этапов дробления $k$.

Модифицируем шаг 4 алгоритма 1. Теперь в случае $j = k$ определим величины ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{s} }_{Q}}$, ${{\bar {s}}_{Q}}$, ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{b} }_{Q}}$, ${{\bar {b}}_{Q}}$, удовлетворяющие условиям

положим исключим куб $Q$ из ${{V}_{k}}$ и перейдем к шагу 2.

В результате такой модификации получим из алгоритма 1 новый алгоритм 2, который будем считать основным алгоритмом. Очевидно, что после выполнения алгоритма 2 снова справедливы неравенства (1.9) и выполняется

Замечание 5. Ранее предполагалось, что $f(x) \geqslant 0$, $\forall x \in \Pi $. Пусть теперь это условие не выполняется. В этом случае модифицируем алгоритм 2 и на шаге 4 этого алгоритма вместо (1.17) положим

Заметим, что выполняются неравенства

из которых следует, что после выполнения модифицированного алгоритма 2, в котором при выполнении шага 4 вместо присвоения по формулам (1.17) осуществляется присвоение по формулам (1.19), (1.20) справедливы неравенства (1.9) и при этом

Замечание 6. В ряде случаев проверка выполнения условия погружения (1.5) или условия отсечения (1.7) оказывается очень трудоемкой (или даже практически невыполнимой). В этих случаях часто оказывается более выгодным заменить условия (1.5), (1.7) на легко проверяемые условия, достаточные для (1.5), (1.7) соответственно. Широкий набор таких условий приведен в работе [4]. Напротив, условия (1.5), (1.7) можно заменить на более слабые: $\mu (Q \cap X) = \mu (Q)$, $\mu (Q \cap X) = 0$ соответственно. Однако указанные условия являются практически проверяемыми лишь в очень редких случаях.

3. СЛУЧАЙ, КОГДА ПОДЫНТЕГРАЛЬНАЯ ФУНКЦИЯ ЯВЛЯЕТСЯ ПЛОТНОСТЬЮ НОРМАЛЬНОГО РАСПРЕДЕЛЕНИЯ

Рассмотрим теперь задачу вычисления величины $I,$ удовлетворяющей (1.1), в случае, когда

Тогда, очевидно, справедливо разложение (1.14). Приведем некоторые уточнения алгоритма 2 для случая (2.1).

Заметим, что в случае (2.1) для практической реализации алгоритма 2 необходимо:

а) достаточно быстро находить величины ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{b} }_{Q}}$, ${{\bar {b}}_{Q}}$, удовлетворяющие (1.16), причем разность ${{\bar {b}}_{Q}} - {{\bar {b}}_{Q}}$ должна быть как можно меньшей;

б) достаточно быстро вычислять $\int_Q {{{\varphi }_{n}}(x)dx} $, где $Q$ – произвольный куб, являющийся элементом многоэтапного дробления куба $\Pi $;

в) достаточно быстро находить величины ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{s} }_{Q}}$, ${{\bar {s}}_{Q}}$, удовлетворяющие (1.15), причем разность ${{\bar {s}}_{Q}} - {{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{s} }_{Q}}$ должна быть как можно меньшей;

г) достаточно быстро проверять выполнение условий (1.5), (1.7) (или условий, достаточных для (1.5), (1.7) соответственно; см. замечание 6).

Опишем методы решения перечисленных задач а)–г).

Рассмотрим сначала задачу нахождения величин ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{b} }_{Q}}$, ${{\bar {b}}_{Q}}$, удовлетворяющих (1.16), где $Q$ – куб k-го (т.е. последнего) этапа дробления. Будем искать величины ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{b} }_{Q}}$, ${{\bar {b}}_{Q}}$, удовлетворяющие равенствам:

Эти величины находятся за $O(n)$ арифметических операций. Пусть $Q = \{ x \in {{{\text{R}}}^{n}}\,|\,{{\alpha }_{i}} \leqslant {{x}_{i}} \leqslant {{\beta }_{i}},$ $i = 1,2, \ldots ,n\} ,$ $c(Q) = (({{\alpha }_{1}} + {{\beta }_{1}}){\text{/}}2,\; \ldots ,\;({{\alpha }_{n}} + {{\beta }_{n}}){\text{/}}2),$$i = 1,2, \ldots ,n$. Тогда ${{\varphi }_{n}}(u) = \min {{\varphi }_{n}}(Q)$, ${{\varphi }_{n}}(\text{v}) = \max {{\varphi }_{n}}(Q)$, т.е. можно положить ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{b} }_{Q}} = {{\varphi }_{n}}(u)$, ${{\bar {b}}_{Q}} = {{\varphi }_{n}}(\text{v})$.

Используя (2.2), уточним для рассматриваемого случая (2.1) оценку погрешности $\delta .$ Предварительно обозначим через $\tau $ половину длины ребер кубов $k$-го (т.е. последнего) этапа дробления куба $\Pi $ (указанное обозначение будем использовать и всюду далее), т.е. $\tau = T{\text{/}}{{2}^{{k + 1}}}$, где $T = {{b}_{1}} - {{a}_{1}}$ – длина ребер куба $\Pi $. Заметим, что (см. утверждение 2 из [7]) для величины $L = {{\sqrt n } \mathord{\left/ {\vphantom {{\sqrt n } {[{{{(2\pi )}}^{{n/2}}}\sqrt e ]}}} \right. \kern-0em} {[{{{(2\pi )}}^{{n/2}}}\sqrt e ]}}$ выполняется $\forall x,x{\text{'}} \in {{{\text{R}}}^{n}}$   $\left| {{{\varphi }_{n}}(x) - {{\varphi }_{n}}(x{\text{'}})} \right| \leqslant L\left\| {x - x{\text{'}}} \right\|$, где $\left\| x \right\| = \max \{ \left| {{{x}_{1}}} \right|,\; \ldots ,\;\left| {{{x}_{n}}} \right|\} $ – норма вектора $x \in {{{\text{R}}}^{n}}$, а следовательно, учитывая (2.2), получаем ${{\bar {b}}_{Q}} \leqslant \max {{\varphi }_{n}}({{{\text{R}}}^{n}}) = 1{\text{/}}{{(2\pi )}^{{n/2}}}$, ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{b} }_{Q}} - {{\bar {b}}_{Q}} \leqslant 2L\tau $, откуда, используя (1.18), (1.15), имеем

Опишем теперь следующий метод, обеспечивающий очень быстрое вычисление (в процессе работы алгоритма 2) интегралов $\int_Q {{{\varphi }_{n}}(x)dx} $, для кубов $Q$, являющихся элементами многоэтапного дробления куба $\Pi $. При этом будем считать, что $\Pi $ является кубом следующего вида:

Рассуждения из замечания 3 показывают, что такое предположение не является ограничительным. Действительно, если первоначально $\Pi $ в (1.2) является не кубом, а n-мерным параллелепипедом, то рассмотрим куб $\Pi {\text{'}}$ вида (2.4), где $a = \min \{ {{a}_{1}}, \ldots ,{{a}_{n}}\} $, $b = \max \{ {{b}_{1}}, \ldots ,{{b}_{n}}\} $. Тогда $\Pi \subseteq \Pi {\text{'}}$, а следовательно, $X = \{ x \in \Pi {\text{'}}\,|\,x \in \Pi ,\;{{g}_{1}}(x) \leqslant 0, \ldots ,{{g}_{r}}(x) \leqslant 0\} $, и для нахождения приближенного значения интеграла $I = \int_X {{{\varphi }_{n}}(x)dx} $ можно использовать многоэтапное дробление куба $\Pi {\text{'}}$. Заметим, однако, что мера куба $\Pi {\text{'}}$ может многократно превосходить меру параллелепипеда $\Pi $. В связи с этим модифицируем шаг 3 алгоритма 1 (а следовательно, и алгоритма 2 в соответствии с замечанием 3).

Итак, пусть теперь куб $\Pi $ имеет вид (2.4). Тогда перед началом работы алгоритма 2 сформируем множество (одномерный массив) $\Omega = \{ \Phi (a + (b - a)i{\text{/}}{{2}^{k}})\,|\,i = 0,1, \ldots ,{{2}^{k}}\} $, где $\Phi (y) = \frac{1}{{\sqrt {2\pi } }}\int_{ - \infty }^y {{{e}^{{ - {{t}^{2}}/2}}}} dt$. Заметим теперь, что если $Q = \{ x \in {{{\text{R}}}^{n}}\,|\,{{\alpha }_{i}} \leqslant {{x}_{i}} \leqslant {{\beta }_{i}},\;i = 1,2, \ldots ,n\} $ – куб, являющийся элементом многоэтапного дробления куба $\Pi $ вида (2.4), то

и при этом $\Phi ({{\beta }_{i}}),$ $\Phi ({{\alpha }_{i}}) \in \Omega $, $i = 1,2,\; \ldots ,\;n$. Таким образом, для вычисления интеграла (2.5) достаточно извлечь $2n$ элементов из заранее заготовленного множества $\Omega $ и произвести над ними n операций вычитания и n – 1 операций умножения.

Аналогичным приемом можно воспользоваться и в случае, когда имеет место более общее разложение подынтегральной функции $f(x) = \prod\nolimits_{i = 1}^n {{{f}_{0}}({{x}_{i}})} $, где ${{f}_{0}}(t)$ – непрерывная на $[a,{\text{ }}b]$ функция. Соответственно, в случае, когда справедливо еще более общее разложение $f(x) = \prod\nolimits_{i = 1}^n {{{f}_{i}}({{x}_{i}})} $, где ${{f}_{i}}(t)$ – непрерывные на $[a,{\text{ }}b]$ функции, нам понадобятся n множеств, аналогичных $\Omega $, а в случае (1.17) – sn множеств. Их можно вычислить заранее для достаточно большого k и хранить на цифровом накопителе для решения быть может не одной, а группы задач.

4. СЛУЧАЙ, КОГДА ОБЛАСТЬ ИНТЕГРИРОВАНИЯ ЯВЛЯЕТСЯ МНОГОГРАННИКОМ

Пусть $X = \{ x \in \Pi \,|\,{{l}_{i}}(x) \leqslant 0,\;i = 1,2, \ldots ,r\} $, ${{l}_{i}}(x) = \left\langle {{{e}^{{(i)}}},x} \right\rangle + {{d}_{i}}$, ${{e}^{{(i)}}} \in {{{\text{R}}}^{n}}$, ${{d}_{i}} \in {\text{R}}$, ${{e}^{{(i)}}} \ne {\mathbf{0}} \in {{{\text{R}}}^{n}}$, $i = 1,2, \ldots ,r$. Будем предполагать, что выполняются условия:

Выполнение условия (3.2) означает, что каждое из ограничений ${{l}_{i}}(x) \leqslant 0$, где $i \in \{ 1,2, \ldots ,r\} $, является существенным, так как в результате отбрасывания этого ограничения получаем множество, отличное от $X$. Проверка выполнения условий (3.1), (3.2) осуществляется в результате решения нескольких задач линейного программирования (ЗЛП).

Для выполнения (3.1) необходимо и достаточно, чтобы значение минимума в ЗЛП

было отрицательным.

Соответственно, для выполнения (3.2) необходима и достаточна строгая положительность значений максимума в ЗЛП

где $j = 1,2, \ldots ,r$.

Рассмотрим вопрос о проверке условия погружения (1.5). Пусть

Тогда для произвольного куба (или параллелепипеда) очевидно, выполняется $Q \subseteq X \Leftrightarrow Q \subseteq {{X}_{j}}$, $j = 1,2, \ldots ,r$, а следовательно, для проверки условия погружения (1.5) достаточно выполнить следующие действия. Для каждого $j \in \{ 1,2, \ldots ,r\} $ находим угловую точку $\text{v}(j)$ куба $Q$, на которой достигается максимум в задаче ${{l}_{j}}(x) \to \max $; $x \in Q$. Очевидно, что $i = 1,2, \ldots ,n$, т.е. точка $\text{v}(j)$ находится за $O(n)$ элементарных операций. Тогда, если ${{l}_{j}}(\text{v}(j)) \leqslant 0$, то $Q \subseteq {{X}_{j}}$, в противном случае В случае выполнения (3.3) условие (1.5) не является выполненным. Если же $Q \subseteq {{X}_{j}}$ при всех $j \in \{ 1,2, \ldots ,r\} $, то условие (1.5) выполняется.

Будем  называть ограничение ${{l}_{j}}(x) \leqslant 0$ такое, что выполняется (3.3), существенным для $Q$, а множество номеров существенных для $Q$ ограничений обозначать через ${{J}_{Q}},$ т.е. ${{J}_{Q}} = \{ j \in \{ 1,2, \ldots ,r\} \,|\,Q{\backslash }{{X}_{j}} \ne \emptyset \} $. Заметим, что если для некоторого куба $Q$, являющегося элементом многоэтапного дробления куба $\Pi $, выполняется $Q \subseteq {{X}_{j}}$ (где $j \in \{ 1,2, \ldots ,r\} $), то для любого куба $Q{\text{'}}$, являющегося элементом дальнейшего дробления куба $Q$, тем более будет выполняться включение $Q{\text{'}} \subseteq {{X}_{j}}$ и поэтому при проверке условия $Q{\text{'}} \subseteq X$ достаточно ограничиться проверкой включений $Q{\text{'}} \subseteq {{X}_{j}}$, $j \in {{J}_{Q}}.$ При таком подходе получаем существенную экономию в вычислениях при проверке условия погружения (1.5). Таким образом, проверка условия погружения (1.5) осуществляется точно и достаточно экономично.

Для проверки выполнения условия отсечения (1.7) достаточно решить задачу

Тогда, если $\gamma \geqslant 0,$ то (1.7) выполняется, а если $\gamma < 0,$ то нет. Заметим, что эта задача сводится к ЗЛП, которая может быть точно решена с помощью симплекс-метода. Однако при этом проверка условия (1.7) оказывается весьма трудоемкой операцией. В связи с этим имеет смысл заменить проверку условия (1.7) на проверку легко проверяемого условия, достаточного для (1.7) (см. замечание 6). Действительно, рассмотрим следующее легко проверяемое условие: (в этом условии участвуют только существенные на $Q$ ограничения). Очевидно, что из (3.4) следует (1.7), т.е. условие (3.4) является достаточным для (1.7). Проверка выполнения условия (3.4) осуществляется следующим образом. Для каждого $j \in {{J}_{Q}}$ находим угловую точку $u(j)$ куба $Q,$ на которой достигается минимум в задаче ${{l}_{j}}(x) \to \min $; $x \in Q$. Очевидно, что $i = 1,2, \ldots ,n$, т.е. эта точка находится за $O(n)$ элементарных операций. Тогда, если ${{l}_{j}}(u(j)) \geqslant 0$, то условие (3.4) (а следовательно, и условие (1.7) ) выполняется. Если же ${{l}_{j}}(u(j)) < 0$ при всех $j \in {{J}_{Q}}$, то условие (3.4) не выполняется.

Обозначим через ${{G}_{\tau }}$ множество кубов $Q$ последнего этапа дробления куба $\Pi $, для которых не выполняется ни условие погружения (1.5), ни условие (3.4) (являющееся достаточным для условия отсечения (1.7)). Тогда, если заменить на шаге 3 алгоритма 2 проверку выполнения условия (1.7) на проверку выполнения условия (3.4) (в соответствии с замечанием 6), то после выполнения алгоритма 2, очевидно, будут справедливы неравенства (1.9), а погрешность $\delta $ будет удовлетворять равенству, аналогичному (1.18), но в котором суммирование производится по кубам $Q \in {{G}_{\tau }}$. Соответственно, в рассматриваемом случае (2.1) будет выполняться неравенство (аналогичное неравенству (2.3))

Замечание 7. В случае $\left| {{{J}_{Q}}} \right| = 1$ справедливо следующее: (3.4) $ \Leftrightarrow $ (1.7), т.е. в этом случае (наиболее частом; см. утверждения 16, 17 из [7]) замена (1.7) на (3.4) тем не менее обеспечивает точную проверку выполнения условия (1.7). В случае же нескольких существенных ограничений может оказаться, что условие (3.4) не справедливо, но условие (1.7) выполняется, и тогда нам придется производить дополнительные дробления куба $Q$ (в случае, если $Q$ не является кубом последнего этапа дробления), или дополнительно вычислять для $Q$ величины ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{s} }_{Q}}$, ${{\bar {s}}_{Q}}$, ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{b} }_{Q}}$, ${{\bar {b}}_{Q}}$ (в случае, если $Q$ является кубом последнего этапа дробления). Практика показывает, что таких дополнительных дроблений обычно оказывается не более двух. На рис. 3 из [7] приведен пример, когда выполняется условие (1.7), но не выполняется условие (3.4). Однако для кубов, являющихся элементами дробления куба, изображенного на этом рисунке, условие (3.4) уже выполняется.

Рассмотрим теперь вопрос о вычислении для кубов $Q \in {{G}_{\tau }}$ величин ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{s} }_{Q}}$, ${{\bar {s}}_{Q}}$, удовлетворяющих (1.15). Предварительно опишем алгоритм вычисления меры множества вида ${{V}^{{(n)}}} = \{ x \in Q\,|\,\left\langle {e,x} \right\rangle + d \leqslant 0\} $, где $Q = \{ x \in {{{\text{R}}}^{n}}\,|\,{{\alpha }_{i}} \leqslant {{x}_{i}} \leqslant {{\beta }_{i}},\;i = 1,2, \ldots ,n\} $ – произвольный n-мерный куб (или параллелепипед), ${{\alpha }_{i}} < {{\beta }_{i}}$, $i = 1,2, \ldots ,n$, $e \in {{{\text{R}}}^{n}},$ $\left| e \right| > 0$ (случай, когда $\left| e \right| = 0$ тривиален).

Алгоритм 3

Шаг 1. Определим точки $u$, $\text{v}$, на которых соответственно достигаются минимум и максимум функции $\left\langle {e,x} \right\rangle + d$ на $Q$:

$i = 1,2, \ldots ,n$. Если $\left\langle {e,\text{v}} \right\rangle + d \leqslant 0$, то ${{V}^{{(n)}}} = Q$ и тогда $\mu ({{V}^{{(n)}}}) = \mu (Q)$. Если $\left\langle {e,u} \right\rangle + d \geqslant 0$, то ${{V}^{{(n)}}} \cap \operatorname{int} Q = \emptyset $ и тогда $\mu ({{V}^{{(n)}}}) = 0$. Если нам еще не удалось определить $\mu ({{V}^{{(n)}}})$, то перейдем к шагу 2.

Шаг 2. Если $n = 1$, то ${{V}^{{(n)}}}$ является прямолинейным отрезком, а $\mu ({{V}^{{(n)}}})$ – разностью его верхней и нижней границ, которые находятся тривиальным образом. Если $n > 1$, то перейдем к шагу 3.

Шаг 3. Определим любую точку $w$ на гиперплоскости, заданной равенством $\left\langle {e,x} \right\rangle + d = 0$ (например, точку пересечения этой гиперплоскости с прямой, проходящей через точки $u$, $\text{v}$). Она найдется, поскольку $\left\langle {e,\text{v}} \right\rangle + d > 0$, $\left\langle {e,u} \right\rangle + d < 0$ (см. шаг 1 алгоритма 3). Пусть

(${{H}_{{ij}}}$ являются ($n - 1$)-мерными гиперплоскостями, проходящими через грани куба $Q$), где $i = 1,2, \ldots ,n$,$j = 1,2$. Определим числа ($\left| {{{h}_{{ij}}}} \right| = \rho (w,{{H}_{{ij}}})$ – расстояние от $w$ до ${{H}_{{ij}}}$), где $i = 1,2, \ldots ,n$, $j = 1,2$.

Шаг 4 (рекурсия по $n$). Осуществим сведение задачи вычисления $\mu ({{V}^{{(n)}}})$ к $2n$ аналогичным задачам размерности $n - 1$ (см. рис. 2 из [7])

где При вычислении $\mu (V_{{ij}}^{{(n - 1)}})$ снова действуем согласно шагам 1–4, присвоив $n: = n - 1$, и т.д.

Рассмотрим теперь вопрос о вычислении для кубов $Q \in {{G}_{\tau }}$ величин ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{s} }_{Q}}$, ${{\bar {s}}_{Q}}$.

В случае $\left| {{{J}_{Q}}} \right| = 1$ положим

где в этом случае $\mu (Q \cap X)$ вычисляется точно с помощью алгоритма 3. А в случае $Q \in {{G}_{\tau }}$, $\left| {{{J}_{Q}}} \right| \geqslant 2$ положим

Оценим погрешность метода. Из утверждений 16, 17 из [7] следует, что при $\tau \to 0 + $ выполняется

откуда, используя (3.5) и выполнение (3.10), (3.11) при $\left| {{{J}_{Q}}} \right| = 1$, $\left| {{{J}_{Q}}} \right| \geqslant 2$ соответственно (где $Q \in {{G}_{\tau }}$), получаем, что $\delta = O({{\tau }^{2}})$ при $\tau \to 0{\kern 1pt} + $.

Таким образом, описан метод второго порядка точности для случая, когда множеством интегрирования является многогранник, а подынтегральная функция является плотностью нормального распределения. Ранее в работе [6] был описан метод второго порядка точности для случая, когда область интегрирования является эллипсоидом. В [6] приведены результаты вычислительных экспериментов, из которых видно, что при увеличении количества этапов дробления на 1 (что соответствует уменьшению $\tau $ ровно в два раза) погрешность $\delta $ уменьшалась примерно в 4 раза.

Автором была составлена программа для вычисления $I = \int_X {{{\varphi }_{n}}(x)dx} $ для случая, когда $X$ – многогранник, реализующая описанный метод второго порядка точности. Рассматривался пример, когда $n = 5$, $X = \{ x \in \Pi \,|\,{{x}_{1}} + {{x}_{2}} - {{x}_{3}} - {{x}_{4}} - {{x}_{5}} - 7 \leqslant 0,$ $2{{x}_{1}} - {{x}_{2}} + 2{{x}_{3}} - {{x}_{4}} + 2{{x}_{5}} - 8 \leqslant 0,$ ${{x}_{1}} - {{x}_{2}} + $ $ + \;2{{x}_{3}} - {{x}_{4}} + 2{{x}_{5}} - 9 \leqslant 0,$ $2{{x}_{1}} + {{x}_{2}} - {{x}_{3}} + {{x}_{4}} - {{x}_{5}} - 7 \leqslant 0\} $, $\Pi = \{ x \in {{{\text{R}}}^{5}}\,|\, - {\kern 1pt} 2 \leqslant {{x}_{i}} \leqslant 2,$ $i = 1,2,3,4,5\} $, т.е. при $r = 4.$ Результаты вычислений величин $\tilde {I}$, $\tilde {I} + \delta $, $\delta $, удовлетворяющих (1.9), приведены в таблице при различных $k = 4,5,6$.

Из этой таблицы видно, что погрешность $\delta $ уменьшается примерно в 4 раза при каждом увеличении $k$ на 1.

5. ОПИСАНИЕ ЧИСЛЕННОГО МЕТОДА ТРЕТЬЕГО ПОРЯДКА ТОЧНОСТИ

В случае, когда множеством интегрирования является многогранник, а подынтегральная функция является плотностью нормального распределения (последнее совсем не обязательно) удается описать метод третьего порядка точности, т.е. при использовании этого метода при увеличении количества этапов дробления на 1 (что соответствует уменьшению $\tau $ ровно в два раза) погрешность $\delta $ должна уменьшаться примерно в 8 раз.

Опишем метод третьего порядка точности. Пусть $Q = \{ x \in {{{\text{R}}}^{n}}\,|\,{{\alpha }_{i}} \leqslant {{x}_{i}} \leqslant {{\beta }_{i}},$ $i = 1,2, \ldots ,n\} \in {{G}_{\tau }}$. Обозначим ${{\tilde {\varphi }}_{n}}(x) = {{\varphi }_{n}}(c(Q)) + \left\langle {\varphi _{n}^{'}(c(Q)),\;x - c(Q)} \right\rangle $; ${{\eta }_{i}} = \max \left\{ {\left| {{{t}^{2}} - 1} \right|:t \in \left[ {{{\alpha }_{i}},{{\beta }_{i}}} \right]} \right\}$, $i = 1,2, \ldots ,n$. Очевидно, что ${{\eta }_{i}} = \max \left\{ {\left| {\alpha _{i}^{2} - 1} \right|;\left| {\beta _{i}^{2} - 1} \right|} \right\}$, если $0 \notin \left[ {{{\alpha }_{i}},{{\beta }_{i}}} \right]$, и ${{\eta }_{i}} = \max \left\{ {\left| {\alpha _{i}^{2} - 1} \right|;\left| {\beta _{i}^{2} - 1} \right|;1} \right\}$, если $0 \in \left[ {{{\alpha }_{i}},{{\beta }_{i}}} \right]$, $i = 1,2, \ldots ,n$. Нетрудно показать (см. утверждение 25 из [7]), что $\forall x \in Q$   $\left| {{{{\tilde {\varphi }}}_{n}}(x) - {{\varphi }_{n}}(x)} \right| \leqslant {{a}_{Q}}{{\tau }^{2}}$, где

(очевидно, что для любого куба $Q \subseteq \Pi $ справедливо ${{a}_{Q}} \subseteq {{a}_{\Pi }},$ т.е. величина ${{a}_{Q}}$ ограничена сверху). Обозначим ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\varphi } }_{n}}(x) = {{\tilde {\varphi }}_{n}}(x) - {{a}_{Q}}{{\tau }^{2}}$, ${{\bar {\varphi }}_{n}}(x) = {{\tilde {\varphi }}_{n}}(x) + {{a}_{Q}}{{\tau }^{2}}$. Тогда и при этом ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\varphi } }_{n}}(x)$, ${{\bar {\varphi }}_{n}}(x)$ являются линейными функциями.

Модифицируем теперь алгоритм 2 следующим образом (в результате получим алгоритм 4). Для каждого куба $Q \in {{G}_{\tau }}$, такого, что $\left| {{{J}_{Q}}} \right| = 1$, присвоим

а с остальными кубами (являющимися элементами многоэтапного дробления куба $\Pi $) поступим согласно алгоритму 2. Очевидно, что после выполнения алгоритма 4 справедливы неравенства (1.9). Оценим погрешность метода (для простоты будем опускать запись $Q \in {{G}_{\tau }}$ под знаком сумм). Заметим, что (аналогично (2.3)) Пусть для кубов $Q \in {{G}_{\tau }}$ таких, что $\left| {{{J}_{Q}}} \right| \leqslant 2$, имеет место (3.10), а при $\left| {{{J}_{Q}}} \right| \geqslant 3$ – (3.11) (что обеспечивается использованием описанного ниже алгоритма 5). Тогда из (4.3), используя утверждения 16–18 из [7], имеем

Таким образом, для описания метода третьего порядка точности необходимо указать метод точного вычисления интеграла (см. (4.1))

где $e{\text{'}} \in {{{\text{R}}}^{n}}$, $d{\text{'}} \in {\text{R}}$, $Q = \{ x \in {{{\text{R}}}^{n}}\,|\,{{\alpha }_{i}} \leqslant {{x}_{i}} \leqslant {{\beta }_{i}},\;i = 1,2, \ldots ,n\} $ – некоторый куб (или в общем случае – параллелепипед), $Q \subseteq \Pi $, $\left| {{{J}_{Q}}} \right| = 1$, и не выполняется условие (3.4).

Выведем формулу для вычисления интеграла (4.4). Обозначим

где $t \in {\text{R}}$, $\nu \in {\text{N}}$, $1 \leqslant \nu \leqslant n$.

Приведем сначала формулы для вычисления интегралов

Очевидно, что ${{I}^{{(n)}}}(t) = 0$, если $t \leqslant 0$ и ${{I}^{{(n)}}}(t) = \frac{1}{{n!}}{{t}^{n}}$, если $t > 0$. Эта формула легко доказывается методом математической индукции. Кроме того, $\forall i \in \{ 1,2, \ldots ,n\} $ в случае $t \leqslant 0$   $I_{i}^{{(n)}}(t) = 0$, а в случае $t > 0$ имеем

Приведем теперь формулы для вычисления интегралов

$i = 1,2, \ldots ,n,$ $t \in {\text{R,}}$ $y \in {{{\text{R}}}^{n}}$. Используя подстановку вида $x = \text{v} + y$, где $\text{v} \in {{{\text{R}}}^{n}}$ – новые переменные, нетрудно показать, что

Приведем теперь формулы для вычисления интегралов

где $i = 1,2, \ldots ,n$, $t \in {\text{R}}$, $Q = \{ x \in {{{\text{R}}}^{n}}\,|\,{{\alpha }_{i}} \leqslant {{x}_{i}} \leqslant {{\beta }_{i}},\;i = 1,2, \ldots ,n\} $ – прямоугольный параллелепипед (${{\alpha }_{i}} < {{\beta }_{i}}$, $i = 1,2, \ldots ,n$). Для любых номеров ${{i}_{1}}, \ldots ,{{i}_{m}}$ таких, что $1 \leqslant {{i}_{1}} < \ldots < {{i}_{m}} \leqslant n$, обозначим через $y({{i}_{1}}, \ldots ,{{i}_{m}}{\text{)}}$ вектор из ${{{\text{R}}}^{n}}$ с координатами где $i = 1,2, \ldots ,n$. Пусть, кроме того, $\alpha = ({{\alpha }_{1}}, \ldots ,{{\alpha }_{n}}) \in {{{\text{R}}}^{n}}$. Тогда по формуле включений и исключений (см., например, [11]) и $\forall i \in \{ 1,2, \ldots ,n\} $

Замечание 8. Множество точек вида $y({{i}_{1}}, \ldots ,{{i}_{m}}{\text{)}}$ в приведенных формулах представляет собой множество всех угловых точек прямоугольного параллелепипеда $Q$ (за исключением $\alpha $). При этом в случае $\sum\nolimits_{i = 1}^n {y({{i}_{1}}, \ldots ,{{i}_{m}}{\text{)}}} \geqslant t$ справедливы равенства ${{I}^{{(n)}}}(y({{i}_{1}}, \ldots ,{{i}_{m}}{\text{)}},t) = 0$, $I_{i}^{{(n)}}(y({{i}_{1}}, \ldots ,{{i}_{m}}{\text{)}},t) = 0$, $i = 1,2, \ldots ,n$, что дает возможность для уменьшения объема вычислений.

Приведем далее формулы для вычисления интегралов

где $i = 1,2, \ldots ,n$, $t \in {\text{R}}$, $1 \leqslant \nu \leqslant n$, $Q = \{ x \in {{{\text{R}}}^{n}}\,|\,{{\alpha }_{i}} \leqslant {{x}_{i}} \leqslant {{\beta }_{i}},\;i = 1,2, \ldots ,n\} $.

Обозначим ${{Q}^{{(\nu )}}} = \{ \text{v} \in {{{\text{R}}}^{\nu }}\,|\,{{\alpha }_{i}} \leqslant {{\text{v}}_{i}} \leqslant {{\beta }_{i}},\;i = 1,2, \ldots ,\nu \} $. Тогда, как нетрудно видеть,

Пусть $\gamma \in {{{\text{R}}}^{n}}$, $t \in {\text{R}}$. Обозначим ${{S}^{{(n)}}}(\gamma ,t) = \left\{ {x \in {{{\text{R}}}^{n}}\,|\,\sum\limits_{i = 1}^n {{{\gamma }_{i}}{{x}_{i}}} \leqslant t} \right\}$. Рассмотрим задачу вычисления интегралов

где $i = 1,2, \ldots ,n$, $t \in {\text{R}}$, $\left| \gamma \right| \ne 0$, $Q = \{ x \in {{{\text{R}}}^{n}}\,|\,{{\alpha }_{i}} \leqslant {{x}_{i}} \leqslant {{\beta }_{i}},\;i = 1,2, \ldots ,n\} $. Не теряя общности, можно считать, что для некоторого $\nu \in \{ 1,2, \ldots ,n\} $ выполняется ${{\gamma }_{1}} \ne 0,\; \ldots ,\;{{\gamma }_{\nu }} \ne 0$, ${{\gamma }_{{\nu + 1}}} = \ldots = {{\gamma }_{n}} = 0$ (иначе перенумеруем переменные). Обозначим $Q{\text{'}} = \{ x \in {{{\text{R}}}^{n}}\,|\,\alpha _{i}^{'} \leqslant {{x}_{i}} \leqslant \beta _{i}^{'}$, $i = 1,2, \ldots ,\nu ;$ ${{\alpha }_{j}} \leqslant {{x}_{j}} \leqslant {{\beta }_{j}},$ $j = \nu + 1, \ldots ,n\} $, где Тогда, используя подстановки ${{x}_{i}} = \frac{1}{{{{\gamma }_{i}}}}{{u}_{i}}$, $i = 1,2, \ldots ,\nu $, нетрудно показать, что

Заметим теперь, что, если $Q \cap X = \{ x \in Q\,|\,\left\langle {e,x} \right\rangle + d \leqslant 0\} $, где $e \in {{{\text{R}}}^{n}}$, $d \in {\text{R}}$, $\left| e \right| \ne 0$, то

или

Для экономии объема вычислений формулу включений и исключений можно применять сразу ко всей сумме в правой части первой формулы или к вычитаемой сумме в правой части второй формулы, совершая последовательный обход угловых точек параллелепипеда $Q{\text{'}}$ (см. замечание 14 из [7]), которые с помощью простых линейных преобразований получаются из угловых точек параллелепипеда $Q$ (угловые точки $w{\text{'}}$ параллелепипеда $Q{\text{'}}$ определяются по формуле $w{\text{'}} = \Lambda w$, где $\Lambda $ – диагональная матрица порядка n с элементами ${{\Lambda }_{{ii}}} = {{\gamma }_{i}}$).

Если количество угловых точек параллелепипеда $Q$, содержащихся в множестве $\{ x \in Q\,|\,\left\langle {e,x} \right\rangle + d < 0\} $, не превосходит ${{2}^{{n - 1}}}$ (что заведомо выполняется в случае $\left\langle {e,c(Q)} \right\rangle + d \geqslant 0$, где $c(Q)$ – центр параллелепипеда $Q$, то для вычисления интеграла (4.4) имеет смысл использовать первую формулу. Соответственно, если количество угловых точек параллелепипеда $Q$, содержащихся в множестве $\{ x \in Q\,|\,\left\langle {e,x} \right\rangle + d > 0\} $, не превосходит ${{2}^{{n - 1}}}$ (что заведомо выполняется в случае $\left\langle {e,c(Q)} \right\rangle + d \leqslant 0$), то для вычисления интеграла (4.4) имеет смысл использовать вторую формулу (см. рис. 7 из [7]).

Итак, для описания алгоритма третьего порядка точности осталось привести алгоритм точного вычисления $\mu (Q \cap X)$ в случае $\left| {{{J}_{Q}}} \right| \leqslant 2$. Нам понадобится следующая

Лемма 1 (см. [12, с. 74). Пусть $V = \left\{ {x \in {{{\text{R}}}^{n}}\,|\,\left\langle {{{c}^{{(i)}}},x} \right\rangle } \right. = {{b}_{i}},\;i = 1,2, \ldots ,{{r}_{1}},{\text{ }}$ $\left. {\left\langle {{{c}^{{(j)}}},x} \right\rangle \leqslant {{b}_{j}},\;j = {{r}_{1}} + 1, \ldots ,r} \right\}$ – выпуклое полиэдральное множество, где ${{c}^{{(i)}}} \in {{{\text{R}}}^{n}}$, ${{b}_{i}} \in {\text{R}}$, $i = 1,2, \ldots ,r$, $0 \leqslant {{r}_{1}} \leqslant r$, $r \geqslant n$. Пусть, далее, ${{\text{v}}^{0}}$ – угловая точка множества $V.$ Тогда найдется множество номеров ${{J}_{0}} \subseteq \{ 1, \ldots ,r\} $ такое, что $\left| {{{J}_{0}}} \right| = n$, точка ${{\text{v}}^{0}}$ удовлетворяет системе равенств $\left\langle {{{c}^{{(i)}}},{{\text{v}}^{0}}} \right\rangle = {{b}_{i}},\;i \in {{J}_{0}}$, и определитель матрицы этой системы отличен от нуля (т.е. ${{\text{v}}^{0}}$ – единственное решение этой системы).

Обозначим

Опишем алгоритм точного вычисления $\mu (Q \cap X)$ для случая $\left\langle {{{J}_{Q}}} \right\rangle \leqslant 2$. Будем для простоты обозначений считать, что ${{J}_{Q}} \subseteq \{ 1,2\} $. Предварительно докажем следующие вспомогательные утверждения.

Лемма 2. Пусть ${{J}_{Q}} = \{ 1,2\} $, $\bar {M} = \max {{l}_{2}}({{P}_{1}}) < 0$. Тогда возможны случаи:

1) $\exists x{\text{'}} \in Q:{{l}_{1}}(x{\text{'}}) > 0,\;{{l}_{2}}(x{\text{'}}) \geqslant 0;$

2) $\forall x \in Q:{{l}_{2}}(x) \geqslant 0 \Rightarrow {{l}_{1}}(x) \leqslant 0,$

и при этом в случае 1) $\mu ({{S}_{1}} \cap {{S}_{2}}) = \mu ({{S}_{1}})$ и $\forall \text{v} \in {\text{Arg}}\max {{l}_{2}}(Q)$ выполнено ${{l}_{1}}(\text{v}) > 0$; а в случае 2) $\mu ({{S}_{1}} \cap {{S}_{2}}) = \mu ({{S}_{1}}) + \mu ({{S}_{2}}) - \mu (Q)$, $\forall \text{v} \in {\text{Arg}}\max {{l}_{2}}(Q)$ выполнено ${{l}_{1}}(\text{v}) < 0$.

Доказательство. Случаи 1), 2), очевидно, взаимно исключают, а также дополняют друг друга. Рассмотрим сначала случай 1). Возможны подслучаи:

Заметим, что из условий 1), (4.6) следует, что $\exists x* \in Q:{{l}_{1}}(x*) = 0,\;{{l}_{2}}(x*) \geqslant 0$, что противоречит условию $\bar {M} < 0.$ Рассмотрим оставшийся подслучай (4.7). Из (4.7) следует, что $\{ x \in Q\,|\,{{l}_{2}}(x) \geqslant 0\} \subseteq \{ x \in Q\,|\,{{l}_{1}}(x) \geqslant 0\} $, а в силу $\max {{l}_{2}}({{P}_{1}}) < 0$ не существует $x \in Q:{{l}_{1}}(x) = 0,\;{{l}_{2}}(x) \geqslant 0$, а следовательно,

Но тогда, используя (4.8), получаем

В случае 2) имеем $\{ x \in Q\,|\,{{l}_{2}}(x) \geqslant 0\} \subseteq \{ x \in Q\,|\,{{l}_{1}}(x) \leqslant 0\} ,$ а в силу $\bar {M} < 0$, не существует $x \in Q:{{l}_{1}}(x) = 0,\;{{l}_{2}}(x) \geqslant 0,$ а следовательно,

Но тогда, используя (4.9 ), получаем при этом пересечение множеств в последнем объединении, очевидно, является пустым, а следовательно, $\mu ({{S}_{1}}) = \mu ({{S}_{1}} \cap {{S}_{2}}) + (\mu (Q) - \mu ({{S}_{2}}))$, откуда $\mu ({{S}_{1}} \cap {{S}_{2}}) = \mu ({{S}_{1}}) + \mu ({{S}_{2}}) - \mu (Q)$.

Лемма 3. Пусть ${{J}_{Q}} = \{ 1,2\} $, $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{M} = \min {{l}_{2}}({{P}_{1}}) > 0$. Тогда возможны следующие случаи:

1) $\exists x{\text{'}} \in Q:{{l}_{1}}(x{\text{'}}) > 0,\;{{l}_{2}}(x{\text{'}}) \leqslant 0;$

2) $\forall x \in Q:{{l}_{2}}(x) \leqslant 0 \Rightarrow {{l}_{1}}(x) \leqslant 0,$

и при этом в случае 1) верно $\mu ({{S}_{1}} \cap {{S}_{2}}) = 0$, и $\forall \text{v} \in {\text{Arg}}\min {{l}_{2}}(Q)$ выполнено ${{l}_{1}}(\text{v}) > 0$; а в случае 2) верно $\mu ({{S}_{1}} \cap {{S}_{2}}) = \mu ({{S}_{2}})$, $\forall \text{v} \in {\text{Arg}}\min {{l}_{2}}(Q)$ выполнено ${{l}_{1}}(\text{v}) < 0$.

Доказательство. Случаи 1), 2), очевидно, взаимно исключают, а также дополняют друг друга. Рассмотрим сначала случай 1). Возможны подслучаи:

Заметим, что из условий 1), (4.10) следует, что $\exists x* \in Q$: ${{l}_{1}}(x*) = 0$, ${{l}_{2}}(x*) \leqslant 0$, что противоречит условию $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{M} > 0.$ Рассмотрим оставшийся подслучай (4.11). Из (4.11) следует, что $\{ x \in Q\,|\,{{l}_{2}}(x) \leqslant 0\} \subseteq \{ x \in Q\,|\,{{l}_{1}}(x) \geqslant 0\} $, а в силу $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{M} > 0$, не существует $x \in Q:{{l}_{1}}(x) = 0,\;{{l}_{2}}(x) \leqslant 0$, а следовательно,

Но тогда, используя (4.12), получаем

В случае 2) имеем $\{ x \in Q\,|\,{{l}_{2}}(x) \leqslant 0\} \subseteq \{ x \in Q\,|\,{{l}_{1}}(x) \leqslant 0\} $, а в силу $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{M} > 0$, не существует $x \in Q:{{l}_{1}}(x) = 0,\;{{l}_{2}}(x) \leqslant 0$, а следовательно,

Но тогда, используя (4.13), получаем

Опишем теперь алгоритм точного вычисления $\mu (Q \cap X)$ в случае ${{J}_{Q}} \subseteq \{ 1,2\} .$ Обозначим ${{V}^{{(n)}}} = Q \cap X$.

Алгоритм 4

Шаг 1. Определим точки $u(1)$, $u(2)$, на которых соответственно достигается минимум функций ${{l}_{1}}(x)$, ${{l}_{2}}(x)$ на $Q$ (см. шаг 1 алгоритма 3), а также точки $\text{v}(1)$, $\text{v}(2)$, на которых соответственно достигается максимум функций ${{l}_{1}}(x)$, ${{l}_{2}}(x)$ на $Q$ (см. шаг 1 алгоритма 3). Если ${{l}_{1}}(\text{v}(1)) \leqslant 0$ или ${{l}_{2}}(\text{v}(2)) \leqslant 0$ (т.е., если $\left| {{{J}_{Q}}} \right| \leqslant 1$), то применим для определения $\mu ({{V}^{{(n)}}})$ алгоритм 3. Если ${{l}_{1}}(u(1)) \geqslant 0$ или ${{l}_{2}}(u(2)) \geqslant 0$, то ${{V}^{{(n)}}} \cap \operatorname{int} Q = \emptyset $, и тогда положим $\mu ({{V}^{{(n)}}}) = 0$. Если нам еще не удалось определить $\mu ({{V}^{{(n)}}})$, то перейдем к шагу 2.

Шаг 2. Если $n = 1$, то ${{V}^{{(n)}}}$ является прямолинейным отрезком, а $\mu ({{V}^{{(n)}}})$ – разностью его верхней и нижней границ, которые находятся тривиальным образом. Если $n > 1$, то перейдем к шагу 3.

Шаг 3. Ищем множество угловых точек ${{\Omega }_{1}}$ многогранника ${{P}_{1}}$. Из леммы 1 следует, что любая угловая точка $w \in {{\Omega }_{1}}$ является либо угловой точкой куба $Q$, либо принадлежит одному из $n{{2}^{{n - 1}}}$ ребер куба $Q$, т.е. удовлетворяет равенству ${{l}_{1}}(w) = 0$ и $n - 1$ равенствам вида ${{w}_{i}} = {{\gamma }_{i}}$, где ${{\gamma }_{i}} \in \{ {{\alpha }_{i}},{{\beta }_{i}}\} $, $i = 1,2, \ldots ,n$. Но тогда для нахождения искомых угловых точек достаточно осуществить перебор (он может оказаться неполным) этих ребер. Для каждого ребра $\Gamma $ куба Q вычисляем значения функции ${{l}_{1}}(x)$ на концах $\Gamma {\text{'}}$, $\Gamma {\text{''}}$ этого ребра, т.е. значения ${{r}_{1}} = {{l}_{1}}(\Gamma {\text{'}})$, ${{r}_{2}} = {{l}_{1}}(\Gamma {\text{''}})$. Если ${{r}_{1}} = 0$, то $\Gamma {\text{'}}$ (угловая точка $Q$) включается в ${{\Omega }_{1}}$. Если ${{r}_{2}} = 0$, то $\Gamma {\text{''}}$ включается в ${{\Omega }_{1}}$. Если одновременно выполняется ${{r}_{1}} > 0$, ${{r}_{2}} > 0$ (или ${{r}_{1}} < 0$, ${{r}_{2}} < 0$), то на ребре $\Gamma $ нет угловых точек множества ${{P}_{1}}$ (т.к. $\Gamma \cap {{P}_{1}} = \emptyset $) и переходим к следующему ребру. Если ${{r}_{1}} < 0$, ${{r}_{2}} > 0$ (или ${{r}_{1}} > 0$, ${{r}_{2}} < 0$), то находим очередную угловую точку $\hat {\Gamma } = \Gamma {\text{''}} + \xi (\Gamma {\text{'}} - \Gamma {\text{''}}) \in {{\Omega }_{1}}$, где $\xi = {{r}_{2}}{\text{/}}({{r}_{2}} - {{r}_{1}}) > 0$ (очевидно, что ${{l}_{1}}(\hat {\Gamma }) = 0$).

Шаг 4. В процессе нахождения множества ${{\Omega }_{1}}$ проверяем выполнение для каждой очередной точки $w \in {{\Omega }_{1}}$ значение ${{l}_{2}}(w)$. Если ${{l}_{2}}(w) = 0$, то переходим к шагу 6. Пусть ${{\hat {l}}_{2}} = {{l}_{2}}(\hat {w})$, где$\hat {w} \in {{\Omega }_{1}}$, – максимальное, а ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{l} }_{2}} = {{l}_{2}}(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{w} )$, где $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{w} \in {{\Omega }_{1}}$, – минимальное из уже вычисленных значений ${{l}_{2}}(w)$ в точках $w \in {{\Omega }_{1}}$. Тогда в первом же случае выполнения ${{\hat {l}}_{2}} > 0,$ ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{l} }_{2}} < 0$ полагаем $w = \hat {w} + \eta (\hat {w} - \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{w} )$, где $\eta = {{\hat {l}}_{2}}{\text{/}}({{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{l} }_{2}} - {{\hat {l}}_{2}})$. Очевидно, что ${{l}_{2}}(w) = 0$, и тогда переходим к шагу 6.

В противном случае переходим к следующей точке из ${{\Omega }_{1}}$. Если все точки из ${{\Omega }_{1}}$ перебраны, то переходим к шагу 5.

Шаг 5. Переход к этому шагу означает, что либо ${{\hat {l}}_{2}} = \max {{l}_{2}}({{\Omega }_{1}}) = \max {{l}_{2}}({{P}_{1}}) < 0$, либо ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{l} }_{2}} = \min {{l}_{2}}({{\Omega }_{1}}) = \min {{l}_{2}}({{P}_{1}}) > 0$. В случае ${{\hat {l}}_{2}} < 0$ находим точку $\text{v} \in Q$ такую , что ${{l}_{2}}(\text{v}) = \max {{l}_{2}}(Q)$ (находится за $O(n)$ арифметических операций). Тогда в силу леммы 2, если ${{l}_{1}}(\text{v}) > 0$, то $\mu ({{S}_{1}} \cap {{S}_{2}}) = \mu ({{S}_{1}})$, а в случае ${{l}_{1}}(\text{v}) < 0$ выполняется $\mu ({{S}_{1}} \cap {{S}_{2}}) = \mu ({{S}_{1}}) + \mu ({{S}_{2}}) - \mu (Q)$. В случае ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{l} }_{2}} > 0$ находим точку $u \in Q$ такую, что ${{l}_{2}}(u) = \min {{l}_{2}}(Q)$. Тогда в силу леммы 3, если ${{l}_{1}}(u) > 0$, то $\mu ({{S}_{1}} \cap {{S}_{2}}) = 0$, а в случае ${{l}_{1}}(u) < 0$ выполняется $\mu ({{S}_{1}} \cap {{S}_{2}}) = \mu ({{S}_{2}})$.

Шаг 6 (рекурсия по n). Переход к этому шагу означает, что найдена точка $w \in Q:$ ${{l}_{1}}(w) = 0$, ${{l}_{2}}(w) = 0$. Пусть ${{H}_{{ij}}}$ удовлетворяют (3.6), где $i = 1,2, \ldots ,n$, $j = 1,2$. Определяем числа ${{h}_{{ij}}}$ по формулам (3.7). Осуществляем сведение задачи вычисления $\mu ({{V}^{{(n)}}})$ к 2n аналогичным задачам меньшей размерности согласно формулам (3.8), (3.9). При вычислении $\mu (V_{{ij}}^{{(n - 1)}})$ снова действуем согласно шагам 1–4, присвоив $n: = n - 1$, и т.д.

Замечание 9. Алгоритм 4, а также леммы 2, 3 (как и алгоритм 3) справедливы и в случае, когда $Q$ – произвольный параллелепипед вида (1.3).

Замечание 10. Отметим, что предложенный в [4], [7] подход позволяет, например, описать простой конечный алгоритм нахождения точного объема любого многогранника (см. алгоритм 6 из [7]).

Замечание 11. Отметим, что при программной реализации метода имеется возможность для распараллеливания вычислительного процесса на этапе дробления кубов (k-1)-го этапа дробления на кубы k-го этапа.

6. ЗАКЛЮЧЕНИЕ

Настоящая работа является продолжением работ [4], [6], [7], которые дают достаточно общие схемы решения для достаточно широкого класса задач многомерного интегрирования, что не отразилось в названии статьи. Это связано с тем, что в практическом смысле эти схемы работают лишь на достаточно узких классах задач. Например, если подынтегральная функция не будет являться “псевдомногочленом”, то проверка условия “погружения” на кубах промежуточных этапов дробления будет представлять из себя достаточно сложную вычислительную задачу, что может не позволить уменьшить размерность задачи на 1. Далее, в случае, когда подынтегральной функцией является гауссовская плотность вероятности, очень легко определяются величины ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{b} }_{Q}}$, ${{\bar {b}}_{Q}}$, используемые в основном алгоритме 2. Для более сложных подынтегральных функций вычисление этих величин может оказаться достаточно трудоемкой задачей. То же относится и к множеству интегрирования. Отметим, однако, что в случае многогранного множества весьма эффективно использовались оценки общей меры кубов (последнего этапа дробления) с небольшим числом “активных” ограничений (0, 1, 2) для получения оценок погрешности вычисления искомого многомерного интеграла (см. (3.12)). Аналогичные результаты при определенных условиях могут быть получены и в случае нескольких гладких ограничений. Большую роль в оценках погрешности в работах [6], [7] сыграли рекурсивные алгоритмы точного вычисления объема координатного куба (или параллелепипеда) с одним или двумя линейными ограничениями. Очевидно, что в случае с подынтегральным множеством, определяемым несколькими гладкими ограничениями, также будут использованы эти алгоритмы. Таким образом, работы [4], [6], [7] вместе с настоящей работой описывают практически реализуемые алгоритмы решения некоторых весьма важных для практического приложения задач, а также несут в себе идеи для некоторых более общих приложений.