Численное интегрирование
Численное интегрирование (историческое название: (численная) квадратура) — вычисление значения определённого интеграла (как правило, приближённое). Под численным интегрированием понимают набор численных методов для нахождения значения определённого интеграла.
Численное интегрирование применяется, когда:
- Сама подынтегральная функция не задана аналитически. Например, она представлена в виде таблицы (массива) значений в узлах некоторой расчётной сетки.
- Аналитическое представление подынтегральной функции известно, но её первообразная не выражается через аналитические функции. Например, <math>f(x) = \exp(-x^2)</math>.
В этих двух случаях невозможно вычисление интеграла по формуле Ньютона — Лейбница. Также возможна ситуация, когда вид первообразной настолько сложен, что быстрее вычислить значение интеграла численным методом.
Одномерный случай
Основная идея большинства методов численного интегрирования состоит в замене подынтегральной функции на более простую, интеграл от которой легко вычисляется аналитически. При этом для оценки значения интеграла получаются формулы вида
- <math>I \approx \sum_{i=1}^{n} w_i\, f(x_i),</math>
где <math>n</math> — число точек, в которых вычисляется значение подынтегральной функции. Точки <math>x_i</math> называются узлами метода, числа <math>w_i</math> — весами узлов. При замене подынтегральной функции на полином нулевой, первой и второй степени получаются соответственно методы прямоугольников, трапеций и парабол (Симпсона). Часто формулы для оценки значения интеграла называют квадратурными формулами.
Частным случаем является метод построения интегральных квадратурных формул для равномерных сеток, известный как формулы Котеса. Метод назван в честь Роджера Котса. Основной идеей метода является замена подынтегральной функции каким-либо интерполяционным многочленом. После взятия интеграла можно написать
- <math>\int\limits_a^b f(x)\,dx = \sum_{i=0}^{n} H_i\, f(x_i) + r_n(f),</math>
где числа <math>H_i</math> называются коэффициентами Котеса и вычисляются как интегралы от соответствующих многочленов, стоящих в исходном интерполяционном многочлене для подынтегральной функции при значении функции в узле <math>x_i = a + i h</math> (<math>h = (b-a)/n</math> — шаг сетки; <math>n</math> — число узлов сетки, а индекс узлов <math>i=0 \ldots n</math>). Слагаемое <math>r_n(f)</math> — погрешность метода, которая может быть найдена разными способами. Для нечетных <math>n \geqslant 1</math> погрешность может быть найдена интегрированием погрешности интерполяционного полинома подынтегральной функции.
Частными случаями формул Котеса являются: формулы прямоугольников (<math>n=0</math>), формулы трапеций (<math>n=1</math>), формула Симпсона (<math>n=2</math>), формула Ньютона (<math>n=3</math>) и т. д.
Метод прямоугольников
Шаблон:Main Пусть требуется определить значение интеграла функции на отрезке <math>\left[ {a},{b} \right]</math>. Этот отрезок делится точками <math>x_0, x_1, \ldots, x_{n-1}, x_n</math> на <math>n</math> равных отрезков длиной <math>\Delta {x} = \frac{b-a}{n}.</math> Обозначим через <math>y_0, y_1, \ldots, y_{n-1}, y_n</math> значение функции <math>f\left(x\right)</math> в точках <math>x_0, x_1, \ldots, x_{n-1}, x_n.</math> Далее составляем суммы <math>y_0 \,\Delta {x} + y_1 \,\Delta {x} + \ldots + y_{n-1} \,\Delta {x}.</math> Каждая из сумм — интегральная сумма для <math>f\left(x\right)</math> на <math>\left[ {a},{b} \right]</math> и поэтому приближённо выражает интеграл
- <math>\int\limits_a^b f(x)\,dx \approx \frac{b-a}{n} (y_0 + y_1 + \ldots + y_{n-1}).</math>
Если заданная функция — положительная и возрастающая, то эта формула выражает площадь ступенчатой фигуры, составленной из «входящих» прямоугольников, также называемая формулой левых прямоугольников, а формула
- <math>\int\limits_a^b f(x)\,dx \approx \frac{b-a}{n} (y_1 + y_2 + \ldots + y_n)</math>
выражает площадь ступенчатой фигуры, состоящей из «выходящих» прямоугольников, также называемая формулой правых прямоугольников. Чем меньше длина отрезков, на которые делится отрезок <math>\left[ {a},{b} \right]</math>, тем точнее значение, вычисляемое по этой формуле, искомого интеграла.
Очевидно, стоит рассчитывать на бо́льшую точность, если брать в качестве опорной точки для нахождения высоты точку посередине промежутка. В результате получаем формулу средних прямоугольников:
- <math>\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),</math>
где <math>h = \frac{b-a}{n}</math>
Учитывая априорно бо́льшую точность последней формулы при том же объёме и характере вычислений её называют формулой прямоугольников
Метод трапеций
Шаблон:Main Если функцию на каждом из частичных отрезков аппроксимировать прямой, проходящей через конечные значения, то получим метод трапеций.
Площадь трапеции на каждом отрезке:
- <math>I_i \approx \frac{f(x_{i-1})+f(x_{i})}{2} (x_{i}-x_{i-1}).</math>
Погрешность аппроксимации на каждом отрезке:
- <math>\left| R_{i} \right| \leqslant \frac{\left( b-a \right)^3}{12n^2} M_{2,i},</math> где <math>M_{2,i}=\max_{x\in[x_{i-1},\,x_i]}{\left| f(x) \right|}.</math>
Полная формула трапеций в случае деления всего промежутка интегрирования на отрезки одинаковой длины <math>h</math>:
- <math>I \approx h\left( \frac{f(x_{0})+f(x_{n})}{2} + \sum_{i=1}^{n-1}f(x_{i})\right),</math> где <math>h=\frac{b-a}{n}.</math>
Погрешность формулы трапеций:
- <math>\left| R \right| \leqslant \frac{\left( b-a \right)^3}{12n^2} M_{2},</math> где <math>M_{2}=\max_{x\in[a,\,b]}{\left| f(x) \right|}.</math>
Метод парабол (метод Симпсона)
Шаблон:Main Использовав три точки отрезка интегрирования, можно заменить подынтегральную функцию параболой. Обычно в качестве таких точек используют концы отрезка и его среднюю точку. В этом случае формула имеет очень простой вид
- <math>I \approx \frac{b-a}{6}\left(f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right)</math>.
Если разбить интервал интегрирования на <math>2N</math> равных частей, то имеем
- <math>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),</math>
где <math>f_i = f\left(a+\frac{(b-a)i}{2N}\right)</math>.
Увеличение точности
Приближение функции одним полиномом на всем отрезке интегрирования, как правило, приводит к большой ошибке в оценке значения интеграла.
Для уменьшения погрешности отрезок интегрирования разбивают на части и применяют численный метод для оценки интеграла на каждой из них.
При стремлении количества разбиений к бесконечности оценка интеграла стремится к его истинному значению для аналитических функций для любого численного метода.
Приведённые выше методы допускают простую процедуру уменьшения шага в два раза, при этом на каждом шаге требуется вычислять значения функции только во вновь добавленных узлах. Для оценки погрешности вычислений используется правило Рунге.
Метод Гаусса
Шаблон:Main Описанные выше методы используют фиксированные точки отрезка (концы и середину) и имеют низкий порядок точности (0 — методы правых и левых прямоугольников, 1 — методы средних прямоугольников и трапеций, 3 — метод парабол (Симпсона)). Если мы можем выбирать точки, в которых мы вычисляем значения функции <math>f(x)</math>, то можно при том же количестве вычислений подынтегральной функции получить методы более высокого порядка точности. Так, для двух (как в методе трапеций) вычислений значений подынтегральной функции можно получить метод уже не второго, а третьего порядка точности:
- <math>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)</math>.
В общем случае, используя <math>n</math> точек, по формуле <math>I \approx \sum_{i=1}^{n} a_i\, f(x_i)</math> можно получить метод с порядком точности <math>2n-1</math>, т. е. получаются точные значения для полиномов степени не выше <math>2n-1</math>.
Значения узлов <math>x_i</math> метода Гаусса по <math>n</math> точкам являются корнями полинома Лежандра степени <math>n</math>. Значения весов вычисляются по формуле <math>a_i = \frac{2}{(1-x_i^2)\,[P_n'(x_i)]^2}</math>, где <math>P_n'</math> - первая производная полинома Лежандра.
Для <math>n=3</math> узлы и веса имеют следующие значения : <math>x_{1,3}=\pm\sqrt{0.6}, x_2=0,</math> веса : <math> a_{1,3}=\frac{5}{9}, a_2=\frac{8}{9}</math>
(полином определен на отрезке <math>[-1,1]</math>).
Наиболее известен метод Гаусса по пяти точкам.
Метод Гаусса — Кронрода
Недостаток метода Гаусса состоит в том, что он не имеет лёгкого (с вычислительной точки зрения) пути оценки погрешности полученного значения интеграла. Использование правила Рунге требует вычисления подынтегральной функции примерно в таком же числе точек, не давая при этом практически никакого выигрыша точности, в отличие от простых методов, где точность увеличивается в несколько раз при каждом новом разбиении. Кронродом был предложен следующий метод оценки значения интеграла
- <math>I \approx \sum_{i=1}^{n} a_i\, f(x_i) + \sum_{i=1}^{n+1} b_i\, f(y_i)</math>,
где <math>x_i</math> — узлы метода Гаусса по <math>n</math> точкам, а <math>3n+2</math> параметров <math>a_i</math>, <math>b_i</math>, <math>y_i</math> подобраны таким образом, чтобы порядок точности метода был равен <math>3n+1</math>.
Тогда для оценки погрешности можно использовать эмпирическую формулу:
- <math>\Delta = \left(200 |I - I_G|\right)^{1.5}</math>,
где <math>I_G</math> — приближённое значение интеграла, полученное методом Гаусса по <math>n</math> точкам. Библиотеки gsl и SLATEC для вычисления определённых интегралов содержат подпрограммы, использующие метод Гаусса — Кронрода по 15, 21, 31, 41, 51 и 61 точкам. Библиотека ALGLIB использует метод Гаусса — Кронрода по 15 точкам.
Метод Чебышёва
Метод Чебышева (или как его иногда называют Гаусса — Чебышева) является одним из представителей методов наивысшей алгебраической точности Гаусса. Его отличительной особенностью является наличие у подынтегральной функции множителя <math>\left( \frac{1}{\sqrt{1-x^2}} \right)</math>. Суть заключается в следующем:
- <math>\int\limits_{-1}^{1} f(x)\left( \frac{1}{\sqrt{1-x^2}} \right)\, dx = \sum_{i=1}^N C_i f(x_i)</math>,
где <math> C_i = \pi/N </math>, <math>x_i=\cos \left( \frac{(2i-1)\pi}{2N} \right)</math>, <math>N</math> — количество узлов метода.
Метод Гаусса-Лагера
Метод Гаусса-Эрмита
Интегрирование при бесконечных пределах
Для интегрирования по бесконечным пределам нужно ввести неравномерную сетку, шаги которой нарастают при стремлении к бесконечности, либо можно сделать такую замену переменных в интеграле, после которой пределы будут конечны. Аналогичным образом можно поступить, если функция особая на концах отрезка интегрирования.
См. в том числе Метод Самокиша.
Методы Монте-Карло
Для определения площади под графиком функции можно использовать следующий стохастический алгоритм:
- ограничим функцию прямоугольником (n-мерным параллелепипедом в случае многих измерений), площадь которого <math>S_{par}</math> можно легко вычислить;
- «набросаем» в этот прямоугольник (параллелепипед) некоторое количество точек (<math>N</math> штук), координаты которых будем выбирать случайным образом;
- определим число точек (<math>K</math> штук), которые попадут под график функции;
- площадь области, ограниченной функцией и осями координат, <math>S</math> даётся выражением <math>S = S_{par}\frac{K}{N}</math>;
Этот алгоритм требует определения экстремумов функции на интервале и не использует вычисленное точное значение функции <math>f(x)</math> кроме как в сравнении, вследствие чего непригоден для практики. Приведённые в основной статье варианты метода Монте-Карло избавлены от этих недостатков.
Для малого числа измерений интегрируемой функции производительность Монте-Карло интегрирования гораздо ниже, чем производительность детерминированных методов. Тем не менее, в некоторых случаях, когда функция задана неявно, а необходимо определить область, заданную в виде сложных неравенств, стохастический метод может оказаться более предпочтительным.
Методы Рунге — Кутты
Шаблон:Main Ме́тоды Ру́нге — Ку́тты — важное семейство численных алгоритмов решения обыкновенных дифференциальных уравнений и их систем — итеративные методы явного и неявного приближённого вычисления, разработанные около 1900 года немецкими математиками К. Рунге и М. В. Куттой. Шаблон:Дополнить раздел
Метод сплайнов
Многомерный случай
В небольших размерностях можно так же применять квадратурные формулы, основанные на интерполяционных многочленах. Интегрирование производится аналогично одномерному интегрированию. Для больших размерностей эти методы становятся неприемлемыми из-за быстрого возрастания числа точек сетки и/или сложной границы области. В этом случае применяется метод Монте-Карло. Генерируются случайные точки в нашей области и усредняются значения функции в них. Так же можно использовать смешанный подход — разбить область на несколько частей, в каждой из которых (или только в тех, где интеграл посчитать не удаётся из-за сложной границы) применить метод Монте-Карло.
Примеры реализации
Ниже приведена реализация на Python 3 метода средних прямоугольников, метода средних трапеций, метода Симпсона и метода Монте-Карло. <source lang="python3"> 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()) </source>
Литература
- Шаблон:Книга
- Шаблон:Книга
- Шаблон:Книга
- Болтачев Г.Ш. Численные методы в теплофизике. Курс лекций Лекция 3: Численное интегрирование