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

Численное интегрирование (историческое название: (численная) квадратура) — вычисление значения определённого интеграла (как правило, приближённое). Под численным интегрированием понимают набор численных методов для нахождения значения определённого интеграла.

Численное интегрирование применяется, когда:

1. Сама подынтегральная функция не задана аналитически. Например, она представлена в виде таблицы (массива) значений в узлах некоторой расчётной сетки.
2. Аналитическое представление подынтегральной функции известно, но её первообразная не выражается через аналитические функции. Например, ${\displaystyle f(x)=\exp(-x^{2})}$.

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

Одномерный случай

Одномерный определённый интеграл как площадь криволинейной трапеции под графиком

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

${\displaystyle I\approx \sum _{i=1}^{n}w_{i}\,f(x_{i}),}$

где ${\displaystyle n}$ — число точек, в которых вычисляется значение подынтегральной функции. Точки ${\displaystyle x_{i}}$ называются узлами метода, числа ${\displaystyle w_{i}}$ — весами узлов. При замене подынтегральной функции на полином нулевой, первой и второй степени получаются соответственно методы прямоугольников, трапеций и парабол (Симпсона). Часто формулы для оценки значения интеграла называют квадратурными формулами.

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

${\displaystyle \int \limits _{a}^{b}f(x)\,dx=\sum _{i=0}^{n}H_{i}\,f(x_{i})+r_{n}(f),}$

где числа ${\displaystyle H_{i}}$ называются коэффициентами Котеса и вычисляются как интегралы от соответствующих многочленов, стоящих в исходном интерполяционном многочлене для подынтегральной функции при значении функции в узле ${\displaystyle x_{i}=a+ih}$ (${\displaystyle h=(b-a)/n}$ — шаг сетки; ${\displaystyle n}$ — число узлов сетки, а индекс узлов ${\displaystyle i=0\ldots n}$). Слагаемое ${\displaystyle r_{n}(f)}$ — погрешность метода, которая может быть найдена разными способами. Для нечетных ${\displaystyle n\geqslant 1}$ погрешность может быть найдена интегрированием погрешности интерполяционного полинома подынтегральной функции.

Частными случаями формул Котеса являются: формулы прямоугольников (${\displaystyle n=0}$), формулы трапеций (${\displaystyle n=1}$), формула Симпсона (${\displaystyle n=2}$), формула Ньютона (${\displaystyle n=3}$) и т. д.

Метод прямоугольников

Основная статья: Метод прямоугольников

Пусть требуется определить значение интеграла функции на отрезке ${\displaystyle \left[{a},{b}\right]}$. Этот отрезок делится точками ${\displaystyle x_{0},x_{1},\ldots ,x_{n-1},x_{n}}$ на ${\displaystyle n}$ равных отрезков длиной ${\displaystyle \Delta {x}={\frac {b-a}{n}}.}$ Обозначим через ${\displaystyle y_{0},y_{1},\ldots ,y_{n-1},y_{n}}$ значение функции ${\displaystyle f\left(x\right)}$ в точках ${\displaystyle x_{0},x_{1},\ldots ,x_{n-1},x_{n}.}$ Далее составляем суммы ${\displaystyle y_{0}\,\Delta {x}+y_{1}\,\Delta {x}+\ldots +y_{n-1}\,\Delta {x}.}$ Каждая из сумм — интегральная сумма для ${\displaystyle f\left(x\right)}$ на ${\displaystyle \left[{a},{b}\right]}$ и поэтому приближённо выражает интеграл

${\displaystyle \int \limits _{a}^{b}f(x)\,dx\approx {\frac {b-a}{n}}(y_{0}+y_{1}+\ldots +y_{n-1}).}$

Если заданная функция — положительная и возрастающая, то эта формула выражает площадь ступенчатой фигуры, составленной из «входящих» прямоугольников, также называемая формулой левых прямоугольников, а формула

${\displaystyle \int \limits _{a}^{b}f(x)\,dx\approx {\frac {b-a}{n}}(y_{1}+y_{2}+\ldots +y_{n})}$

выражает площадь ступенчатой фигуры, состоящей из «выходящих» прямоугольников, также называемая формулой правых прямоугольников. Чем меньше длина отрезков, на которые делится отрезок ${\displaystyle \left[{a},{b}\right]}$, тем точнее значение, вычисляемое по этой формуле, искомого интеграла.

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

${\displaystyle \int \limits _{a}^{b}f(x)\,dx\approx h\sum _{i=1}^{n}f\left(x_{i-1}+{\frac {h}{2}}\right)=h\sum _{i=1}^{n}f\left(x_{i}-{\frac {h}{2}}\right),}$

где ${\displaystyle h={\frac {b-a}{n}}}$

Учитывая априорно бо́льшую точность последней формулы при том же объёме и характере вычислений её называют формулой прямоугольников

Метод трапеций

Основная статья: Метод трапеций

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

Площадь трапеции на каждом отрезке:

${\displaystyle I_{i}\approx {\frac {f(x_{i-1})+f(x_{i})}{2}}(x_{i}-x_{i-1}).}$

Погрешность аппроксимации на каждом отрезке:

${\displaystyle \left|R_{i}\right|\leqslant {\frac {\left(b-a\right)^{3}}{12n^{2}}}M_{2,i},}$ где ${\displaystyle M_{2,i}=\max _{x\in [x_{i-1},\,x_{i}]}{\left|f''(x)\right|}.}$

Полная формула трапеций в случае деления всего промежутка интегрирования на отрезки одинаковой длины ${\displaystyle h}$:

${\displaystyle I\approx h\left({\frac {f(x_{0})+f(x_{n})}{2}}+\sum _{i=1}^{n-1}f(x_{i})\right),}$ где ${\displaystyle h={\frac {b-a}{n}}.}$

Погрешность формулы трапеций:

${\displaystyle \left|R\right|\leqslant {\frac {\left(b-a\right)^{3}}{12n^{2}}}M_{2},}$ где ${\displaystyle M_{2}=\max _{x\in [a,\,b]}{\left|f''(x)\right|}.}$

Метод парабол (метод Симпсона)

Основная статья: Метод Симпсона

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

${\displaystyle I\approx {\frac {b-a}{6}}\left(f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right)}$.

Если разбить интервал интегрирования на ${\displaystyle 2N}$ равных частей, то имеем

${\displaystyle I\approx {\frac {b-a}{6N}}\left(f_{0}+4\left(f_{1}+f_{3}+\ldots +f_{2N-1}\right)+2\left(f_{2}+f_{4}+\ldots +f_{2N-2}\right)+f_{2N}\right),}$

где ${\displaystyle f_{i}=f\left(a+{\frac {(b-a)i}{2N}}\right)}$.

Увеличение точности

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

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

При стремлении количества разбиений к бесконечности оценка интеграла стремится к его истинному значению для аналитических функций для любого численного метода.

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

Метод Гаусса

Основная статья: Метод Гаусса (численное интегрирование)

Описанные выше методы используют фиксированные точки отрезка (концы и середину) и имеют низкий порядок точности (0 — методы правых и левых прямоугольников, 1 — методы средних прямоугольников и трапеций, 3 — метод парабол (Симпсона)). Если мы можем выбирать точки, в которых мы вычисляем значения функции ${\displaystyle f(x)}$, то можно при том же количестве вычислений подынтегральной функции получить методы более высокого порядка точности. Так, для двух (как в методе трапеций) вычислений значений подынтегральной функции можно получить метод уже не второго, а третьего порядка точности:

${\displaystyle I\approx {\frac {b-a}{2}}\left(f\left({\frac {a+b}{2}}-{\frac {b-a}{2{\sqrt {3}}}}\right)+f\left({\frac {a+b}{2}}+{\frac {b-a}{2{\sqrt {3}}}}\right)\right)}$.

В общем случае, используя ${\displaystyle n}$ точек, по формуле ${\displaystyle I\approx \sum _{i=1}^{n}a_{i}\,f(x_{i})}$ можно получить метод с порядком точности ${\displaystyle 2n-1}$, т. е. получаются точные значения для полиномов степени не выше ${\displaystyle 2n-1}$.

Значения узлов ${\displaystyle x_{i}}$ метода Гаусса по ${\displaystyle n}$ точкам являются корнями полинома Лежандра степени ${\displaystyle n}$. Значения весов вычисляются по формуле ${\displaystyle a_{i}={\frac {2}{(1-x_{i}^{2})\,[P_{n}'(x_{i})]^{2}}}}$, где ${\displaystyle P_{n}'}$ - первая производная полинома Лежандра.

Для ${\displaystyle n=3}$ узлы и веса имеют следующие значения : ${\displaystyle x_{1,3}=\pm {\sqrt {0.6}},x_{2}=0,}$ веса : ${\displaystyle a_{1,3}={\frac {5}{9}},a_{2}={\frac {8}{9}}}$

(полином определен на отрезке ${\displaystyle [-1,1]}$).

Наиболее известен метод Гаусса по пяти точкам.

Метод Гаусса — Кронрода

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

${\displaystyle I\approx \sum _{i=1}^{n}a_{i}\,f(x_{i})+\sum _{i=1}^{n+1}b_{i}\,f(y_{i})}$,

где ${\displaystyle x_{i}}$ — узлы метода Гаусса по ${\displaystyle n}$ точкам, а ${\displaystyle 3n+2}$ параметров ${\displaystyle a_{i}}$, ${\displaystyle b_{i}}$, ${\displaystyle y_{i}}$ подобраны таким образом, чтобы порядок точности метода был равен ${\displaystyle 3n+1}$.

Тогда для оценки погрешности можно использовать эмпирическую формулу:

${\displaystyle \Delta =\left(200|I-I_{G}|\right)^{1.5}}$,

где ${\displaystyle I_{G}}$ — приближённое значение интеграла, полученное методом Гаусса по ${\displaystyle n}$ точкам. Библиотеки gsl и SLATEC для вычисления определённых интегралов содержат подпрограммы, использующие метод Гаусса — Кронрода по 15, 21, 31, 41, 51 и 61 точкам. Библиотека ALGLIB использует метод Гаусса — Кронрода по 15 точкам.

Метод Чебышёва

Метод Чебышева (или как его иногда называют Гаусса — Чебышева) является одним из представителей методов наивысшей алгебраической точности Гаусса. Его отличительной особенностью является наличие у подынтегральной функции множителя ${\displaystyle \left({\frac {1}{\sqrt {1-x^{2}}}}\right)}$. Суть заключается в следующем:

${\displaystyle \int \limits _{-1}^{1}f(x)\left({\frac {1}{\sqrt {1-x^{2}}}}\right)\,dx=\sum _{i=1}^{N}C_{i}f(x_{i})}$,

где ${\displaystyle C_{i}=\pi /N}$, ${\displaystyle x_{i}=\cos \left({\frac {(2i-1)\pi }{2N}}\right)}$, ${\displaystyle N}$ — количество узлов метода.

Метод Гаусса-Лагера

Метод Гаусса-Эрмита

Интегрирование при бесконечных пределах

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

См. в том числе Метод Самокиша.

Методы Монте-Карло

Основная статья: Метод Монте-Карло § Интегрирование методом Монте-Карло

Рисунок 3. Численное интегрирование функции методом Монте-Карло

Для определения площади под графиком функции можно использовать следующий стохастический алгоритм:

• ограничим функцию прямоугольником (n-мерным параллелепипедом в случае многих измерений), площадь которого ${\displaystyle S_{par}}$ можно легко вычислить;
• «набросаем» в этот прямоугольник (параллелепипед) некоторое количество точек (${\displaystyle N}$ штук), координаты которых будем выбирать случайным образом;
• определим число точек (${\displaystyle K}$ штук), которые попадут под график функции;
• площадь области, ограниченной функцией и осями координат, ${\displaystyle S}$ даётся выражением ${\displaystyle S=S_{par}{\frac {K}{N}}}$;

Этот алгоритм требует определения экстремумов функции на интервале и не использует вычисленное точное значение функции ${\displaystyle f(x)}$ кроме как в сравнении, вследствие чего непригоден для практики. Приведённые в основной статье варианты метода Монте-Карло избавлены от этих недостатков.

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

Методы Рунге — Кутты

Основная статья: Метод Рунге — Кутты

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

Метод сплайнов

Многомерный случай

Пример узлов интегрирования на тетраэдре

В небольших размерностях можно так же применять квадратурные формулы, основанные на интерполяционных многочленах. Интегрирование производится аналогично одномерному интегрированию. Для больших размерностей эти методы становятся неприемлемыми из-за быстрого возрастания числа точек сетки и/или сложной границы области. В этом случае применяется метод Монте-Карло. Генерируются случайные точки в нашей области и усредняются значения функции в них. Так же можно использовать смешанный подход — разбить область на несколько частей, в каждой из которых (или только в тех, где интеграл посчитать не удаётся из-за сложной границы) применить метод Монте-Карло.

Примеры реализации

Ниже приведена реализация на Python 3 метода средних прямоугольников, метода средних трапеций, метода Симпсона и метода Монте-Карло.

import math, random
from numpy import arange

def get_i():
    return math.e ** 1 - math.e ** 0

def method_of_rectangles(func, min_lim, max_lim, delta):
    def integrate(func, min_lim, max_lim, n):
        integral = 0.0
        step = (max_lim - min_lim) / n
        for x in arange(min_lim, max_lim-step, step):
            integral += step * func(x + step / 2)
        return integral

    d, n = 1, 1
    while abs(d) > delta:
        d = (integrate(func, min_lim, max_lim, n * 2) - integrate(func, min_lim, max_lim, n)) / 3
        n *= 2

    a = abs(integrate(func, min_lim, max_lim, n))
    b = abs(integrate(func, min_lim, max_lim, n)) + d
    if a > b:
        a, b = b, a
    print('Rectangles:')
    print('\t%s\t%s\t%s' % (n, a, b))

def trapezium_method(func, min_lim, max_lim, delta):
    def integrate(func, min_lim, max_lim, n):
        integral = 0.0
        step = (max_lim - min_lim) / n
        for x in arange(min_lim, max_lim-step, step):
            integral += step*(func(x) + func(x + step)) / 2
        return integral

    d, n = 1, 1
    while abs(d) > delta:
        d = (integrate(func, min_lim, max_lim, n * 2) - integrate(func, min_lim, max_lim, n)) / 3
        n *= 2

    a = abs(integrate(func, min_lim, max_lim, n))
    b = abs(integrate(func, min_lim, max_lim, n)) + d
    if a > b:
        a, b = b, a
    print('Trapezium:')
    print('\t%s\t%s\t%s' % (n, a, b))

def simpson_method(func, min_lim, max_lim, delta):
    def integrate(func, min_lim, max_lim, n):
        integral = 0.0
        step = (max_lim - min_lim) / n
        for x in arange(min_lim + step / 2, max_lim - step / 2, step):
            integral += step / 6 * (func(x - step / 2) + 4 * func(x) + func(x + step / 2))
        return integral

    d, n = 1, 1
    while abs(d) > delta:
        d = (integrate(func, min_lim, max_lim, n * 2) - integrate(func, min_lim, max_lim, n)) / 15
        n *= 2

    a = abs(integrate(func, min_lim, max_lim, n))
    b = abs(integrate(func, min_lim, max_lim, n)) + d
    if a > b:
        a, b = b, a
    print('Simpson:')
    print('\t%s\t%s\t%s' % (n, a, b))

def monte_karlo_method(func, n):
    in_d, out_d = 0., 0.
    for i in arange(n):
        x, y = random.uniform(0, 1), random.uniform(0, 3)
        if y < func(x): in_d += 1

    print('M-K:')
    print('\t%s\t%s' % (n, abs(in_d/n * 3)))

method_of_rectangles(lambda x: math.e ** x, 0.0, 1.0, 0.001)
trapezium_method(lambda x: math.e ** x, 0.0, 1.0, 0.001)
simpson_method(lambda x: math.e ** x, 0.0, 1.0, 0.001)
monte_karlo_method(lambda x: math.e ** x, 100)
print('True value:\n\t%s' % get_i())

Литература

• Каханер Д., Моулер К., Нэш С. Численные методы и программное обеспечение (пер. с англ.).. — Изд. второе, стереотип.. — М.: Мир, 2001. — 575 с. — ISBN 5-03-003392-0.
• Самарский А. А., Гулин А. В. Численные методы: Учеб. пособие для вузов. — М.: Наука. Гл. ред. физ-мат. лит., 1989. — 432 с. — ISBN 5-02-013996-3.
• Пискунов Н. С. Дифференциальное и интегральное исчисления : Учеб. пособие для вузов: В 2 т. — 13-е изд.. — М.: Наука. Гл. ред. физ-мат. лит., 1985. — 432 с.
• Болтачев Г.Ш. Численные методы в теплофизике. Курс лекций Лекция 3: Численное интегрирование

См. также

• Численное дифференцирование
• Вычислительные методы
• Список квадратурных формул