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

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

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}$ точкам. Библиотеки и для вычисления определённых интегралов содержат подпрограммы, использующие метод Гаусса — Кронрода по 15, 21, 31, 41, 51 и 61 точкам. Библиотека использует метод Гаусса — Кронрода по 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}$ — количество узлов метода.

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

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

Этот раздел не завершён . Вы поможете проекту, исправив и дополнив его.

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

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

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

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

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

• ограничим функцию прямоугольником ( 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 с.

См. также

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

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