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

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

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

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

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

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

Одномерный случай[править | править вики-текст]

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

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

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

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

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

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

Метод прямоугольников[править | править вики-текст]

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

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

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

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

где

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

Метод трапеций[править | править вики-текст]

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

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

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

где

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

где

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

где

Метод парабол (метод Симпсона)[править | править вики-текст]

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

.

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

где .

Увеличение точности[править | править вики-текст]

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

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

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

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

Метод Гаусса[править | править вики-текст]

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

.

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

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

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

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

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

Метод Гаусса — Кронрода[править | править вики-текст]

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

,

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

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

,

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

Метод Чебышёва[править | править вики-текст]

Интегрирование при бесконечных пределах[править | править вики-текст]

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

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

Методы Монте-Карло[править | править вики-текст]

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

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

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

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

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

Методы Рунге — Кутты[править | править вики-текст]

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

Метод сплайнов[править | править вики-текст]

Многомерный случай[править | править вики-текст]

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

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

Примеры реализации[править | править вики-текст]

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

from matplotlib import mlab
import math, random

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

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

    d, n = 1, 1
    while math.fabs(d) > delta:
        d = (integrate(func, mim_lim, max_lim, n * 2) - integrate(func, mim_lim, max_lim, n)) / 3
        n *= 2

    print 'Rect'
    print ' '.join([
        '\t',
        str(n),
        str(math.fabs(integrate(func, mim_lim, max_lim, n))),
        str(math.fabs(integrate(func, mim_lim, max_lim, n)) + d)])

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

    d, n = 1, 1
    while math.fabs(d) > delta:
        d = (integrate(func, mim_lim, max_lim, n * 2) - integrate(func, mim_lim, max_lim, n)) / 3
        n *= 2

    print 'Rect'
    print ' '.join([
        '\t',
        str(n),
        str(math.fabs(integrate(func, mim_lim, max_lim, n))),
        str(math.fabs(integrate(func, mim_lim, max_lim, n)) + d)])

def simpson_method(func, mim_lim, max_lim, delta):
    def integrate(func, mim_lim, max_lim, n):
        integral = 0.0
        step = (max_lim - mim_lim) / n
        for x in mlab.frange(mim_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 math.fabs(d) > delta:
        d = (integrate(func, mim_lim, max_lim, n * 2) - integrate(func, mim_lim, max_lim, n)) / 15
        n *= 2

    print 'Rect'
    print ' '.join([
        '\t',
        str(n),
        str(math.fabs(integrate(func, mim_lim, max_lim, n))),
        str(math.fabs(integrate(func, mim_lim, max_lim, n)) +d )])

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

    print 'Mk'
    print '\t' + str(n) + ' ' + str(math.fabs(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 get_i()


Литература[править | править вики-текст]

См. также[править | править вики-текст]