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

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

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

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

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

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

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

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

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

I \approx \sum_{i=1}^{n} w_i\, f(x_i),

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

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

\int\limits_a^b f(x)\,dx = \sum_{i=0}^{n} H_i\, f(x_i) + r_n(f),

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

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

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

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

\int\limits_a^b f(x)\,dx \approx \frac{b-a}{n} (y_0 + y_1 + \ldots + y_{n-1}).

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

\int\limits_a^b f(x)\,dx \approx \frac{b-a}{n} (y_1 + y_2 + \ldots + y_n)

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

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

\int\limits_a^b f(x)\,dx \approx h \sum_{i=1}^{n}f(x_{i-1} + \frac{h}{2}) = h \sum_{i=1}^{n}f(x_i - \frac{h}{2}),

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

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

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

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

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

~I_i \approx \frac{f(x_{i-1})+f(x_{i})}{2} (x_{i}-x_{i-1})

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

~\left| R_{i} \right| \leqslant \frac{\left( b-a \right)^3}{12n^2} M_{2,i}\,, где M_{2,i}=\max_{x\mathcal{2}[x_{i-1},x_i]} \left| f''(x) \right|

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

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

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

~\left| R \right| \leqslant \frac{\left( b-a \right)^3}{12n^2} M_{2} \,, где M_{2}=\max_{x\mathcal{2}[a,b]}^{\color{White}[} \left| f''(x) \right|.

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

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

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

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

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 + ... +f_{2N-2}\right) + f_{2N}\right),

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

Увеличение точности[править | править исходный текст]

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

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

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

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

Метод Гаусса[править | править исходный текст]

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

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).

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

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

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

(Полином определен на отрезке [-1,1]).

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

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

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

I \approx \sum_{i=1}^{n} a_i\, f(x_i) + \sum_{i=1}^{n+1} b_i\, f(y_i),

где x_i — узлы метода Гаусса по n точкам, а 3n+2 параметров a_i, b_i, y_i подобраны таким образом, чтобы порядок точности метода был равен 3n+1.

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

\Delta = \left(200 |I - I_G|\right)^{1.5},

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

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

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

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

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

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

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

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

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

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

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

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

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

Метод Симпсона, Средних прямоугольников, Монте-карло и средних трапеций 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()

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

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

Литература[править | править исходный текст]

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

См. также[править | править исходный текст]