Численное интегрирование

Как ищут интегралы?

Как многим известно, поиск значения определённого интеграла (или кратного, если речь идёт об интегрировании по поверхности) может происходить:

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

Метод суммы Римана;
Метод интерполирующего многочлена.

Метод суммы Римана известен всем, у кого была математика на первом курсе в универе или в школе был учитель, у которого нашлось время это рассказать в рамках программы старших классов.

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

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

В рамках данной статьи я покажу два способа нахождения интеграла по треугольнику от заранее заданной функции.

Локальные координаты треугольника и полиномы.

Пусть дан невырожденный треугольник с вершинами . Тогда любую точку данного треугольника можно представить через “локальные” координаты $\alpha, \beta$ :

$\left(\begin{matrix}x_p\\y_p\end{matrix}\right)=\left(\begin{matrix}a_x-c_x&b_x-c_x\\a_y-c_y&b_y-c_y\end{matrix}\right)\left(\begin{matrix}\alpha\\\beta\end{matrix}\right)+\left(\begin{matrix}c_x\\ c_y\end{matrix}\right).$

Данное представление возможно по той причине, что треугольник невырожденный, а следовательно векторы, построенные по его сторонам будут линейно независимыми, из чего следует, что

$\left|\begin{matrix}a_x-c_x&b_x-c_x\\a_y-c_y&b_y-c_y\end{matrix}\right|\ne0,$

что гарантирует единственность $\alpha, \beta$ , а ещё для любой точки треугольника будет верно, что

$\alpha\ge0, \beta\ge0, \alpha+\beta\le1.$

Таким образом при вычислении кратного интеграла можно сделать следующий переход:

$\iint\limits_{\triangle} f(x, y)dxdy=J\int\limits_{0}^{1}d\alpha\int\limits_{0}^{1-\alpha}f(x(\alpha, \beta), y(\alpha, \beta))d\beta,$

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

Теперь, если принять интерполяционный полином $P(\alpha, \beta)$ приблизительно равным интерполируемой функции $f(\alpha, \beta)$ , то можно с некоторым отклонением считать интегралы по треугольнику у функции и полинома равными:

$J\int\limits_{0}^{1}d\alpha\int\limits_{0}^{1-\alpha}f(x(\alpha, \beta), y(\alpha, \beta))d\beta=J\int\limits_{0}^{1}d\alpha\int\limits_{0}^{1-\alpha}P(\alpha, \beta)d\beta.$

Если учесть, что полином имеет вид

$P=\sum c_i \alpha^{n_i}\beta^{m_i}; n_i, m_i \in\mathbb{N},$

и тот факт, что интеграл от суммы равен сумме интегралов, то

$J\int\limits_{0}^{1}d\alpha\int\limits_{0}^{1-\alpha}P(\alpha, \beta)d\beta = J\sum c_i \left[\int\limits_{0}^{1}d\alpha\int\limits_{0}^{1-\alpha} \alpha^{n_i}\beta^{m_i}d\beta\right].$

Остаётся рассмотреть, чему равен интеграл от монома. Для начала нужно решить внутренний интеграл:

$\int\limits_{0}^{1}d\alpha\int\limits_{0}^{1-\alpha} \alpha^{n_i}\beta^{m_i}d\beta = \frac{1}{m_i+1} \int\limits_{0}^{1}\alpha^{n_i} (1-\alpha)^{m_i+1} d\alpha .$

Теперь, чтобы найти интеграл, нужно решить небольшую вспомогательную задачу: найти интеграл вида

$I_{n, m}=\int\limits_{0}^{1}\alpha^n (1-\alpha)^m d\alpha, n, m \in \mathbb{N}.$

Для этого нужно сделать интегрирование по частям:

$\int\limits_{0}^{1}\alpha^n (1-\alpha)^m d\alpha=\frac{1}{n+1}\int\limits_{0}^{1} (1-\alpha)^m d\alpha^{n+1} =\\ \frac{1}{n+1}\left[\underbrace{\left. \alpha^{n+1}(1-\alpha)^m \right|^{1}_0}_{0-0=0} - \int\limits_{0}^{1} \alpha^{n+1} d(1-\alpha)^m \right] = \frac{-1}{n+1}\int\limits_{0}^{1} \alpha^{n+1} d(1-\alpha)^m = \frac{m}{n+1}\int\limits_{0}^{1} \alpha^{n+1}(1-\alpha)^{m-1} d\alpha =\\=\ldots =\frac{m(m-1)(m-2)\ldots(m-(q-1))}{(n+1)(n+2)\ldots(n+q)}\int\limits_{0}^{1} \alpha^{n+q}(1-\alpha)^{m-q} d\alpha=\ldots=\\= \frac{m!}{(n+1)(n+2)\ldots(n+m)}\int\limits_{0}^{1} \alpha^{n+m}(1-\alpha)^{0}. d\alpha$

Далее если умножить полученное выражение на $inline $1=\frac{n!}{n!}, $inline$ то выходит следующее:

$\frac{m!}{(n+1)(n+2)\ldots(n+m)}\int\limits_{0}^{1} \alpha^{n+m}(1-\alpha)^{0} = \frac{m!n!}{(n+m)!}\int\limits_{0}^{1} \alpha^{n+m}(1-\alpha)^{0} = \frac{m!n!}{(n+m)!}\int\limits_{0}^{1} \alpha^{n+m} = \frac{m!n!}{(n+m+1)!} \left. \alpha^{n+m+1}\right|^1_0.$

А тогда выходит, что

$I_{n, m} = \frac{m!n!}{(n+m+1)!}.$

Тогда выходит так, что

$\frac{1}{m_i+1} \int\limits_{0}^{1}\alpha^{n_i} (1-\alpha)^{m_i+1} d\alpha = \frac{1}{m_i+1} I_{n_i, m_i+1} = \frac{1}{m_i+1} \frac{{n_i}!(m_i +1)!}{(n_i + m_i + 2)!} = \frac{{n_i}!{m_i}!}{(n_i + m_i + 2)!}.$

Интеграл полинома теперь находится как сумма:

$J\sum c_i \left[\int\limits_{0}^{1}d\alpha\int\limits_{0}^{1-\alpha} \alpha^{n_i}\beta^{m_i}d\beta\right] = J\sum c_i \frac{{n_i}!{m_i}!}{(n_i + m_i + 2)!}.$

Упрощение через фиксацию точек

Если на треугольнике зафиксировать точки и степень интерполяционного полинома, то выходит, что при вычислении интеграла понадобится лишь определять коэффициенты , а вот такие замечательные штуки, как дроби с факториалами становятся постоянными и их достаточно вычислить один раз.

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

$\sum c_i \frac{{n_i}!{m_i}!}{(n_i + m_i + 2)!} = \left(\begin{matrix} I_{n_1,m_1}&I_{n_2,m_2}&\ldots&I_{n_i, m_i}&\ldots&I_{n_q, m_q}\end{matrix}\right) \left(\begin{matrix}c_1\\c_2\\\ldots\\c_i\\\ldots\\c_q\end{matrix}\right).$

Коэффициенты находятся в зависимости от способа построения интерполяционного полинома. В рамках данной статьи строится интерполяционный полином Лагранжа, для которого коэффициенты находятся как решение СЛАУ вида

$\Lambda C = F,$

где это матрица-столбец, где на -ой строке стоит значение функции в точке , а $\Lambda$ это такая матрица, где элемент на позиции есть значение монома $x^{n_j} y^{m_j}$ в точке .

Если матрица $\Lambda$ квадратная и $det\Lambda\ne 0$ , то $С =\Lambda^{-1}F$ , а тогда можно переписать то «страшное» выражение с матрицами:

$\left(\begin{matrix} I_{n_1,m_1}&I_{n_2,m_2}&\ldots&I_{n_i, m_i}&\ldots&I_{n_q, m_q}\end{matrix}\right) \left(\begin{matrix}c_1\\c_2\\\ldots\\c_i\\\ldots\\c_q\end{matrix}\right) = \underbrace{\left(\begin{matrix} I_{n_1,m_1}&I_{n_2,m_2}&\ldots&I_{n_i, m_i}&\ldots&I_{n_q, m_q}\end{matrix}\right)\Lambda^{-1} }_{\Phi}\left(\begin{matrix}f_1\\f_2\\\ldots\\f_i\\\ldots\\f_q\end{matrix}\right) = \\ = \Phi F$

Описанное выражение с матрицами выше позволяет теперь записать следующее тождество

$J\sum c_i \frac{{n_i}!{m_i}!}{(n_i + m_i + 2)!} = J \Phi F,$

из чего следует

$J\int\limits_{0}^{1}d\alpha\int\limits_{0}^{1-\alpha}P(\alpha, \beta)d\beta = J \Phi F.$

Таким образом достаточно вычислить только один раз матрицу-строку $\Phi$ и для каждой новой функции или треугольника вычислять значения якобиана и значений в «локальных» координатах(так как изначально полином строился в локальных координатах).

Эксперимент

Подбор разных полиномов и локальных точек.

Линейный полином

Линейный полином имеет вид

$P = c_1 + c_2 \alpha + c_3 \beta.$

Для данного полинома автором данной статьи предлагается брать такие «локальные» координаты, которые соответствуют вершинам треугольника:

№	$\alpha$	$\beta$
1
2
3

В результате проделанных манипуляций, получается, что

$\Phi = \left(\begin{matrix} \frac{1}{6}&\frac{1}{6}&\frac{1}{6} \end{matrix}\right).$

Тогда интеграл функции по треугольнику, если в качестве интерполяционного полинома брать линейный полином и интерполировать по вершинам, будет равен(приблизительно) как шестая часть от суммы значений исходной функции в вершинах, умноженной на якобиан:

$\iint\limits_{\triangle} f(x,y)dxdy\approx J\frac{f(a_x, a_y)+f(b_x, b_y)+f(c_x, c_y)}{6}$

Полином второй степени

Полином второй степени имеет вид

$P = c_1 + c_2 \alpha + c_3 \beta + с_4 \alpha^2 + c_5 \alpha\beta + c_6 \beta^2 .$

Используемые “локальные” точки:

№	$\alpha$	$\beta$
1
2
3
4	$\frac{1}{2}$	$\frac{1}{2}$
5		$\frac{1}{2}$
6	$\frac{1}{2}$

Полученная матрица $\Phi$ имеет следующий вид:

$\Phi = \left(\begin{matrix} 0 & 0 & 0 & \frac{1}{6}&\frac{1}{6}&\frac{1}{6} \end{matrix}\right).$

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

Отдельный полином

Рассмотрим полином, который представляется из тех и только тех мономов, для которых степени соответствующих переменных принимают натуральные значения от 0 до 3:

$P = \sum{ [c_k m_k]} ; m_k = x^{i_k} y^{j_k}; i_k, i_k\in\overline{0,3}.$

Такой полином имеет вид

$P = c_1 + c_2 \alpha + c_3 \alpha^2 + c_4 \alpha^3 + c_5 \beta + c_6 \alpha\beta +c_7 \alpha^2 \beta + c_8 \alpha^3 \beta + c_9 \beta^2 + c_{10} \alpha \beta^2 + c_{11} \alpha^2 \beta^2 + c_{12} \alpha^3 \beta^2 + c_{13} \beta^3 + c_{14} \alpha \beta^3 + c_{15} \alpha^2 \beta^3 + c_{16} \alpha^3 \beta^3.$

Локальные точки:

№	$\alpha$	$\beta$
1	$\frac{2}{3}$	$\frac{1}{6}$
2	$\frac{1}{6}$	$\frac{2}{3}$
3	$\frac{1}{6}$	$\frac{1}{6}$
4	$\frac{1}{3}$	$\frac{1}{3}$
5	$\frac{1}{2}$
6		$\frac{1}{2}$
7	$\frac{1}{2}$	$\frac{1}{2}$
8	$\frac{1}{4}$
9		$\frac{1}{4}$
10
11		$\frac{1}{4}$
12
13	$\frac{1}{4}$
14	$\frac{1}{2}$	$\frac{1}{4}$
15	$\frac{1}{4}$	$\frac{1}{2}$
16	$\frac{1}{4}$	$\frac{1}{4}$

После расчётов получается, что

$\Phi = \left(\begin{matrix} \frac{1}{2} & \frac{1}{6} & \frac{1}{12} & \frac{1}{20}& \frac{1}{6}& \frac{1}{24} &\frac{1}{60}&\frac{1}{120}& \frac{1}{12}& \frac{1}{60} & \frac{1}{180} & \frac{1}{420} & \frac{1}{20} & \frac{1}{120} & \frac{1}{420} & \frac{1}{1120} \end{matrix}\right)$

Программная реализация

Для реализации задуманного способа интегрирования был написан код, в рамках которого написаны функции integrate_p3(реализация интегрирования с использованием линейного полинома), integrate_p6(реализация интегрирования с использованием полинома второй степени) и integrate_p16(реализация интегрирования с использованием отдельного полинома). Если есть желание, можно ознакомиться с кодом и даже потыкать самостоятельно, чтобы увидеть, какой результат выйдет при тех или иных условиях.

Вообще, в результате экспериментов было обнаружено, что чем меньше площадь у треугольника, тем больше расчёт через выведенные формулы близок по значению с тем, что даёт метод суммы Римана, а иногда даже совпадает с разными аналитическими решениями.

А зачем это вообще надо?

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

А что дальше?

Если произойдёт чудо и эта статья будет выпущена и хоть кому-то понравится, то я продолжу эту тематику. Кроме этого у меня была мысль написать небольшой цикл по некоторым иным темам, но до таких дней, когда у меня до этого дойдут руки, ещё надо дожить.

ссылка на оригинал статьи https://habr.com/ru/articles/1057290/