Численное Интегрирование Многомерных Функций: От Кубатурных Формул до Методов Монте-Карло

#ЧисленноеИнтегрирование #КубатурныеФормулы #МногомерныеИнтегралы #МетодыМонтеКарло #ОшибкаИнтегрирования #ВычислительнаяМатематика

Введение в Многомерное Интегрирование

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

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

Кубатурные Формулы для Прямоугольных Областей

Рассмотрим интеграл по прямоугольной области $[\alpha, \beta] \times [a, b]$:

$$I(f) = \int_{\alpha}^{\beta} \int_{a}^{b} f(x, y) \, dx \, dy$$

Для начала предположим, что стороны этого прямоугольника малы, то есть $(\beta - \alpha)$ и $(b - a)$ малы. Затем мы можем разбить большую область на множество таких маленьких прямоугольников, применяя ту же идею, что и в одномерном случае. Фактически, у нас есть два малых параметра — длины сторон прямоугольника, по которым мы будем раскладывать функцию.

Мы будем строить аналог формулы средних прямоугольников. Для этого обозначим координаты середины прямоугольника как $x_c = \frac{\alpha + \beta}{2}$ и $y_c = \frac{a + b}{2}$. Площадь прямоугольника обозначим как $S = (\beta - \alpha)(b - a)$.

В формуле средних прямоугольников мы заменяем подынтегральную функцию её значением в середине отрезка. Здесь мы можем заменить функцию $f(x, y)$ её значением в центральной точке прямоугольника $f(x_c, y_c)$. Тогда простейшая кубатурная формула (для размерности больше единицы такие конструкции называются кубатурными формулами, а не квадратурными) будет выглядеть так:

$$I(f) \approx S \cdot f(x_c, y_c)$$

Эта формула точна для функций вида $f(x, y) = Px + Qy + D$, то есть для линейных функций двух переменных.

Вывод Ошибки Кубатурной Формулы

Для оценки ошибки воспользуемся разложением функции $f(x, y)$ в ряд Тейлора относительно центральной точки $(x_c, y_c)$:

$$f(x, y) = f(x_c, y_c) + \frac{\partial f}{\partial x}(x_c, y_c)(x - x_c) + \frac{\partial f}{\partial y}(x_c, y_c)(y - y_c) + \frac{1}{2!} \left( \frac{\partial^2 f}{\partial x^2}(x_c, y_c)(x - x_c)^2 + 2 \frac{\partial^2 f}{\partial x \partial y}(x_c, y_c)(x - x_c)(y - y_c) + \frac{\partial^2 f}{\partial y^2}(x_c, y_c)(y - y_c)^2 \right) + O(h^3)$$

где $h$ — малый параметр, характеризующий размеры прямоугольника.

Интегрируем это разложение по нашей области. Первый член: $\int_{\alpha}^{\beta} \int_{a}^{b} f(x_c, y_c) \, dx \, dy = S \cdot f(x_c, y_c)$. Это и есть наша кубатурная формула. Второй и третий члены (с первыми производными): $\int_{\alpha}^{\beta} (x - x_c) \, dx = 0$ и $\int_{a}^{b} (y - y_c) \, dy = 0$, так как интегрирование ведется по симметричным относительно середины отрезкам. Следовательно, эти члены обнуляются. Член с $\frac{\partial^2 f}{\partial x \partial y}(x_c, y_c)(x - x_c)(y - y_c)$ также обнуляется из-за симметрии. Остаются члены со вторыми производными по одной переменной:

$$\int_{\alpha}^{\beta} (x - x_c)^2 \, dx = \frac{(\beta - \alpha)^3}{12}$$ $$\int_{a}^{b} (y - y_c)^2 \, dy = \frac{(b - a)^3}{12}$$

Тогда главный член погрешности будет:

$$R(f) = \frac{1}{24} \left( (\beta - \alpha)^3 \frac{\partial^2 f}{\partial x^2}(x_c, y_c) (b - a) + (b - a)^3 \frac{\partial^2 f}{\partial y^2}(x_c, y_c) (\beta - \alpha) \right)$$ $$R(f) = \frac{S}{24} \left( (\beta - \alpha)^2 \frac{\partial^2 f}{\partial x^2}(x_c, y_c) + (b - a)^2 \frac{\partial^2 f}{\partial y^2}(x_c, y_c) \right)$$

Обозначим $h_1 = \beta - \alpha$ и $h_2 = b - a$. Тогда погрешность имеет порядок $O(S(h_1^2 + h_2^2))$.

Композитная Кубатурная Формула

Для интегрирования по большой области, мы разбиваем её на $N_1 \times N_2$ маленьких прямоугольников. Пусть $h_1 = \frac{\beta - \alpha}{N_1}$ и $h_2 = \frac{b - a}{N_2}$. Общая погрешность будет суммой погрешностей по каждому прямоугольнику. Если $N_1$ и $N_2$ — количество разбиений по каждой координате, то суммарная погрешность $R_{total}$ будет:

$$R_{total} = \sum_{i,j} R_{ij} \approx \sum_{i,j} \frac{S_{ij}}{24} \left( h_1^2 \frac{\partial^2 f}{\partial x^2}(x_{c,ij}, y_{c,ij}) + h_2^2 \frac{\partial^2 f}{\partial y^2}(x_{c,ij}, y_{c,ij}) \right)$$

Если предположить, что вторые производные ограничены, то $R_{total} \approx O(N_1 N_2 \cdot (h_1^2 + h_2^2))$. Подставляя $h_1 = \frac{\beta - \alpha}{N_1}$ и $h_2 = \frac{b - a}{N_2}$, получаем:

$$R_{total} \approx O\left( N_1 N_2 \left( \frac{(\beta - \alpha)^2}{N_1^2} + \frac{(b - a)^2}{N_2^2} \right) \right) = O\left( \frac{N_2}{N_1} (\beta - \alpha)^2 + \frac{N_1}{N_2} (b - a)^2 \right)$$

Если $N_1 \approx N_2 \approx N$, то $R_{total} \approx O(h_1^2 + h_2^2)$. Общее число узлов $N_{total} = N_1 N_2$. Тогда $h_1 \sim 1/N_1$ и $h_2 \sim 1/N_2$. Если $N_1 = N_2 = \sqrt{N_{total}}$, то $h_1 \sim 1/\sqrt{N_{total}}$ и $h_2 \sim 1/\sqrt{N_{total}}$. Тогда погрешность $R_{total} \approx O(1/N_{total})$. Это означает, что порядок сходимости по общему числу точек $N_{total}$ равен $1$. В одномерном случае для формулы средних прямоугольников погрешность была $O(1/N^2)$ по длине шага $h$, или $O(1/N^2)$ по числу точек $N$. Здесь же, при том же порядке по шагу, порядок по числу точек ухудшается. Для двумерного случая это $O(1/N_{total})$, а для $k$-мерного — $O(1/N_{total}^{1/k})$. Это существенное ухудшение эффективности.

Последовательное Интегрирование

Другой подход — последовательное интегрирование, когда многомерный интеграл сводится к повторному применению одномерных интегралов.

$$I(f) = \int_{C}^{D} \left( \int_{A}^{B} f(x, y) \, dx \right) \, dy$$

Обозначим внутренний интеграл как $F(y) = \int_{A}^{B} f(x, y) \, dx$. Мы можем вычислить $F(y)$ для каждой точки $y_j$ с помощью одномерной квадратурной формулы (например, формулы средних прямоугольников или трапеций):

$$F(y_j) \approx \sum_{i} c_i f(x_i, y_j)$$

Затем внешний интеграл $\int_{C}^{D} F(y) \, dy$ также вычисляется с помощью одномерной квадратурной формулы:

$$I(f) \approx \sum_{j} d_j F(y_j) \approx \sum_{j} d_j \sum_{i} c_i f(x_i, y_j)$$

Если использовать формулу средних прямоугольников для обоих интегралов, мы получим ту же формулу, что и при прямом применении кубатурной формулы средних прямоугольников. Если использовать формулу трапеций, то узлы будут иметь веса, зависящие от их положения: $h_x h_y$ для внутренних точек, $1/2 h_x h_y$ для точек на границе (но не в углах), и $1/4 h_x h_y$ для угловых точек. Порядок погрешности при этом будет таким же, как и для формулы средних прямоугольников.

Интегрирование по Непрямоугольным Областям: Триангуляция

Если область интегрирования не является прямоугольной, разбиение на прямоугольники становится неудобным, особенно вблизи границ. В таких случаях гораздо более гибким решением является триангуляция области, то есть разбиение её на треугольники. Для трехмерных областей используются тетраэдры (пирамидки).

Современные алгоритмы умеют автоматически выполнять триангуляцию, создавая "хорошие" треугольники (например, избегая слишком тонких или вырожденных треугольников). После триангуляции функцию можно интерполировать на каждом треугольнике кусочно-линейными или кусочно-квадратичными функциями, а затем складывать результаты интегрирования по этим элементам.

Проклятие Размерности и Методы Монте-Карло

Проблема численного интегрирования становится особенно острой в высоких размерностях. В 1950-60-х годах, например, при расчетах ядерных реакторов, возникала необходимость интегрировать функции с размерностью до десятка, при этом вычислительные мощности были крайне ограничены.

Пусть у нас есть $k$-мерная область, и мы разбиваем каждую сторону на $N_i$ частей. Общее число узлов сетки $N_{total} = N_1 \cdot N_2 \cdot \ldots \cdot N_k$. Если все $N_i = N$, то $N_{total} = N^k$. Предположим, что погрешность нашей кубатурной формулы имеет порядок $O(h^p)$, где $h$ — характерный шаг сетки. В терминах общего числа узлов $N_{total}$ погрешность будет $O(N_{total}^{-p/k})$. Если мы хотим достичь точности $\varepsilon$, то $N_{total}^{-p/k} \sim \varepsilon$, откуда $N_{total} \sim \varepsilon^{-k/p}$.

Рассмотрим пример:

• Требуемая точность $\varepsilon = 0.01$ (два знака после запятой).
• Размерность $k = 10$.
• Порядок формулы $p = 2$ (например, формула средних прямоугольников).

Тогда $N_{total} \sim (0.01)^{-10/2} = (0.01)^{-5} = (10^{-2})^{-5} = 10^{10}$. Это означает, что для достижения такой точности нам потребуется $10^{10}$ точек, что соответствует $10^{10}$ операций умножения и сложения. Это совершенно непосильная задача для современных компьютеров, не говоря уже о вычислительных мощностях 60-х годов. Это явление известно как "проклятие размерности".

Для борьбы с этой проблемой были разработаны вероятностные методы, известные как методы Монте-Карло.

Методы Монте-Карло для Интегрирования

Основная идея методов Монте-Карло заключается в следующем: Пусть у нас есть интеграл по области $\Omega$, мера которой равна единице (если мера не равна единице, можно нормировать).

$$I(f) = \int_{\Omega} f(x) \, d\mu(x)$$

Мы выбираем $N$ случайных точек $T_1, T_2, \ldots, T_N$ из области $\Omega$. Эти точки должны быть попарно независимыми и равномерно распределенными. Технически это реализуется с помощью генераторов псевдослучайных чисел.

Тогда интеграл аппроксимируется следующей кубатурной формулой:

$$I(f) \approx \frac{1}{N} \sum_{i=1}^{N} f(T_i)$$

Погрешность этой формулы оценивается с помощью неравенства Чебышева. Утверждается, что с вероятностью $1 - \theta$ (для любого $\theta \in (0, 1)$) погрешность будет порядка $O\left(\frac{1}{\sqrt{N}}\right)$. Точнее, она будет меньше или равна $\frac{C}{\sqrt{N}}$ с вероятностью $1-\theta$. Важно, что порядок погрешности $O(1/\sqrt{N})$ не зависит от размерности $k$. Это ключевое преимущество методов Монте-Карло в высоких размерностях.

Сравним с предыдущим примером:

• Требуемая точность $\varepsilon = 0.01$.
• $1/\sqrt{N} \sim 0.01 \Rightarrow N \sim (1/0.01)^2 = 100^2 = 10^4$. Для достижения той же точности в десятимерном пространстве нам потребуется всего $10^4$ точек, вместо $10^{10}$ точек, необходимых для детерминированных кубатурных формул. Это делает методы Монте-Карло незаменимыми для интегрирования в высоких размерностях.

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

На этом мы завершаем рассмотрение численного интегрирования.

---

Generated by AI-powered TranscribeLecture.com • 14.03.2026