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

Основные авторы описания: А.Б.Кукаркин

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Задача одномерного численного интегрирования состоит в приближенном вычислении определенного интеграла по отрезку: требуется вычислить

[math]I=\int\limits_a^b f(x)\,dx,[/math]

где [math]f(x)[/math] — интегрируемая функция, определенная на отрезке [math][a,b][/math].

Численное интегрирование оcуществляется с помощью квадратурных формул — приближенных равенств вида

[math]I\approx\sum\limits_{i=1}^n C_i f(x_i).[/math]

Сумма в правой части называется квадратурной суммой; различные точки [math]x_i[/math] отрезка [math][a,b][/math] называются узлами, а числа [math]C_i[/math] — коэффициентами квадратурной формулы.

Значение квадратурной суммы, принимаемое за приближенное значение интеграла, зависит от выбора узлов [math]x_i[/math] и коэффициентов [math]C_i[/math]. При вычислении значения квадратурной суммы основные временные затраты приходятся на вычисление значений подынтегральной функции в узлах.

Погрешностью квадратурной формулы называется разность

[math]R_n(f)=\int\limits_a^b\!f(x)\,dx-\sum\limits_{i=1}^n C_if(x_i).[/math]

Если функция [math]f(x)[/math] такова, что [math]R_n(f)=0[/math], то говорят,что квадратурная формула точна для функции [math]f(x)[/math].

Целое число [math]k\ge0[/math] называется алгебраической степенью точности квадратурной формулы, если квадратурная формула точна для любого многочлена степени не выше [math]k[/math] и не точна для многочлена степени [math]k+1[/math].

Если квадратурная формула с [math]n[/math] узлами имеет алгебраическую степень точности [math]k[/math], то [math]k\le 2n-1[/math].

Широко известны квадратурные формулы, полученные посредством замены подынтегральной функции [math]f(x)[/math] алгебраическим интерполяционным многочленом, значения которого совпадают со значениями функции в [math]n[/math] узлах [math]x_i[/math]; такие квадратурные формулы называются интерполяционными. Для погрешности интерполяционной квадратурной формулы справедлива оценка

[math]|R_n(f)|\le\frac{M_n}{n!}\int\limits_a^b|\,\omega_n(x)|\,dx,[/math]

где [math]M_n=\max\limits_{x\in[a,b]}|f^{(n)}(x)|[/math], [math]\omega_n(x)=\prod\limits_{i=1}^n (x-x_i)[/math].

Если интерполяционная квадратурная формула с [math]n[/math] узлами имеет алгебраическую степень точности [math]k[/math], то [math]k\ge n-1[/math].

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

Формулами Ньютона-Котеса называются интерполяционные квадратурные формулы, [math]n[/math] узлов которых заданы равноотстоящими: [math]x_1=\frac{a+b}2[/math] при [math]n=1[/math] и [math]x_i=a+(i-1)\frac{b-a}{n-1}, 1\le i \le n[/math] при [math]n\gt 1[/math].

Свое название эти формулы получили в память того, что они в достаточно общей форме были рассмотрены Исааком Ньютоном, а их коэффициенты при [math]1\le n\le 10[/math] были найдены Роджером Котесом.

Наиболее известными формулами Ньютона-Котеса являются формула средних прямоугольников ([math]n=1[/math])

[math]I\approx(b-a)f\left(\frac{a+b}2\right)[/math]

и формула трапеций ([math]n=2[/math])

[math]I\approx\frac{b-a}2\,(f(a)+f(b));[/math]

алгебраическая степень точности каждой из них равна 1, а для их погрешностей справедливы оценки

[math]|R_1(f)|\le\frac{M_2}{24}\,(b-a)^3,\quad |R_2(f)|\le\frac{M_2}{12}\,(b-a)^3,[/math]

где [math]M_2=\max\limits_{x\in[a,b]}|f^{\prime\prime}(x)|[/math].

Формулами Ньютона-Котеса при [math]n=3[/math] и [math]n=4[/math] являются формула Симпсона

[math]I\approx\frac{b-a}6\,\left(f(a)+4f\left(\frac{a+b}2\right)+f(b)\right),[/math]

которую называют также формулой парабол, и формула трех восьмых

[math]I\approx\frac{b-a}{8}\,\left(f(a)+3f\left(\frac{2a+b}3\right)+3f\left(\frac{a+2b}3\right)+f(b)\right),[/math]

название которой происходит от коэффициента 3/8; алгебраическая степень точности каждой из них равна 3, а для их погрешностей справедливы оценки

[math]|R_3(f)|\le\frac{M_4}{2880}\,(b-a)^5,\quad |R_4(f)|\le\frac{M_4}{6480}\,(b-a)^5,[/math]

где [math]M_4=\max\limits_{x\in[a,b]}|f^{(4)}(x)|[/math].

Алгебраическая степень точности формулы Ньютона-Котеса с [math]n[/math] узлами равна [math]n-1[/math] при четном [math]n[/math] и равна [math]n[/math] при нечетном [math]n[/math].

Коэффициенты формул Ньютона-Котеса положительны при [math]1\le n\le8[/math] и [math]n=10[/math], а при [math]n=9[/math] и [math]n\ge11[/math] среди коэффициентов имеются как положительные так и отрицательные.

Положительность коэффициентов квадратурной формулы важна для ее практического применения. Дело в том, что при вычислении интегральной суммы влияние погрешностей округления на точность результата тем сильнее, чем больше [math]\sum\limits_{i=1}^n |C_i|[/math]. Для формул Ньютона-Котеса эта сумма неограниченно возрастает при [math]n\to\infty[/math]. Поэтому при больших [math]n[/math] формулы Ньютона-Котеса оказываются практически непригодными.

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

1.1.2 Составные квадратурные формулы

Составными называются квадратурные формулы, построение которых осуществляется следующим образом. Отрезок интегрирования [math][a,b][/math] разбивается на [math]m[/math] отрезков равной длины [math]h=\dfrac{b-a}m[/math]. Тогда по свойству аддитивности интеграл может быть вычислен как сумма интегралов по отрезкам разбиения. Для приближенного вычисления каждого из слагаемых этой суммы применяется одна и та же квадратурная формула, называемая исходной.

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

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

составная формула средних прямоугольников

[math]I\approx h\,(\sum\limits_{i=1}^m f(a+ih-\frac{h}{2}))[/math]

с числом узлов [math]m[/math] и оценкой погрешности

[math]|R_m(f)|\le\frac{M_2}{24}\,(b-a)h^2,[/math]

составная формула трапеций

[math]I\approx \frac{h}{2}\,(f(a)+f(b)+2\sum\limits_{i=1}^{m-1}f(a+ih))[/math]

с числом узлов [math]m+1[/math] и оценкой погрешности

[math]|R_{m+1}(f)|\le\frac{M_2}{12}\,(b-a)h^2,[/math]

составная формула Симпсона

[math]I\approx \frac{h}{6}\,(f(a)+f(b)+2\sum\limits_{i=1}^{m-1}f(a+ih)+4\sum\limits_{i=1}^m f(a+ih-\frac{h}{2}))[/math]

с числом узлов [math]2m+1[/math] и оценкой погрешности

[math]|R_{2m+1}(f)|\le\frac{M_4}{2880}\,(b-a)h^4.[/math]

1.1.3 Квадратурные Формулы Гаусса

Как отмечалось выше, если квадратурная формула с [math]n[/math] узлами является интерполяционной, то ее алгебраическая степень точности не меньше [math]n-1[/math]; при любых заданных узлах построение интерполяционной квадратурной формулы осуществляется за счет выбора ее [math]n[/math] коэффициентов. За счет же выбора [math]n[/math] узлов интерполяционной квадратурной формулы можно добиться того, чтобы она имела возможно более высокую алгебраическую степень точности, а именно [math]2n-1[/math].

Задача построения такой квадратурной формулы рассматривалась Карлом Фридрихом Гауссом; им была доказана ее разрешимость.

Квадратурная формула с [math]n[/math] узлами, алгебраическая степень точности которой равна [math]2n-1[/math], называется квадратурной формулой Гаусса или квадратурной формулой наивысшей алгебраической степени точности.

В случае отрезка [math][-1,1][/math] квадратурная формула

[math]\int\limits_{-1}^1 f(t)\,dt\approx\sum\limits_{i=1}^n c_i f(t_i)[/math]

является квадратурной формулой Гаусса тогда и только тогда, когда она является интерполяционной, а ее узлы [math]t_i[/math] являются корнями многочлена Лежандра

[math]P_n(t)=\frac1{2^nn!}\frac{d^n}{dx^n}((x^2-1)^n).[/math]

В частности, при [math]n=2[/math] и [math]n=3[/math] квадратурные формулы Гаусса для отрезка [math][-1,1][/math] таковы:

[math]\int\limits_{-1}^1 f(t)\,dt\approx f\left(-\frac1{\sqrt3}\right)+f\left(\frac1{\sqrt3}\right),[/math]

[math]\int\limits_{-1}^1 f(t)\,dt\approx\frac59\;f\left(-\sqrt{\frac35}\right)+\frac89\;f(0)+ \frac59\;f\left(\sqrt{\frac35}\right).[/math]

Чтобы получить квадратурную формулу Гаусса для произвольного отрезка [math][a,b][/math], следует сделать замену переменной

[math]x=\frac{a+b}2+\frac{b-a}2\,t,[/math]

в результате которой

[math]\int\limits_a^b\!f(x)\,dx=\frac{b-a}2\int\limits_{-1}^1\!f\left(\frac{a+b}2+\frac{b-a}2\,t\right)\,dt.[/math]

Воспользовавшись квадратурной формулой Гаусса для отрезка [math][-1,1][/math], получим квадратурную формулу Гаусса для отрезка [math][a,b][/math]:

[math]I\approx\frac{b-a}2\sum\limits_{i=1}^n c_i f(x_i),[/math]

где [math]x_i=0.5(a+b+(b-a)t_i)[/math], [math]t_i[/math] — корни многочлена Лежандра [math]P_n(t)[/math].

Для погрешности квадратурной формулы Гаусса c [math]n[/math] узлами справедлива оценка

[math]|R_n(f)|\le\frac{M_{2n}(n!)^4}{((2n)!)^3(2n+1)}(b-a)^{2n+1},[/math]

где [math]M_{2n}=\max\limits_{x\in[a,b]}|f^{(2n)}(x)|[/math].

Коэффициенты квадратурных формул Гаусса положительны. Поэтому использование квадратурных формул Гаусса с большим числом узлов не приводит к тем осложнениям, которые возникают при использовании формул Ньютона-Котеса.

1.2 Математическое описание алгоритма

Если квадратурная формула выбрана, то алгоритм приближенного вычисления интеграла состоит в вычислении квадратурной суммы, т.е. в вычислении значений функции в [math]n[/math] узлах, умножении их на соответствующие коэффициенты и суммировании полученных чисел.

Исходные данные: функция [math]f(x)[/math] и два одномерных массива [math]n[/math] чисел — массив узлов и массив коэффициентов.

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