5.3. Формулы Гаусса

В Разделе 5.1 мы видели, что веса квадратуры \(\{w_i\}\) с известными узлами \(\{x_i\}\) могут быть найдены решением линейной системы уравнений

\[\sum\limits_{i=1}^{N} w_i x_i^p = \frac{1}{p + 1} \quad\forall p \in \{0, \dots, N-1\}.\]

Хотя она была записана для формул Ньютона-Котса, её можно решить относительно весов для любого набора различных узлов \(\{x_i\}\), не обязательно равномерно-распределённых.

Но эти же условия можно рассматривать как условия на веса и узлы квадратуры, т.е. рассмотреть нелинейную систему уравнений

(5.19)\[\sum\limits_{i=1}^{N} w_i x_i^p = \frac{1}{p + 1}, \quad\forall p \in \{0, \dots, M\}\]

c \(2N\) неизвестными: \(\{w_i\}_{i=1}^N\) и \(\{x_i\}_{i=1}^N\). Тогда можно надеяться, что система будет разрешима вплоть до \(M = 2N-1\). Если это так, то мы получим квадратурную формулу с \(N\) узлами, точную для многочленов вплоть до степени \(2N-1\).

Решение системы «в лоб» из-за её нелинейности проблематично, но оказывается, что для её решения есть элегантный подход.

Вывод квадратуры Гаусса

Рассмотрим интегрирование функций на отрезке \([-1, 1]\). Введём для квадратурной формулы с узлами \(\{x_1, x_2, \dots, x_N\}\) многочлен степени \(N\) с нулями в этих узлах

\[\Omega_{N}(x) = \text{const} \cdot (x - x_1)(x - x_2)\dots(x - x_N).\]

Выразим через \(\Omega_{N}\) произвольный многочлен \(P_M\) степени \(M \ge N\) так

\[P_M(x) = \Omega_{N}(x) q_{M-N}(x) + r_{N-1}(x),\]

где \(q_{M-N}\) – многочлен степени \(M - N\), \(r_{N-1}\) – многочлен степени не выше \(N - 1\).

Запишем ошибку интегрирования \(P_M\) при условии, что веса квадратуры определены \(w_i\) из системы (5.4), т.е. квадратура точна для всех многочленов вплоть до степени \(N-1\):

(5.20)\[\begin{split}\begin{split} \sum\limits_{i} w_i P_M(x_i) &- \int\limits_{-1}^{1} P_M(x) \diff x = \\ & = \left(\sum\limits_{i} w_i \Omega_{N}(x_i) q_{M-N}(x_i) - \int\limits_{-1}^{1} \Omega_{N}(x) q_{M-N}(x) \diff x \right) + \left(\sum\limits_{i} w_i r_{N-1}(x_i) - \int\limits_{-1}^{1} r_{N-1}(x) \diff x \right) \end{split}\end{split}\]

Вторая скобка равна нулю благодаря выбору весов. Сумма в первой скобке равна нулю, т.к. по своему определению \(\Omega_{N}(x_i) = 0\). Остаётся лишь интеграл в первой скобке. Отсюда получаем вывод: ошибка интегрирования равна нулю, если

(5.21)\[\int\limits_{-1}^{1} \Omega_{N}(x) q_{M-N}(x) \diff x = 0,\quad \forall q_{M-N}(x).\]

Т.е. многочлен \(\Omega_{N}(x)\) ортогонален всем многочленам вплоть до степени \(M - N\), если скалярное произведение определено как

\[f \cdot g = \int\limits_{-1}^{1} f(x) g(x) \diff x.\]

Как известно из матанализа, для многочленов ортонормированный базис составляют многочлены Лежандра

\[L_k(x) = \frac{1}{2^k\ k!} \frac{d^k}{dx^k}(x^2 - 1)^k.\]

Значит, в качестве \(\Omega_N\) возьмём \(L_N\), чтобы удовлетворить (5.21) для полиномов \(q_{M-N}\) как можно более высокой степени. Тогда, если степень \(q_{M-N}\) не превосходит \(M-N \le N-1\), т.е. \(M \le 2N-1\), то (5.21) выполнится, т.е.

\[\int\limits_{-1}^{1} \Omega_{N}(x) q_{M-N}(x) \diff x = 0,\quad \forall M \le 2N-1,\]

а для \(M > 2N-1\) получим наименьшее скалярное произведение.

Поскольку \(\Omega_N \equiv L_N\), то узлы интегрирования совпадают с нулями \(L_N(x_i) = 0\), \(i=1,\dots,N\). Зная узлы, мы решаем систему на веса \(\{w_i\}_{i=1}^N\) (5.4) из \(N\) уравнений, получая квадратурную формулу.

Из выражения для ошибки (5.20) и такого способа выбора узлов и весов приходим к выводу, что итоговая квадратурная формула будет точна для полиномов \(P_M\) степени не выше \(M \le 2N-1\). Хотя узлов мы используем всего \(N\).

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

Определение 5.13 (Квадратурная формула Гаусса)

Квадратурные формулы с узлами в нулях многочленов Лежандра называются квадратурными формулами Гаусса. Они имеют при заданном числе узлов (т.е. количестве вычислений подынтегральной функции) максимальный алгебраический порядок точности из всех квадратур вида (5.2).

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

1. Узлы и веса квадратур Гаусса на отрезке \([-1, 1]\) симметричны относительно нуля, поэтому в таблицах чаще всего приводится только половина их;

2. В некоторых таблицах приводятся узлы и веса на отрезке \([-1, 1]\), в некоторых – на \([0, 1]\). В зависимости от этого формулы пересчёта на произвольный отрезок \([a, b]\) будут слегка отличаться.

5.3.1. Квадратуры гауссова типа с частично зафиксированными узлами

На практике иногда возникает потребность получить максимально точные квадратурные формулы при условии, что часть узлов \(\{{\hat x}_1, \dots, {\hat x}_p\}\)зафиксирована, а остальные \(N - p\) узлов можно выбрать свободно. В таком случае те же рассуждения, что и при выводе формул Гаусса, приводят к следующему условию максимальной степени квадратуры:

(5.23)\[\int\limits_{-1}^{1} {\hat\Omega}_p(x) \Omega_{N-p}(x) q_{M-N}(x) \diff x = 0\]

5.3.1.1. Квадратуры Лобатто

Гауссовы квадратуры относятся к открытому типу, т.к. все нули многочленов Лежандра лежат внутри интервала \((-1, 1)\). Если нужно построить формулу закрытого типа, то положения \(x_1 = -1\) и \(x_N = 1\) фиксируются. Выражение (5.23) приобретает вид

(5.24)\[\int\limits_{-1}^{1} (1 - x^2) \Omega_{N-2}(x) q_{M-N}(x) \diff x = 0\]

Поскольку функция \({\hat\Omega}_p(x) = (1 - x^2)\) положительна на интервале \((-1, 1)\), выражение (5.24) можно рассматривать как скалярное произведение, а \(\Omega_{N-2}(x)\) должен принадлежать семейству ортогональных многочленов для такого скалярного произведения. Поиск коэффициентов таких многочленов можно производить через ортогонализацию Грама-Шмидта. Когда многочлен найден, каким-либо численным методом находятся его нули, которые будут дополнительными узлами, а затем веса.

Для квадратур Лобатто дополнительные узлы определяются нулями \(L_{N-1}'(x)\). Алгебраический порядок точности квадратуры равен \(2N - 3\) (выбор \(N - 2\) узлов и \(N\) весов позволяют получить точные решения для многочленов с \(2N-2\) коэффициентами).

5.3.1.2. Квадратуры Кронрода

Гауссовы квадратуры обладают неприятным свойством – при повышении порядка аппроксимации полностью меняется набор узлов интегрирования. Это означает, что оценка фактической погрешности численного интегрирования требует расчёта подынтегральной функции в \(N+1\) точке, при этом уточнённая квадратура будет иметь порядок лишь на 2 выше, чем исходная. В случае, например, интегрирования методом Ромберга аналогичное удвоение количества точек удваивает и порядок квадратуры.

В 1960-х Александр Семёнович Кронрод предложил «уточняющие» квадратуры, получаемые добавлением к \(N\) узлам гауссовой квадратуры ещё \(N+1\) узлов так, чтобы порядок получившейся квадратуры был максимален.

Поскольку в (5.23) теперь берётся не всюду положительный многочлен \({\hat\Omega}_p(x) = L_p(x)\), этот интеграл уже нельзя рассматривать как скалярное произведение, и поиск дополняющего многочлена нужно делать иначе.

Записав многочлен в виде

(5.25)\[K(x) = x^{N+1} + K_N x^N + \dots + K_1 x + K_0,\]

и подставив в (5.23) \(q(x) = 1\), получаем

\[\int\limits_{-1}^{1} L_N(x) (x^{N+1} + K_N x^N + \dots) \diff x = \int\limits_{-1}^{1} L_N(x) (x^{N+1} + K_N x^N) \diff x = 0,\]

т.к. многочлен Лежандра \(L_N(x)\) ортогонален всем многочленам степени ниже \(N\). Отсюда

\[K_N = \frac{\int\limits_{-1}^{1} L_N(x) x^{N+1} \diff x}{\int\limits_{-1}^{1} L_N(x) x^N \diff x}.\]

По аналогии, подставив \(q(x) = x\), получаем

\[\int\limits_{-1}^{1} L_N(x) (x^{N+2} + K_N x^{N+1} + K_{N-1} x^N) \diff x = 0,\]

откуда уже можно найти \(K_{N-1}\) и т.д.

Квадратура с \(N\) гауссовыми узлами и \(N+1\) узлами Кронрода точна для многочленов степени до \(3N\) при чётных \(N\) и \(3N + 1\) при нечётных \(N\) (как и для квадратур Ньютона-Котса, при нечётное число узлов даёт дополнительный порядок за счёт симметрии).

5.3.2. Альтернативная система уравнений для вычисления весов квадратур

При выводе системы (5.4) требовалось точное выполнение квадратурной формулы для степеней \(x\). Можно вместо этого потребовать точное выполнение квадратурной формулы для многочленов Лежандра. Система в предположении, что узлы заданы на отрезке \([-1, 1]\), принимает вид

\[\begin{split}\left( \begin{array}{cccc} 1 & 1 & \dots & 1 \\ x_1 & x_2 & \dots & x_N \\ L_2(x_1) & L_2(x_2) & \dots & L_2(x_N) \\ & & \ddots & \\ L_{N-1}(x_1) & L_{N-1}(x_2) & \dots & L_{N-1}(x_N) \end{array} \right) \left( \begin{array}{c} w_1 \\ w_2 \\ w_3 \\ \vdots \\ w_N \end{array} \right) = \left( \begin{array}{c} 1 \\ 0 \\ 0 \\ \vdots \\ 0 \end{array} \right)\end{split}\]

Правая часть состоит в основном из нулей, т.к. \(\int\limits_{-1}^{1} L_M(x) \diff x = \int\limits_{-1}^{1} L_M(x) L_0(x)\diff x = 0\) при всех \(M > 0\).

Учитывая, что

\[L_{M+1}(x) = \frac{2M + 1}{M + 1} x L_M(x) - \frac{M}{M + 1} L_{M-1}(x),\]

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

5.3.3. Упражения

1. Ниже приведена заготовка функции, находящей узлы и коэффициенты формулы Ньютона-Котса, точной для многочлена степени \(p\). Заполните пропуски ⟨...⟩

   """
   Находит узлы и веса формулы Ньютона-Котса с числом узлов `nnodes` на отрезке `[0; 1]`.
   Если аргумент `closed=true` (по умолчанию), то расчет делается для 
   закрытых квадратурных формул, иначе - для открытых
   """
   function newton_cotes_nodes_weights(nnodes::Integer; closed::Bool=true)
       if closed
           nodes = range(⟨...⟩)
       else
           nodes = range(⟨...⟩)
       end
       rhp = zero.(nnodes)
       rhp[1] = 1
       w = legendre_matrix(nodes) \ rhp
       return (nodes=nodes, weights=w)
   end
   
   function legendre_matrix(xs::AbstractVector)
       n = length(xs)
       m = similar(xs, Float64, (n, n))
       for c in 1:n
           m[1,c] = 1
           ⟨...⟩
           for r in 2:n
               ⟨...⟩
           end
       end
       return m
   end
   

2. Вычислите численное приближение для интеграла

   \[\int\limits_{-5}^{5} \frac{\diff x}{1 + x^2}\]

   при помощи формул Ньютона-Котса различных порядков. Сходятся ли значения к точному значению \(2 \text{arctg}(5) = 2.7468015338900317\dots\)?

3. Вычислите численное приближение для интеграла

   \[\int\limits_{-5}^{5} \frac{\diff x}{1 + x^2}\]

   при помощи формул Гаусса по 7 точкам (\(G_7\)) и Гаусса-Кронрода по 15 точкам (\(K_{15}\)). Даёт ли величина \(\vert G_7 - K_{15}\vert\) хорошую оценку погрешности численного интегрирования?