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

Базовая идея численного интегрирования состоит в том, чтобы аппроксимировать интеграл следующим образом: где выражение справа называется квадратурой. Как и в случае численного дифференцирования, простейшие квадратуры можно вывести, аппроксимируя $f(x)$ интерполяционным многочленом Лагранжа.

Формулы Ньютона-Котеса

Рассмотрим функцию $f(x) \in C^{n}[a; b]$ и $n$ различных узлов $x_1, \dots, x_n$. Тогда, раскладывая $f(x)$ в базисные многочлены Лагранжа, мы получаем где $\xi(x) \in (a; b)$. Тогда интегрирование $f(x)$ на интервале $[a; b]$ дает: что дает выражение для коэффициентов квадратуры $a_i$: В случае, когда интерполяционные узлы $x_1, \dots, x_n$ распределены равномерно, мы получаем формулы Ньютона-Котеса для численного интегрирования. Оценим остаточный член формул Ньютона-Котеса для $n$ узлов $x_i = x_1 + (i-1)h$, где $i = 1, \dots, n$ и $x_1 = a, x_n = b$: Представим переменную $x$ как $x = x_1 + (t-1) h$, что дает $dx = h dt$ и $t = \dfrac{x-x_1}{h} + 1$. Тогда, произведя замену в выражении выше, имеем: Таким образом, порядок точности можно оценить как $O(h^n)$. Чтобы получить более точное выражение для остаточного члена, необходимо заметить, что полином в подинтегральном выражении меняет знак только в узлах $x_i$. Тогда, пользуя первой теоремой о среднем значении, мы получаем: где $\xi_j \in [x_j; x_{j+1}]$. Эта форма $R$ приводит к теореме об остаточном члене формул Ньютона-Котеса, которую мы рассмотрим без доказательства. Основным следствием теоремы является тот факт, что предпочтительным является нечетное число узлов. Это связано с формой полинома в подинтегральном выражении в \eqref{eq:R_newton_cotes}, который в случае нечетного числа узлов имеет равное количество отрезков, на которых функция всюду положительна и всюду отрицальна, что снижает погрешность метода.

Формулы трапеций и Симпсона

Частными случаями формул Ньютона-Котеса являются формула трапеций ($n=2$) и формула Симпсона ($n=3$). В связи с их распространенностью рассмотрим вывод этих формул с явными выражением для остаточного члена. Так, для $n=2$ выражение \eqref{eq:general_numer_quad_lagrange} принимает вид: что после интегрирования линейных полиномов дает: Так как полином внутри подинтегрального выражения не меняет знак для $x \in [a; b]$, первая теорема о среднем значении позволяет записать выражение где $\xi \in (a; b)$. Таким образом мы получаем выражение для формулы трапеций: где, как мы помним, $x_1 = a$, $x_2 = b$ и $h = b - a$. Для вывода формулы Симпсона с $n=3$ мы воспользуемся альтернативным подходом, который сделает вывод остаточного члена проще, а также в очередной раз продемонстрирует связь между интерполяционным многочленом Лагранжа и рядом Тейлора. Рассмотрим разложение функции $f(x)$ в ряд Тейлора в промежуточном узле $x_2$: Тогда интеграл от $f(x)$ принимает вид: В силу четности степеней члены, включающие в себя нечетные производные, обращаются в ноль. Более того, так как $(x-x_2)^4 \geq 0$, мы можем применить первую теорему о среднем значении, что в результате дает: где $\xi \in (x_1; x_3)$. Для того, чтобы избавиться от второй производной, воспользуемся формулой для численного дифференцирования с остаточным членом \eqref{eq:num_diff_second_deriv}. После подстановки равенство выше принимает вид: где $\xi_1, \xi_2 \in (x_1; x_3)$. Так как в общем случае $\xi_1 \neq \xi_2$ нам необходимо каким-то образом скомбинировать два последних члена. Для этого рассмотрим остаточный член в общем случае и предположим, что существует такое $\xi \in (x_1; x_3)$, что: где $\alpha \in \mathbb{R}$ – неопределенный коэффициент. Для того чтобы найти его, заметим, что $f^{(4)}(\xi) = 24$ для любого нормированного многочлена 4-й степени. Тогда рассмотрим в качестве $f(x)$ многочлен $(x-x_2)^4$: Таким образом мы получаем формулу Симпсона с явным выражением для остаточного члена: где $x_1 = a$, $x_2 = x_1 + h$ и $x_3 = x_1 + 2h = b$. Малость погрешности формулы Симпсона, $O(h^5)$, при использовании всего лишь трех узлов объясняет ее частое использование в реальных приложениях. Необходимо отметить, что такая малая погрешность сохраняется только при достаточной гладкости функции $f(x)$. Для того, чтобы численно исследовать зависимость остаточного члена формулы Симпсона от гладкости интегрируемой функции, мы проведет следующий вычислительный эксперимент. Рассмотрим две функции $f_1(x) = e^x$ и $f_2(x) = |x|$, первая из которых является бесконечно гладкой, а вторая лишь непрерывной. Построим на основе последовательности шагов $\{h_i\}_{i=1}^n$ соответствующую последовательность интегралов: где $x = 1/2$ и $x = 0$ являются центральными узлами численного интегрирования для $f_1(x)$ и $f_2(x)$ соответственно. Рассчитав те же интегралы с помощью формулы Симпсона и найдя полную погрешность $E$, включающую в себя как остаточный член, так и вычислительную погрешность, для каждого случая, мы получаем рисунок ниже. Для начала заметим, что численное интегрирование устойчиво с вычислительной точки зрения, и полная погрешность полностью соответствует остаточному члену вплоть до тех пор, пока $E$ не достигнет машинного эпсилон. Сам остаточный член, как мы видим, пропорционален $O(h^5)$ в случае бесконечно гладкой функции $f_1(x)$, что и предполагается равенством \eqref{eq:simpson}, в то время как в случае функции $f_2(x)$ погрешность метода пропорциональна лишь $O(h^2)$, что связано с тем, что уже первая производная от $f_2(x)$ имеет разрыв, что делает формулу \eqref{eq:simpson} нерабочей.

Формула средних

Отдельно рассмотрим случай, когда на отрезке $[a; b]$ мы имеем только один узел, расположенный в центре отрезка, т.е. $x_1 = \frac{a + b}{2}$. Разложим функцию $f(x) \in C^2[a; b]$ в ряд Тейлора в этом узле: Тогда интеграл от $f(x)$ вычисляется следующим образом: Так как выражение $(x-x_1)^2$ всюду неотрицательно, мы можем воспользователься первой теоремой о среднем значении, что дает: что после интегрирования приводит к формуле средних: Заметим, что остаточный член в формуле средних пропорционален $O(h^3)$, где $h = b-a$, т.е. имеет тот же порядок погрешности, что и формула трапеций. Однако формула средних использует только один узел и, соответственно, только одно значение функции $f(x)$, что делает ее предпочтительной с вычислительной точки зрения.

Степень точности численного интегрирования

Несмотря на то, что оценка остаточного члена с точки зрения его зависимости от шага интегрирования $h$ является удобной для анализа и сравнения различных формул интегрирования, на практике часто пользуются другой оценкой, напрямую следующей из формы остаточного члена. Этой оценкой является степень точности численного интегрирования. Степенью точности квадратуры называется такое наибольшее $n \in \mathbb{N}$, что формула квадратуры дает точный результат для всех $x^i, i = 0, \dots n$. Легко проверить, что вследствие линейности операции интегрирования, это определение автоматически приводит к следующей теореме: Форма остаточного члена в формулах трапеций и Симпсона позволяет заключить, что они имеют степени точности $1$ и $3$ соответственно. Например, интегрируя любой полином, имеющий степень $0, 1, 2$ или $3$, формула Симпсона будет давать точный результат. В дальнейшем мы рассмотрим способ построения квадратур, имеющих наибольшую степень точности при наименьшем числе используемых узлов.

Составные формулы численного интегрирования

Как и в случае с интерполяцией, использование формул Ньютона-Котеса большого порядка для увелечения точности интегрирования на больших отрезках приводит к паразитным решениям. В таком случае предпочтительным является кусочный подход к интегрированию. Мы построим составные формулы из формул Симпсона, трапеций и средних, хотя подобный подход, очевидно, можно применить к любым другим квадратурам. Рассмотрим интеграл $\int_a^b f(x) dx$. Разделим отрезок $[a; b]$ на четное число подотрезков $n$ и применим формулу Симпсона на каждом из них. Тогда для $f(x) \in C^{4}[a;b]$, $h = \frac{b - a}{n}$ и $x_i = a + (i-1)h$, где $i = 1, \dots, n+1$ мы имеем: где $\xi_i \in (x_{2i - 1}; x_{2i + 1})$. Заметим, что все нечетные узлы за исключением $x_1$ и $x_{n+1}$ дважды повторяются в сумме, что позволяет упростить выражение: Для упрощения выражения для остаточного члена заметим, что так как $f(x) \in C^{4}[a; b]$, мы имеем: Тогда по теореме о промежуточном значении мы получаем: где $\xi \in (a; b)$. Остаточный член в таком случае принимает следующую форму: Заметим, что по сравнению со стандартной формулой Симпсона, остаточный член составной формулы Симпсона пропорционален $O(h^4)$. Результирующее выражение для составной формулы представлено в следующей теореме. Аналогично можно вывести составные формулы трапеций и средних, в которых $n$ может быть как четным, так и нечетным. Мы приведет соответствующие теоремы без доказательства.

Вычислительная устойчивость операции интегрирования

По сравнению с дифференцированием, операция интегрирования, как мы вскоре покажем, способна к стабилизации вычислительной погрешности. Интуитивное объяснение этого эффект заключается в том, что интегрирование предполагает суммирование близких значений, в то время как дифференцирование вычисляет их разность. Рассмотрим составную формулу Симпсона и предположим, что значение $f(x)$ в точке $x_i$ вычисляется с погрешностью округления $e_i$: Тогда полная погрешность округления, накапливаемая составной формулой Симпсона, может быть оценена следующим образом: Предположим, что погрешность округления ограничена, например, машинным эпсилон: $|e_i| \leq \epsilon$, $i = 1, \dots, n+1$. Тогда полная погрешность погрешность оценивается как: Этот результат показывает, что верхняя грань для накопленной погрешности округления не зависит от $n$ или $h$, что означает, что увеличение числа подотрезков не приводит к дестабилизации полной погрешности. Действительно, из графика можно заметить, что полная погрешность интегрирования падает до тех пор, пока она не достигнет значения, сравнимого с машинным эпсилон, после чего уменьшение погрешности становится невозможным, и она стабилизируется на уровне машинного эпсилон.

Квадратуры Гаусса

Формулы Ньютона-Котеса были построены, исходя из равномерного распределения узлов. Как и в случае интерполяции, логично предположить, что существует иное, неравномерное распределение узлов, которое позволило бы максимизировать точность. Как уже было сказано ранее, на практике часто возникает необходимость максимизировать степень точности, т.е. максимизировать степень полинома, интегрирование которого с помощью данной квадратуры дает точный результат при заданном количестве узлов. Эту задачу решают квадратуры Гаусса, суть которых сводится к нахождению таких $x_1, \dots, x_n$ и $c_1, \dots, c_n$, что приближение максимизирует степень точности. Так как всего мы имеем $2n$ оптимизируемых параметров, логично предположить, что полином $2n - 1$ степени, имеющий так же $2n$ параметров, может быть интегрирован точно при правильно подобранных параметрах. Очевидно, что если квадратура дает точный результат для любого полинома этой степени, то все полиномы низших степеней автоматически интегрируются точно как частные случаи. Рассмотрим случай $n=2$ и интервал интегрирования $[-1; 1]$. Тогда квадратура принимает вид: Мы ожидаем, что эта квадратура дает точный результат при интегрировании полинома третьей степени: где $a_0, a_1, a_2, a_3 \in \mathbb{R}$ – произвольные константы. Можно заметить, что квадратура будет точно вычислять интеграл этого полинома тогда, когда точно будут вычисляться интегралы от функций $f(x) = 1$, $f(x) = x$, $f(x) = x^2$, $f(x) = x^3$. В таком случае мы получаем систему уравнений: Из первых двух уравнений мы получаем: что при подстановке в третье дает: Тогда четвертое уравнение становится: что в результате дает следующее решение: Таким образом квадратура Гаусса со степенью точности 3 имеет следующий вид: Подобным образом мы можем построить квадратуры Гаусса произвольной степени точности, однако из-за нелинейности решаемой системы уравнений это становится алгебраически утомительным. Для рассмотрения более удобного подхода нам необходимо ознакомиться с понятием линейно-независимых функций и многочленами Лежандра.

Ортогональные полиномы и многочлены Лежандра

По аналогии с векторами мы введем понятие линейно независимых функций. Система функций $\{\phi_1, \dots, \phi_n\}$ называется линейно независимой на $[a; b]$, если

\begin{equation} \sum_{i=1}^{n} c_i \phi_1(x) = 0 \quad \Longleftrightarrow \quad c_1 = \dots = c_n = 0. \end{equation}

В противном случае система функций называется линейно зависимой. Важное свойство полиномов состоит в том, что любые полиномы различных степеней являются линейно независимыми. Это свойство доказывается в следующей теореме. Также добавим, что если система полиномов $\{\phi_i\}_{i=0}^{n}$ является линейно независимой, то любой полином степени меньшей или равной $n$ можно представить в виде единственной линейной комбинации полиномов $\{\phi_i\}_{i=0}^{n}$. Более того, если такая система полиномов является к тому же ортогональной, то верно следующее равенство: где $P_k(x)$ – полиномом любой степени $k < n$. Действительно, мы может разложить $P_k(x)$ в линейную комбинацию полиномов $\phi(x)$: что дает Линейно независимая система функций необязательно является ортогональной. Для случая полиномов мы можем построить ортогональную систему полиномов, адаптировав процесс Грама-Шмидта, хорошо известный из курса линейной алгебры. Многочлены Лежандра могут быть получены различными способами, однако мы рассмотрим их как ортогональную систему полиномов, построенную с помощью процесса Грама-Шмидта. Действительно, для случая $\omega(x) = 1$ и интервала $[-1; 1]$ мы имеем: Остальные многочлены Лежандра могут быть получены аналогичным образом. Наконец, мы имеем возможность доказать следующую теорему, которая постулирует, что квадратуры Гаусса могут быть построены, если в качестве узлов выбраны корни соответствующего многочлена Лежандра. Интегрирование на произвольном интервале $[a; b]$ с помощью квадратуры Гаусса реализуется с помощью замены переменных: где $x \in [a; b]$ и $t \in [-1; 1]$. Тогда квадратура Гаусса вычисляется следующим образом: