Хабрахабр

Вычисление определённых интегралов: базовые алгоритмы

Обычно эти методы реализованы в стандартных математических библиотеках, таких как GNU Scientific Library для C, SciPy для Python и других. image
В этой публикации описаны простейшие методы вычисления интегралов функций от одной переменной на отрезке, также называемые квадратурными формулами. Также хотелось бы отметить связь квадратурных формул и методов численного интегрирования обыкновенных дифференциальных уравнений, о которых хочу написать ещё одну публикацию. Публикация имеет целью продемонстрировать, как эти методы работают "под капотом", и обратить внимание на некоторые вопросы точности и производительности алгоритмов.

Определение интеграла

Интегралом (по Риману) от функции $f(x)$ на отрезке $[a;b]$ называется следующий предел:

$\int_a^b f(x)dx = \lim_ \sum_{i=0}^{n-1} f(\xi_i)(x_{i+1} - x_i),~~(1)$

где $\Delta x = \max\lbrace x_{i+1} - x_i\rbrace$ — мелкость разбиения, $x_0 = a$, $x_n = b$, $\xi_i$ — произвольное число на отрезке $[x_i; x_{i+1}]$.

Если интеграл от функции существует, то значение предела одно и то же вне зависимости от разбиения, лишь бы оно было достаточно мелким.
image
Более наглядно геометрическое определение — интеграл равен площади криволинейной трапеции, ограниченной осью 0x, графиком функции и прямыми x = a и x = b (закрашенная область на рисунке).

Квадратурные формулы

Определение интеграла (1) можно переписать в виде

$I = \int_a^b f(x)dx \approx I_n = (b - a)\sum_{i=0}^{n-1} w_i f(\xi_i),~~(2)$

где $w_i$ — весовые коэффициенты, сумма которых должна быть равна 1, а сами коэффициенты — стремиться к нулю при увеличении числа $n$ точек, в которых вычисляется функция.

формул для приближенного вычисления интеграла). Выражение (2) — основа всех квадратурных формул (т.е. Задача состоит в том, чтобы выбрать точки $\lbrace \xi_i \rbrace$ и веса $w_i$ таким образом, чтобы сумма в правой части приближала требуемый интеграл как можно точнее.

Вычислительная задача

Задана функция $f(x)$, для которой есть алгоритм вычисления значений в любой точке отрезка $[a;b]$ (имеются в виду точки, представимые числом с плавающей точкой — никаких там функций Дирихле!).

6. Требуется найти приближённое значение интеграла $\int_a^b f(x)dx$.
Решения будут реализованы на языке Python 3.

Для проверки методов используется интеграл $\int_0^{3/2} \left[ 2x + \frac{1}{\sqrt{x + 1/16}}\right]dx = 17/4$.

Кусочно-постоянная аппроксимация

Идейно простейшие квадратурные формулы возникают из применения выражения (1) "в лоб":

$I_n = \sum_{i=0}^{n-1} f(\xi_i) (x_{i+1} - x_i)$

от метода разбиения отрезка точками $\lbrace x_i\rbrace$ и выбора точек $\lbrace\xi_i\rbrace$ значение предела не зависит, то выберем их так, чтобы они удобно вычислялись — например, разбиение возьмём равномерным, а для точек вычисления функции рассмотрим варианты: 1) $\xi_i = x_i$; 2) $\xi_i = x_{i+1}$; 3) $\xi_i = (x_i + x_{i+1}) / 2$. Т.к.

Получаем методы левых прямоугольников, правых прямоугольников и прямоугольников со средней точкой, соответственно.

Реализация

def _rectangle_rule(func, a, b, nseg, frac): """Обобщённое правило прямоугольников.""" dx = 1.0 * (b - a) / nseg sum = 0.0 xstart = a + frac * dx # 0 <= frac <= 1 задаёт долю смещения точки, # в которой вычисляется функция, # от левого края отрезка dx for i in range(npoints): sum += func(xstart + i * dx) return sum * dx def left_rectangle_rule(func, a, b, nseg): """Правило левых прямоугольников""" return _rectangle_rule(func, a, b, nseg, 0.0) def right_rectangle_rule(func, a, b, nseg): """Правило правых прямоугольников""" return _rectangle_rule(func, a, b, npoints, 1.0) def midpoint_rectangle_rule(func, a, b, nseg): """Правило прямоугольников со средней точкой""" return _rectangle_rule(func, a, b, nseg, 0.5)

Для анализа производительности квадратурных формул построим график погрешности в координатах "число точек — отличие численного результата от точного".

image

Что можно заметить:

  1. Формула со средней точкой гораздо точнее, чем с правой или левой точками
  2. Погрешность формулы со средней точкой падает быстрее, чем у двух остальных
  3. При очень мелком разбиении погрешность формулы со средней точкой начинает возрастать
    Первые два пункта связаны с тем, что формула прямоугольников со средней точкой имеет второй порядок аппроксимации, т.е. $|I_n - I| = O( 1/n^2)$, а формулы правых и левых прямоугольников — первый порядок, т.е. $|I_n - I| = O(1/n)$.
    Возрастание погрешности при измельчении шага интегрирования связано с нарастанием погрешности округления при суммировании большого числа слагаемых. Эта ошибка растёт как $|I_n - I| = O(1/n)$, что не даёт при интегрировании достигнуть машинной точности.
    Вывод: методы прямоугольников с правой и левой точками имеют низкую точность, которая к тому же медленно растёт с измельчением разбиения. Поэтому они имеют смысл разве что в демонстрационных целях. Метод прямоугольников со средней точкой имеет более высокий порядок аппроксимации, что даёт ему шансы на использование в реальных приложениях (об этом чуть ниже).

Кусочно-линейная аппроксимация

Следующий логический шаг — аппроксимировать интегрируемую функцию на каждом из подотрезков линейной функцией, что даёт квадратурную формулу трапеций:

$I_n = \sum_{i=0}^{n-1} \frac{f(x_i) + f(x_{i+1})}{2}(x_{i+1} - x_i)~~(3)$

image
Иллюстрация метода трапеций для n=1 и n=2.

В случае равномерной сетки длины всех отрезков разбиения равны, и формула имеет вид

$I_n = h\left(\frac{f(a) + f(b)}{2} + \sum_{i=1}^{n-1} f(a + ih)\right),~h=\frac{b-a}{n}~~(3a)$

Реализация

def trapezoid_rule(func, a, b, nseg): """Правило трапеций nseg - число отрезков, на которые разбивается [a;b]""" dx = 1.0 * (b - a) / nseg sum = 0.5 * (func(a) + func(b)) for i in range(1, nseg): sum += func(a + i * dx) return sum * dx

Построив график ошибки от числа точек разбиения, убеждаемся, что метод трапеций тоже имеет второй порядок аппроксимации и вообще даёт результаты, слабо отличающиеся от метода прямоугольников со средней точкой (в дальнейшем — просто метод прямоугольников).
image

Контроль точности вычисления

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

Один из простых методов оценки погрешности — правило Рунге — разность значений интегралов, рассчитанных по n и 2n точкам, даёт оценку погрешности: $\Delta_{2n} \approx |I_{2n} - I_n|$. Как это реализовать? При расчёте методом трапеций для удвоения числа точек нужны новые значения функции только в серединах отрезков предыдущего разбиения, т.е. Метод трапеций удобнее для удвоения мелкости разбиения, чем метод прямоугольников с центральной точкой. предыдущее приближение интеграла можно использовать для вычисления следующего.

Чем ещё хорош метод прямоугольников

Это означает, что его можно использовать для функций, имеющих на краях отрезка интегрируемые особенности (например, sinx/x или x-1/2 от 0 до 1). Метод прямоугольников не требует вычислять значения функции на концах отрезка. Отличие от метода трапеций лишь в том, что при уменьшении шага вдвое отбрасывается результат предыдущих вычислений, однако можно утроить число точек, и тогда предыдущее значение интеграла также можно использовать для вычисления нового. Поэтому показанный далее метод экстраполяции будет работать точно так же и для метода прямоугольников. Формулы для экстраполяции в этом случае необходимо скорректировать на другое соотношение шагов интегрирования.

Отсюда получаем следующий код для метода трапеций с контролем точности:

def trapezoid_rule(func, a, b, rtol = 1e-8, nseg0 = 1): """Правило трапеций rtol - желаемая относительная точность вычислений nseg0 - начальное число отрезков разбиения""" nseg = nseg0 old_ans = 0.0 dx = 1.0 * (b - a) / nseg ans = 0.5 * (func(a) + func(b)) for i in range(1, nseg): ans += func(a + i * dx) ans *= dx err_est = max(1, abs(ans)) while (err_est > abs(rtol * ans)): old_ans = ans ans = 0.5 * (ans + midpoint_rectangle_rule(func, a, b, nseg)) # новые точки для уточнения интеграла # добавляются ровно в середины предыдущих отрезков nseg *= 2 err_est = abs(ans - old_ans) return ans

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

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

Кусочно-параболическая аппроксимация

Для этого требуется, чтобы число отрезков разбиения было чётным, тогда параболы могут быть проведены через тройки точек с абсциссами {(x0=a, x1, x2), (x2, x3, x4), ..., (xn-2, xn-1, xn=b)}. Следующим шагом аппроксимируем функцию элементами парабол.


Иллюстрация кусочно-параболического приближения на 3 и 5 точках (n=2 и n=3).

Приближая интеграл от функции на каждом из отрезков [xk;xk+2] интегралом от параболической аппроксимации на этом отрезке и считая точки равномерно распределенными (xk+1=xk+h), получаем формулу Симпсона:

+4f(b-h) + f(b)]~~(4)$" data-tex="display"/> <img src="https://habrastorage.org/getpro/habr/formulas/f92/668/4c2/f926684c22abdead9a95c1d7e987618b.svg" alt="$I_{Simps, n} = \sum_{i=0}^{n/2-1}\frac{h}{3}[f(x_{2i})+4f(x_{2i+1})+f(x_{2i + 2})] = \\ =\frac{h}{3}[f(a)+4f(a+h)+2f(a+2h) + ...

Из формулы (4) напрямую получается "наивная" реализация метода Симпсона:

Заголовок спойлера

def simpson_rule(func, a, b, nseg): """Правило трапеций nseg - число отрезков, на которые разбивается [a;b]""" if nseg%2 = 1: nseg += 1 dx = 1.0 * (b - a) / nseg sum = (func(a) + 4 * func(a + dx) + func(b)) for i in range(1, nseg / 2): sum += 2 * func(a + (2 * i) * dx) + 4 * func(a + (2 * i + 1) * dx) return sum * dx / 3

Для оценки погрешности можно использовать точно так же вычисление интеграла с шагами h и h/2 — но вот незадача, при вычислении интеграла с более мелким шагом результат предыдущего вычисления придётся отбросить, хотя половина новых вычислений функции будет в тех же точках, что и раньше.

Присмотревшись повнимательнее, заметим, что интеграл по формуле Симпсона может быть представлен через два интеграла по формуле трапеций с разными шагами. Бесполезной траты машинного времени, к счастью, можно избежать, если реализовать метод Симпсона более хитроумным образом. Яснее всего это видно на базовом случае аппроксимации интеграла по трём точкам $(a, f_0),~(a+h, f_1),~(a+2h, f_2)$:

$I_{Simps,2} = \frac{h}{3}(f_0 + 4f_1 + f_2) = \frac{4}{3}h\left(\frac{f_0 + f_1}{2} + \frac{f_1 + f_2}{2}\right) - \frac{1}{3}\cdot2h\frac{f_0 + f_2}{2} = \\ =\frac{4I_{trap,2} - I_{trap,1}}{3}$

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

Как-то так...

class Quadrature: """Базовые определения для квадратурных формул""" __sum = 0.0 __nseg = 1 # число отрезков разбиения __ncalls = 0 # считает число вызовов интегрируемой функции def __restart(func, x0, x1, nseg0, reset_calls = True): """Обнуление всех счётчиков и аккумуляторов. Возвращает интеграл методом трапеций на начальном разбиении""" if reset_calls: Quadrature.__ncalls = 0 Quadrature.__nseg = nseg0 # вычисление суммы для метода трапеций с начальным числом отрезков разбиения nseg0 Quadrature.__sum = 0.5 * (func(x0) + func(x1)) dx = 1.0 * (x1 - x0) / nseg0 for i in range(1, nseg0): Quadrature.__sum += func(x0 + i * dx) Quadrature.__ncalls += 1 + nseg0 return Quadrature.__sum * dx def __double_nseg(func, x0, x1): """Вдвое измельчает разбиение. Возвращает интеграл методом трапеций на новом разбиении""" nseg = Quadrature.__nseg dx = (x1 - x0) / nseg x = x0 + 0.5 * dx i = 0 AddedSum = 0.0 for i in range(nseg): AddedSum += func(x + i * dx) Quadrature.__sum += AddedSum Quadrature.__nseg *= 2 Quadrature.__ncalls += nseg return Quadrature.__sum * 0.5 * dx def trapezoid(func, x0, x1, rtol = 1e-10, nseg0 = 1): """Интегрирование методом трапеций с заданной точностью. rtol - относительная точность, nseg0 - число отрезков начального разбиения""" ans = Quadrature.__restart(func, x0, x1, nseg0) old_ans = 0.0 err_est = max(1, abs(ans)) while (err_est > abs(rtol * ans)): old_ans = ans ans = Quadrature.__double_nseg(func, x0, x1) err_est = abs(old_ans - ans) print("Total function calls: " + str(Quadrature.__ncalls)) return ans def simpson(func, x0, x1, rtol = 1.0e-10, nseg0 = 1): """Интегрирование методом парабол с заданной точностью. rtol - относительная точность, nseg0 - число отрезков начального разбиения""" old_trapez_sum = Quadrature.__restart(func, x0, x1, nseg0) new_trapez_sum = Quadrature.__double_nseg(func, x0, x1) ans = (4 * new_trapez_sum - old_trapez_sum) / 3 old_ans = 0.0 err_est = max(1, abs(ans)) while (err_est > abs(rtol * ans)): old_ans = ans old_trapez_sum = new_trapez_sum new_trapez_sum = Quadrature.__double_nseg(func, x0, x1) ans = (4 * new_trapez_sum - old_trapez_sum) / 3 err_est = abs(old_ans - ans) print("Total function calls: " + str(Quadrature.__ncalls)) return ans

Сравним эффективность метода трапеций и парабол:

>>> import math
>>> Quadrature.trapezoid(lambda x: 2 * x + 1 / math.sqrt(x + 1 / 16), 0, 1.5, rtol=1e-9)
Total function calls: 65537
4.250000001385811
>>> Quadrature.simpson(lambda x: 2 * x + 1 / math.sqrt(x + 1 / 16), 0, 1.5, rtol=1e-9)
Total function calls: 2049
4.2500000000490985

Как видим, обоими методами ответ можно получть с достаточно высокой точностью, но количество вызовов подынтегральной функции разительно отличается — метод более высокого порядка эффективнее в 32 раза!

ошибка численного интегрирования $|I_{Simps, n} - I| = O(1/n^4)$ (а интегралы от кубических многочленов с помощью этой формулы вычисляются с точностью до ошибок округления при любом чётном n>0!).
image
Отсюда и возникает такой рост эффективности по сравнению с простой формулой трапеций. Построив график погрешности интегрирования от числа шагов, можно убедиться, что порядок аппроксимации формулы Симпсона равен четырём, т.е.

Что дальше?

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

Ошибка квадратурной формулы представляется в виде ряда по степеням шага интегрирования h. Попробуем пойти другим путём. Замечательное свойство метода трапеций (и прямоугольников со средней точкой!) в том, что для неё этот ряд состоит только из чётных степеней:

$I_{trap,n}[f, a, b] = \int_a^b f(x)dx + C_2h^2 + C_4h^4 + C_6h^6 + ...,~h = \frac{b-a}{n}~~~(5)$

На нахождении последовательных приближений к этому разложению основана экстраполяция Ричардсона: вместо того, чтобы приближать подынтегральную функцию многочленом, по рассчитанным приближениям интеграла $I(h)$ строится полиномиальная аппроксимация, которая при h=0 должна давать наилучшее приближение к истинному значению интеграла.

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

Поскольку приведённая реализация легко позволяет дробить разбиение вдвое, удобно рассматривать формулы для шагов h и h/2. Если считать, что каждое последующее слагаемое меньше предыдущего, то можно последовательно исключать степени h, имея приближения интеграла, рассчитанные с разными шагами.

$I_{trap,n} - I \approx C_2h^2;~I_{trap,2n} - I \approx C_2\left(\frac{h}{2}\right)^2$

Легко показать, что исключение старшего члена погрешности формулы трапеций в точности даст формулу Симпсона:

$I = I_{trap,2n} - C_2\left(\frac{h}{2}\right)^2 + O(h^4) \approx I_{trap,2n} - \frac{I_{trap,2n}-I_{trap,n}}{1-2^2} = I_{Simps,2n}$

Повторяя аналогичную процедуру для формулы Симпсона, получаем:

$I_{Simps,2n} - I \approx C_4\left(\frac{h}{2}\right)^4;~I_{Simps,n} - I \approx C_4h^4$

$I = I_{Simps,2n} - C_4\left(\frac{h}{2}\right)^4 + O(h^6) \approx I_{Simps,2n} - \frac{I_{Simps,2n}-I_{Simps,n}}{1-2^4}$

Если продолжить, вырисовывается такая таблица:

2 порядок

4 порядок

6 порядок

...

I0,0

I1,0

I1,1

I2,0

I2,1

I2,2

...

...

...

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

Элементы таблицы, как можно вывести из разложения (5), связаны рекуррентным соотношением:

$I_{i,j} = I_{i,j-1} - \frac{I_{i,j-1}-I_{i-1,j-1}}{1-\left(\frac{h_{i-j}}{h_i}\right)^2} = I_{i,j-1} - \frac{I_{i,j-1}-I_{i-1,j-1}}{1-2^{2j}}~~~(6)$

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

$\Delta_{i,j} \approx I_{i,j} - I_{i,j-1}$

Если метод Симпсона учитывает два предыдущих значения по методу трапеций, то метод Ромберга использует все ранее вычисленные методом трапеций значения для получения более точной оценки интеграла. Применение экстраполяции Ричардсона вместе с интегрированием методом трапеций называется методом Ромберга.

Реализация

Дополнительный метод добавляется в класс Quadrature

class Quadrature: """Базовые определения для квадратурных формул""" __sum = 0.0 __nseg = 1 # число отрезков разбиения __ncalls = 0 # считает число вызовов интегрируемой функции def __restart(func, x0, x1, nseg0, reset_calls = True): """Обнуление всех счётчиков и аккумуляторов. Возвращает интеграл методом трапеций на начальном разбиении""" if reset_calls: Quadrature.__ncalls = 0 Quadrature.__nseg = nseg0 # вычисление суммы для метода трапеций с начальным разбиением на nseg0 отрезков Quadrature.__sum = 0.5 * (func(x0) + func(x1)) dx = 1.0 * (x1 - x0) / nseg0 for i in range(1, nseg0): Quadrature.__sum += func(x0 + i * dx) Quadrature.__ncalls += 1 + nseg0 return Quadrature.__sum * dx def __double_nseg(func, x0, x1): """Вдвое измельчает разбиение. Возвращает интеграл методом трапеций на новом разбиении""" nseg = Quadrature.__nseg dx = (x1 - x0) / nseg x = x0 + 0.5 * dx i = 0 AddedSum = 0.0 for i in range(nseg): AddedSum += func(x + i * dx) Quadrature.__sum += AddedSum Quadrature.__nseg *= 2 Quadrature.__ncalls += nseg return Quadrature.__sum * 0.5 * dx def romberg(func, x0, x1, rtol = 1e-10, nseg0 = 1, maxcol = 5, reset_calls = True): """Интегрирование методом Ромберга nseg0 - начальное число отрезков разбиения maxcol - максимальный столбец таблицы""" # инициализация таблицы Itable = [[Quadrature.__restart(func, x0, x1, nseg0, reset_calls)]] i = 0 maxcol = max(0, maxcol) ans = Itable[i][i] error_est = max(1, abs(ans)) while (error_est > abs(rtol * ans)): old_ans = ans i += 1 d = 4.0 ans_col = min(i, maxcol) Itable.append([Quadrature.__double_nseg(func, x0, x1)] * (ans_col + 1)) for j in range(0, ans_col): diff = Itable[i][j] - Itable[i - 1][j] Itable[i][j + 1] = Itable[i][j] + diff / (d - 1.0) d *= 4.0 ans = Itable[i][ans_col] if (maxcol <= 1): # методы трапеций и парабол обрабатываются отдельно error_est = abs(ans - Itable[i - 1][-1]) elif (i > maxcol): error_est = abs(ans - Itable[i][min(i - maxcol - 1, maxcol - 1)]) else: error_est = abs(ans - Itable[i - 1][i - 1]) print("Total function calls: " + str(Quadrature.__ncalls)) return ans

Проверим, как работает аппроксимация высокого порядка:

>>> Quadrature.romberg(lambda x: 2 * x + 1 / math.sqrt(x + 1/16), 0, 1.5, rtol=1e-9, maxcol = 0) # трапеции
Total function calls: 65537
4.250000001385811
>>> Quadrature.romberg(lambda x: 2 * x + 1 / math.sqrt(x + 1/16), 0, 1.5, rtol=1e-9, maxcol = 1) # параболы
Total function calls: 2049
4.2500000000490985
>>> Quadrature.romberg(lambda x: 2 * x + 1 / math.sqrt(x + 1/16), 0, 1.5, rtol=1e-9, maxcol = 4)
Total function calls: 257
4.250000001644076

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

Некоторые замечания

Уменьшение числа вычислений подынтегрального выражения не только экономит вычислительные ресурсы (хотя при более оптимизированной реализации и это тоже), но и уменьшает влияние погрешностей округления на результат. Замечание 1. Количество вызовов функции в этих задачах характеризует число суммирований при вычислении интеграла. То есть, повышается не только точность интегрирования при заданном числе вызовов функции, но и предельно достижимая точность вычисления интеграла. Так, при попытке вычислить интеграл тестовой функции метод трапеций зависает при попытке достигнуть относительной точности 5×10-15, метод парабол — при желаемой точности 2×10-16(что является пределом для чисел в двойной точности), а метод Ромберга справляется с вычислением тестового интеграла вплоть до машинной точности (с ошибкой в младшем бите).

В первую очередь, это относится к случаям, когда задаваемая погрешность близка к машинной точности. Замечание 2. Если метод сходится при задании некоторой точности, это не означает, что вычисленное значение интеграла имеет ту же самую точность.

Это значит, что для функций с изломами или разрывами он может оказаться хуже простых методов. Замечание 3. Хотя метод Ромберга для ряда функций работает почти магическим образом, он предполагает наличие у подынтегральной функции ограниченных производных высоких порядков. Например, проинтегрируем f(x)=|x|:

>>> Quadrature.trapezoid(abs, -1, 3, rtol=1e-5)
Total function calls: 9
5.0
>>> Quadrature.simpson(abs, -1, 3, rtol=1e-5)
Total function calls: 17
5.0
>>> Quadrature.romberg(abs, -1, 3, rtol=1e-5, maxcol = 2)
Total function calls: 17
5.0
>>> Quadrature.romberg(abs, -1, 3, rtol=1e-5, maxcol = 3)
Total function calls: 33
5.0
>>> Quadrature.romberg(abs, -1, 3, rtol=1e-5, maxcol = 4)
Total function calls: 33
5.000001383269357

На самом деле, лучше ограничить число столбцов таблицы Ромберга на уровне 4-6. Замечание 4. Может показаться, что чем выше порядок аппроксимации, тем лучше. Второе слагаемое представляет собой разность двух последовательных элементов j-1-го столбца, поделенную на примерно 4j. Чтобы понять это, посмотрим на формулу (6). в j-1-м столбце находятся аппроксимации интеграла порядка 2j, то сама разность имеет порядок (1/ni)2j ~ 4-ij. Т.к. Т.е. C учётом деления получается ~4-(i+1)j ~ 4-j2. при j~7 второе слагаемое в (6) теряет точность после приведения порядков при сложении чисел с плавающей точкой, и повышение порядка аппроксимации может вести к накоплению ошибки округления.

Как говорится, почувствуйте разницу. Замечание 5. Желающие могут ради интереса применить описанные методы для нахождения интеграла $\int_0^1\sqrt{x}\sin{x}dx$ и эквивалентного ему $\int_0^1 2t^2\sin{t^2}dt$.

Заключение

Продемонстрировано, как с помощью несложной модификации получить на базе метода трапеций класс квадратурных формул по методу Ромберга, что значительно ускоряет сходимость численного интегрирования. Представлено описание и реализация базовых методов численного интегрирования функций на равномерной сетке. слабо меняющихся на отрезке интегрирования, не имеющих особенностей на краях отрезка (см. Метод хорошо работает для интегрирования "обычных" функций, т.е. Замечание 5), быстрых осцилляций и т.д.

Продвинутые методы численного интегрирования для более сложных случаев можно найти в книгах из списка литературы (в [3] — с примерами реализации на C++).

Литература

  1. А.А. Самарский, А.В. Гулин. Численные методы. М.: Наука. 1989.
  2. J. Stoer, R. Bulirsch. Introduction to Numerical Analysis: Second Edition. Springer-Verlag New York. 1993.
  3. W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery. Numerical Recipes: Third Edition. Cambridge University Press. 2007.
Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть