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

2026-07-09 в 7:24, admin, рубрики: интегрирование, интерполяционные полиномы, метод конечных элементов, прикладная математика, численное интегрирование

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

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

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

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

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

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

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

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

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

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

$left(begin{matrix}x_p\y_pend{matrix}right)=left(begin{matrix}a_x-c_x&b_x-c_x\a_y-c_y&b_y-c_yend{matrix}right)left(begin{matrix}alpha\betaend{matrix}right)+left(begin{matrix}c_x\ c_yend{matrix}right).$

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

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

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

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

$iintlimits_{triangle} f(x, y)dxdy=Jintlimits_{0}^{1}dalphaintlimits_{0}^{1-alpha}f(x(alpha, beta), y(alpha, beta))dbeta,$

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

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

$Jintlimits_{0}^{1}dalphaintlimits_{0}^{1-alpha}f(x(alpha, beta), y(alpha, beta))dbeta=Jintlimits_{0}^{1}dalphaintlimits_{0}^{1-alpha}P(alpha, beta)dbeta.$

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

$P=sum c_i alpha^{n_i}beta^{m_i}; n_i, m_i inmathbb{N},$

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

$Jintlimits_{0}^{1}dalphaintlimits_{0}^{1-alpha}P(alpha, beta)dbeta=Jsum c_i left[intlimits_{0}^{1}dalphaintlimits_{0}^{1-alpha} alpha^{n_i}beta^{m_i}dbetaright].$

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

$intlimits_{0}^{1}dalphaintlimits_{0}^{1-alpha} alpha^{n_i}beta^{m_i}dbeta=frac{1}{m_i+1} intlimits_{0}^{1}alpha^{n_i} (1-alpha)^{m_i+1} dalpha .$

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

$I_{n, m}=intlimits_{0}^{1}alpha^n (1-alpha)^m dalpha, n, m in mathbb{N}.$

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

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

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

$frac{m!}{(n+1)(n+2)ldots(n+m)}intlimits_{0}^{1} alpha^{n+m}(1-alpha)^{0}=frac{m!n!}{(n+m)!}intlimits_{0}^{1} alpha^{n+m}(1-alpha)^{0}=frac{m!n!}{(n+m)!}intlimits_{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} intlimits_{0}^{1}alpha^{n_i} (1-alpha)^{m_i+1} dalpha=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)!}.$

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

$Jsum c_i left[intlimits_{0}^{1}dalphaintlimits_{0}^{1-alpha} alpha^{n_i}beta^{m_i}dbetaright]=Jsum 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_qend{matrix}right).$

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

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

Если матрица квадратная и , то $С=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_qend{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_qend{matrix}right)=\=Phi F$

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

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

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

$Jintlimits_{0}^{1}dalphaintlimits_{0}^{1-alpha}P(alpha, beta)dbeta=J Phi F.$

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

Эксперимент

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

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

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

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

№
1
2
3

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

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

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

$iintlimits_{triangle} f(x,y)dxdyapprox Jfrac{f(a_x, a_y)+f(b_x, b_y)+f(c_x, c_y)}{6}$

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

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

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

№
1
2
3
4	$frac{1}{2}$	$frac{1}{2}$
5		$frac{1}{2}$
6	$frac{1}{2}$

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

$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_kinoverline{0,3}.$

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

$P=c_1 + c_2 alpha + c_3 alpha^2 + c_4 alpha^3 + c_5 beta + c_6 alphabeta +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.$

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

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

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

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

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

А что дальше?

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

Автор: Dahnny_142

Источник

Информация

Комментарии

Рекомендуем