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

Материал из Википедии — свободной энциклопедии
Перейти к навигации Перейти к поиску

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

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

  1. Сама подынтегральная функция не задана аналитически. Например, она представлена в виде таблицы (массива) значений в узлах некоторой расчётной сетки.
  2. Аналитическое представление подынтегральной функции известно, но её первообразная не выражается через аналитические функции. Например, .

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

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

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

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

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

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

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

Частными случаями формул Котеса являются: формулы прямоугольников (), формулы трапеций (), формула Симпсона (), формула Ньютона () и т. д.

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

[править | править код]

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

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

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

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

где

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

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

[править | править код]

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

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

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

где

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

где

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

где

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

[править | править код]

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

.

Если разбить интервал интегрирования на равных частей, то имеем

где .

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

[править | править код]

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

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

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

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

Метод Гаусса

[править | править код]

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

.

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

Значения узлов метода Гаусса по точкам являются корнями полинома Лежандра степени . Значения весов вычисляются по формуле , где - первая производная полинома Лежандра.

Для узлы и веса имеют следующие значения : веса :

(полином определен на отрезке ).

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

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

[править | править код]

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

,

где  — узлы метода Гаусса по точкам, а параметров , , подобраны таким образом, чтобы порядок точности метода был равен .

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

,

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

Метод Чебышева (или как его иногда называют Гаусса — Чебышева) является одним из представителей методов наивысшей алгебраической точности Гаусса. Его отличительной особенностью является наличие у подынтегральной функции множителя . Суть заключается в следующем:

,

где , , — количество узлов метода.

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

[править | править код]

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

[править | править код]

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

[править | править код]

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

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

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

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

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

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

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

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

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

[править | править код]

Ме́тоды Ру́нге — Ку́тты — важное семейство численных алгоритмов решения обыкновенных дифференциальных уравнений и их систем — итеративные методы явного и неявного приближённого вычисления, разработанные около 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: Численное интегрирование