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

Приближение определенного интеграла функции

Сравнение между 2-точечной гауссовой и трапециевидной квадратурой.
Синяя кривая показывает функцию, определенный интеграл которой на интервале [−1, 1] должен быть вычислен (интегральное выражение). Формула трапеции аппроксимирует функцию линейной функцией, которая совпадает с подынтегральным выражением на концах интервала и представлена ​​оранжевой пунктирной линией. Аппроксимация, по-видимому, не очень хорошая, поэтому ошибка велика ( формула трапеции дает аппроксимацию интеграла, равную y (–1) + y (1) = –10 , тогда как правильное значение равно 2 ⁄ 3 ). Чтобы получить более точный результат, интервал необходимо разбить на множество подынтервалов, а затем использовать составную
формулу трапеции, что требует гораздо больше вычислений. Вместо этого гауссова квадратура выбирает более подходящие точки, поэтому даже линейная функция аппроксимирует функцию лучше (черная пунктирная линия). Поскольку подынтегральное выражение представляет собой многочлен 3-й степени ( y ( x ) = 7 x ³ – 8 x ² – 3 x + 3 ), двухточечное квадратурное уравнение Гаусса даже возвращает точный результат.

В численном анализе квадратурное правило Гаусса для n точек , названное в честь Карла Фридриха Гаусса , ^[1] представляет собой квадратурное правило, построенное для получения точного результата для многочленов степени 2n − 1 или меньше путем подходящего выбора узлов x _i и весов w _i для i = 1, ..., n .

Современная формулировка с использованием ортогональных многочленов была разработана Карлом Густавом Якоби в 1826 году. ^[2] Наиболее распространенная область интегрирования для такого правила берется как [−1, 1] , поэтому правило формулируется как $${\displaystyle \int _{-1}^{1}f(x)\,dx\approx \sum _{i=1}^{n}w_{i}f(x_{i}),}$$

что является точным для многочленов степени 2 n − 1 или меньше. Это точное правило известно как квадратурное правило Гаусса–Лежандра . Квадратурное правило будет точным приближением к интегралу выше только в том случае, если f ( x ) хорошо приближается многочленом степени 2 n − 1 или меньше на [−1, 1] .

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

$${\displaystyle f(x)=\left(1-x\right)^{\alpha }\left(1+x\right)^{\beta }g(x),\quad \alpha ,\beta >-1,}$$

где g ( x ) хорошо аппроксимируется полиномом низкой степени, то альтернативные узлы x _i ' и веса w _i ' обычно дают более точные квадратурные правила. Они известны как квадратурные правила Гаусса-Якоби , т.е.

$${\displaystyle \int _{-1}^{1}f(x)\,dx=\int _{-1}^{1}\left(1-x\right)^{\alpha }\left(1+x\right)^{\beta }g(x)\,dx\approx \sum _{i=1}^{n}w_{i}'g\left(x_{i}'\right).}$$

Обычные веса включают ( Чебышева–Гаусса ) и . Также может потребоваться интегрировать по полубесконечным ( квадратура Гаусса–Лагерра ) и бесконечным интервалам ( квадратура Гаусса–Эрмита ). ${\textstyle {\frac {1}{\sqrt {1-x^{2}}}}}$ ${\textstyle {\sqrt {1-x^{2}}}}$

Можно показать (см. Press et al. или Stoer and Bulirsch), что квадратурные узлы x _i являются корнями многочлена, принадлежащего классу ортогональных многочленов (классу, ортогональному относительно взвешенного внутреннего произведения). Это ключевое наблюдение для вычисления квадратурных узлов и весов Гаусса.

Квадратура Гаусса–Лежандра

Графики полиномов Лежандра (до n = 5)

Для простейшей задачи интеграции, указанной выше, то есть, f ( x ) хорошо аппроксимируется полиномами от , соответствующие ортогональные полиномы являются полиномами Лежандра , обозначаемыми как P _n ( x ) . С n -м полиномом, нормализованным так, чтобы дать P _n (1) = 1 , i -й узел Гаусса, x _i , является i -м корнем P _n , а веса задаются формулой ^[3] ${\displaystyle [-1,1]}$ $${\displaystyle w_{i}={\frac {2}{\left(1-x_{i}^{2}\right)\left[P'_{n}(x_{i})\right]^{2}}}.}$$

Некоторые квадратурные правила низкого порядка приведены в таблице ниже (на интервале [−1, 1] , другие интервалы см. в разделе ниже).

Изменение интервала

Интеграл по [ a , b ] должен быть преобразован в интеграл по [−1, 1] перед применением квадратурного правила Гаусса. Это изменение интервала можно выполнить следующим образом: $${\displaystyle \int _{a}^{b}f(x)\,dx=\int _{-1}^{1}f\left({\frac {b-a}{2}}\xi +{\frac {a+b}{2}}\right)\,{\frac {dx}{d\xi }}d\xi }$$

с ${\displaystyle {\frac {dx}{d\xi }}={\frac {b-a}{2}}}$

Применение квадратурного правила Гаусса приводит к следующему приближению: ${\displaystyle n}$ ${\displaystyle (\xi ,w)}$ $${\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {b-a}{2}}\sum _{i=1}^{n}w_{i}f\left({\frac {b-a}{2}}\xi _{i}+{\frac {a+b}{2}}\right).}$$

Пример двухточечной квадратурной формулы Гаусса

Используйте квадратурное правило Гаусса для двух точек, чтобы приблизительно рассчитать расстояние в метрах, пройденное ракетой от до, как указано в формуле ${\displaystyle t=8\mathrm {s} }$ ${\displaystyle t=30\mathrm {s} ,}$ $${\displaystyle s=\int _{8}^{30}{\left(2000\ln \left[{\frac {140000}{140000-2100t}}\right]-9.8t\right){dt}}}$$

Измените пределы так, чтобы можно было использовать веса и абсциссы, указанные в Таблице 1. Также найдите абсолютную относительную истинную ошибку. Истинное значение дано как 11061,34 м.

Решение

Во-первых, изменение пределов интегрирования с на дает ${\displaystyle \left[8,30\right]}$ ${\displaystyle \left[-1,1\right]}$

$${\displaystyle {\begin{aligned}\int _{8}^{30}{f(t)dt}&={\frac {30-8}{2}}\int _{-1}^{1}{f\left({\frac {30-8}{2}}x+{\frac {30+8}{2}}\right){dx}}\\&=11\int _{-1}^{1}{f\left(11x+19\right){dx}}\end{aligned}}}$$

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

• ${\displaystyle c_{1}=1.000000000}$
• ${\displaystyle x_{1}=-0.577350269}$
• ${\displaystyle c_{2}=1.000000000}$
• ${\displaystyle x_{2}=0.577350269}$

Теперь мы можем использовать квадратурную формулу Гаусса, поскольку $${\displaystyle {\begin{aligned}11\int _{-1}^{1}{f\left(11x+19\right){dx}}&\approx 11\left[c_{1}f\left(11x_{1}+19\right)+c_{2}f\left(11x_{2}+19\right)\right]\\&=11\left[f\left(11(-0.5773503)+19\right)+f\left(11(0.5773503)+19\right)\right]\\&=11\left[f(12.64915)+f(25.35085)\right]\\&=11\left[(296.8317)+(708.4811)\right]\\&=11058.44\end{aligned}}}$$ $${\displaystyle {\begin{aligned}f(12.64915)&=2000\ln \left[{\frac {140000}{140000-2100(12.64915)}}\right]-9.8(12.64915)\\&=296.8317\end{aligned}}}$$ $${\displaystyle {\begin{aligned}f(25.35085)&=2000\ln \left[{\frac {140000}{140000-2100(25.35085)}}\right]-9.8(25.35085)\\&=708.4811\end{aligned}}}$$

Учитывая, что истинное значение равно 11061,34 м, абсолютная относительная истинная погрешность составляет ${\displaystyle \left|\varepsilon _{t}\right|}$ $${\displaystyle \left|\varepsilon _{t}\right|=\left|{\frac {11061.34-11058.44}{11061.34}}\right|\times 100\%=0.0262\%}$$

Другие формы

Проблему интегрирования можно выразить немного более общим образом, введя положительную весовую функцию ω в подынтегральное выражение и допустив интервал, отличный от [−1, 1] . То есть, задача состоит в вычислении для некоторых выборов a , b и ω . Для a = −1 , b = 1 и ω ( x ) = 1 задача та же, что и рассмотренная выше. Другие выборы приводят к другим правилам интегрирования. Некоторые из них приведены в таблице ниже. Номера уравнений даны для Абрамовица и Стегуна (A & S). $${\displaystyle \int _{a}^{b}\omega (x)\,f(x)\,dx}$$

Основная теорема

Пусть p _n — нетривиальный многочлен степени n такой, что $${\displaystyle \int _{a}^{b}\omega (x)\,x^{k}p_{n}(x)\,dx=0,\quad {\text{for all }}k=0,1,\ldots ,n-1.}$$

Обратите внимание, что это будет верно для всех ортогональных многочленов выше, поскольку каждый p _n построен так, чтобы быть ортогональным другим многочленам p _j для j < n , а x ^k находится в пределах этого множества.

Если мы выберем n узлов x _i в качестве нулей p _n , то существуют n весов w _i , которые делают вычисленный гауссовой квадратурой интеграл точным для всех полиномов h ( x ) степени 2 n − 1 или меньше. Более того, все эти узлы x _i будут лежать в открытом интервале ( a , b ) . ^[4]

Чтобы доказать первую часть этого утверждения, пусть h ( x ) будет любым многочленом степени 2 n − 1 или меньше. Разделим его на ортогональный многочлен p _n , чтобы получить , где q ( x ) — частное степени n − 1 или меньше (потому что сумма его степени и степени делителя p _n должна быть равна степени делимого), а r ( x ) — остаток, также степени n − 1 или меньше (потому что степень остатка всегда меньше степени делителя). Поскольку p _n по предположению ортогонален всем одночленам степени меньше n , он должен быть ортогонален частному q ( x ) . Следовательно $${\displaystyle h(x)=p_{n}(x)\,q(x)+r(x).}$$ $${\displaystyle \int _{a}^{b}\omega (x)\,h(x)\,dx=\int _{a}^{b}\omega (x)\,{\big (}\,p_{n}(x)q(x)+r(x)\,{\big )}\,dx=\int _{a}^{b}\omega (x)\,r(x)\,dx.}$$

Поскольку остаток r ( x ) имеет степень n − 1 или меньше, мы можем интерполировать его точно, используя n точек интерполяции с полиномами Лагранжа l _i ( x ) , где $${\displaystyle l_{i}(x)=\prod _{j\neq i}{\frac {x-x_{j}}{x_{i}-x_{j}}}.}$$

У нас есть $${\displaystyle r(x)=\sum _{i=1}^{n}l_{i}(x)\,r(x_{i}).}$$

Тогда его интеграл будет равен $${\displaystyle \int _{a}^{b}\omega (x)\,r(x)\,dx=\int _{a}^{b}\omega (x)\,\sum _{i=1}^{n}l_{i}(x)\,r(x_{i})\,dx=\sum _{i=1}^{n}\,r(x_{i})\,\int _{a}^{b}\omega (x)\,l_{i}(x)\,dx=\sum _{i=1}^{n}\,r(x_{i})\,w_{i},}$$

где w _i , вес, связанный с узлом x _i , определяется как равный взвешенному интегралу l _i ( x ) (см. ниже другие формулы для весов). Но все x _i являются корнями p _n , поэтому формула деления выше говорит нам, что для всех i . Таким образом, мы в конечном итоге имеем $${\displaystyle h(x_{i})=p_{n}(x_{i})\,q(x_{i})+r(x_{i})=r(x_{i}),}$$ $${\displaystyle \int _{a}^{b}\omega (x)\,h(x)\,dx=\int _{a}^{b}\omega (x)\,r(x)\,dx=\sum _{i=1}^{n}w_{i}\,r(x_{i})=\sum _{i=1}^{n}w_{i}\,h(x_{i}).}$$

Это доказывает, что для любого многочлена h ( x ) степени 2 n − 1 или меньше его интеграл в точности равен сумме квадратур Гаусса.

Чтобы доказать вторую часть утверждения, рассмотрим факторизованную форму многочлена p _n . Любые комплексно-сопряженные корни дадут квадратный множитель, который либо строго положителен, либо строго отрицателен на всей действительной прямой. Любые множители для корней вне интервала от a до b не изменят знака на этом интервале. Наконец, для множителей, соответствующих корням x _i внутри интервала от a до b , которые имеют нечетную кратность, умножьте p _n еще на один множитель, чтобы получить новый многочлен $${\displaystyle p_{n}(x)\,\prod _{i}(x-x_{i}).}$$

Этот многочлен не может менять знак на интервале от a до b , поскольку все его корни там теперь имеют четную кратность. Поэтому интеграл, поскольку весовая функция ω ( x ) всегда неотрицательна. Но p _n ортогонален всем многочленам степени n -1 или меньше, поэтому степень произведения должна быть не менее n . Следовательно, p _n имеет n различных корней, все действительные, на интервале от a до b . $${\displaystyle \int _{a}^{b}p_{n}(x)\,\left(\prod _{i}(x-x_{i})\right)\,\omega (x)\,dx\neq 0,}$$ $${\displaystyle \prod _{i}(x-x_{i})}$$

Общая формула для весов

Веса можно выразить как

где — коэффициент в . Чтобы доказать это, отметим, что с помощью интерполяции Лагранжа можно выразить r ( x ) через as , поскольку r ( x ) имеет степень меньше n и, таким образом, фиксируется значениями, которых она достигает в n различных точках. Умножение обеих сторон на ω ( x ) и интегрирование от a до b дает ${\displaystyle a_{k}}$ ${\displaystyle x^{k}}$ ${\displaystyle p_{k}(x)}$ ${\displaystyle r(x_{i})}$ $${\displaystyle r(x)=\sum _{i=1}^{n}r(x_{i})\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}{\frac {x-x_{j}}{x_{i}-x_{j}}}}$$ $${\displaystyle \int _{a}^{b}\omega (x)r(x)dx=\sum _{i=1}^{n}r(x_{i})\int _{a}^{b}\omega (x)\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}{\frac {x-x_{j}}{x_{i}-x_{j}}}dx}$$

Веса w _i таким образом определяются как $${\displaystyle w_{i}=\int _{a}^{b}\omega (x)\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}{\frac {x-x_{j}}{x_{i}-x_{j}}}dx}$$

Это интегральное выражение для можно выразить через ортогональные многочлены следующим образом. ${\displaystyle w_{i}}$ ${\displaystyle p_{n}(x)}$ ${\displaystyle p_{n-1}(x)}$

Мы можем написать $${\displaystyle \prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}\left(x-x_{j}\right)={\frac {\prod _{1\leq j\leq n}\left(x-x_{j}\right)}{x-x_{i}}}={\frac {p_{n}(x)}{a_{n}\left(x-x_{i}\right)}}}$$

где - коэффициент в . Принимая предел x к получаем с помощью правила Лопиталя ${\displaystyle a_{n}}$ ${\displaystyle x^{n}}$ ${\displaystyle p_{n}(x)}$ ${\displaystyle x_{i}}$ $${\displaystyle \prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}\left(x_{i}-x_{j}\right)={\frac {p'_{n}(x_{i})}{a_{n}}}}$$

Таким образом, мы можем записать интегральное выражение для весов как

В подынтегральном выражении записываем $${\displaystyle {\frac {1}{x-x_{i}}}={\frac {1-\left({\frac {x}{x_{i}}}\right)^{k}}{x-x_{i}}}+\left({\frac {x}{x_{i}}}\right)^{k}{\frac {1}{x-x_{i}}}}$$

урожайность $${\displaystyle \int _{a}^{b}\omega (x){\frac {x^{k}p_{n}(x)}{x-x_{i}}}dx=x_{i}^{k}\int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx}$$

при условии , что является многочленом степени k − 1 , который затем ортогонален . Таким образом, если q ( x ) является многочленом не более n-й степени, то мы имеем ${\displaystyle k\leq n}$ $${\displaystyle {\frac {1-\left({\frac {x}{x_{i}}}\right)^{k}}{x-x_{i}}}}$$ ${\displaystyle p_{n}(x)}$ $${\displaystyle \int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx={\frac {1}{q(x_{i})}}\int _{a}^{b}\omega (x){\frac {q(x)p_{n}(x)}{x-x_{i}}}dx}$$

Мы можем оценить интеграл в правой части следующим образом. Поскольку является многочленом степени n − 1 , мы имеем где s ( x ) является многочленом степени . Поскольку s ( x ) ортогонален, мы имеем ${\displaystyle q(x)=p_{n-1}(x)}$ ${\displaystyle {\frac {p_{n}(x)}{x-x_{i}}}}$ $${\displaystyle {\frac {p_{n}(x)}{x-x_{i}}}=a_{n}x^{n-1}+s(x)}$$ ${\displaystyle n-2}$ ${\displaystyle p_{n-1}(x)}$ $${\displaystyle \int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx={\frac {a_{n}}{p_{n-1}(x_{i})}}\int _{a}^{b}\omega (x)p_{n-1}(x)x^{n-1}dx}$$

Затем мы можем написать $${\displaystyle x^{n-1}=\left(x^{n-1}-{\frac {p_{n-1}(x)}{a_{n-1}}}\right)+{\frac {p_{n-1}(x)}{a_{n-1}}}}$$

Член в скобках — это многочлен степени , который, следовательно, ортогонален . Таким образом, интеграл можно записать как ${\displaystyle n-2}$ ${\displaystyle p_{n-1}(x)}$ $${\displaystyle \int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx={\frac {a_{n}}{a_{n-1}p_{n-1}(x_{i})}}\int _{a}^{b}\omega (x)p_{n-1}(x)^{2}dx}$$

Согласно уравнению ( 2 ), веса получаются путем деления этого значения на , что дает выражение в уравнении ( 1 ). ${\displaystyle p'_{n}(x_{i})}$

${\displaystyle w_{i}}$также может быть выражено через ортогональные многочлены и теперь . В 3-членном рекуррентном соотношении член с исчезает, поэтому в уравнении (1) его можно заменить на . ${\displaystyle p_{n}(x)}$ ${\displaystyle p_{n+1}(x)}$ ${\displaystyle p_{n+1}(x_{i})=(a)p_{n}(x_{i})+(b)p_{n-1}(x_{i})}$ ${\displaystyle p_{n}(x_{i})}$ ${\displaystyle p_{n-1}(x_{i})}$ ${\textstyle {\frac {1}{b}}p_{n+1}\left(x_{i}\right)}$

Доказательство того, что веса положительны

Рассмотрим следующий многочлен степени , где, как и выше, x _j являются корнями многочлена . Очевидно , . Поскольку степень меньше , применяется квадратурная формула Гаусса, включающая веса и узлы, полученные из . Поскольку для j не равно i , имеем ${\displaystyle 2n-2}$ $${\displaystyle f(x)=\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}{\frac {\left(x-x_{j}\right)^{2}}{\left(x_{i}-x_{j}\right)^{2}}}}$$ ${\displaystyle p_{n}(x)}$ ${\displaystyle f(x_{j})=\delta _{ij}}$ ${\displaystyle f(x)}$ ${\displaystyle 2n-1}$ ${\displaystyle p_{n}(x)}$ ${\displaystyle f(x_{j})=0}$ $${\displaystyle \int _{a}^{b}\omega (x)f(x)dx=\sum _{j=1}^{n}w_{j}f(x_{j})=\sum _{j=1}^{n}\delta _{ij}w_{j}=w_{i}>0.}$$

Поскольку обе функции и неотрицательны, то отсюда следует, что . ${\displaystyle \omega (x)}$ ${\displaystyle f(x)}$ ${\displaystyle w_{i}>0}$

Вычисление квадратурных правил Гаусса

Существует множество алгоритмов для вычисления узлов x _i и весов w _i квадратурных правил Гаусса. Наиболее популярными являются алгоритм Голуба-Велша, требующий O ( n ² ) операций, метод Ньютона для решения с использованием трехчленной рекуррентности для оценки, требующий O ( n ² ) операций, и асимптотические формулы для больших n, требующие O ( n ) операций. ${\displaystyle p_{n}(x)=0}$

Рекуррентное соотношение

Ортогональные многочлены с для скалярного произведения , степенью и старшим коэффициентом один (т.е. монические ортогональные многочлены) удовлетворяют рекуррентному соотношению ${\displaystyle p_{r}}$ ${\displaystyle (p_{r},p_{s})=0}$ ${\displaystyle r\neq s}$ ${\displaystyle (\cdot ,\cdot )}$ ${\displaystyle (p_{r})=r}$ $${\displaystyle p_{r+1}(x)=(x-a_{r,r})p_{r}(x)-a_{r,r-1}p_{r-1}(x)\cdots -a_{r,0}p_{0}(x)}$$

и скалярное произведение определено $${\displaystyle (f(x),g(x))=\int _{a}^{b}\omega (x)f(x)g(x)dx}$$

для где n — максимальная степень, которая может быть принята за бесконечность, и где . Прежде всего, многочлены, определяемые рекуррентным соотношением, начинающимся с, имеют ведущий коэффициент один и правильную степень. Учитывая начальную точку , ортогональность может быть показана по индукции. Для одного есть ${\displaystyle r=0,1,\ldots ,n-1}$ ${\textstyle a_{r,s}={\frac {\left(xp_{r},p_{s}\right)}{\left(p_{s},p_{s}\right)}}}$ ${\displaystyle p_{0}(x)=1}$ ${\displaystyle p_{0}}$ ${\displaystyle p_{r}}$ ${\displaystyle r=s=0}$ $${\displaystyle (p_{1},p_{0})=(x-a_{0,0})(p_{0},p_{0})=(xp_{0},p_{0})-a_{0,0}(p_{0},p_{0})=(xp_{0},p_{0})-(xp_{0},p_{0})=0.}$$

Теперь, если ортогональны, то также , поскольку во всех скалярных произведениях обращаются в нуль, за исключением первого и того, где встречается тот же ортогональный многочлен. Следовательно, ${\displaystyle p_{0},p_{1},\ldots ,p_{r}}$ ${\displaystyle p_{r+1}}$ $${\displaystyle (p_{r+1},p_{s})=(xp_{r},p_{s})-a_{r,r}(p_{r},p_{s})-a_{r,r-1}(p_{r-1},p_{s})\cdots -a_{r,0}(p_{0},p_{s})}$$ ${\displaystyle p_{s}}$ $${\displaystyle (p_{r+1},p_{s})=(xp_{r},p_{s})-a_{r,s}(p_{s},p_{s})=(xp_{r},p_{s})-(xp_{r},p_{s})=0.}$$

Однако, если скалярное произведение удовлетворяет (что имеет место для квадратуры Гаусса), рекуррентное соотношение сводится к трехчленному рекуррентному соотношению: Для — многочлен степени, меньшей или равной r − 1. С другой стороны, ортогонален каждому многочлену степени, меньшей или равной r − 1. Следовательно, имеем и для s < r − 1. Тогда рекуррентное соотношение упрощается до ${\displaystyle (xf,g)=(f,xg)}$ ${\displaystyle s<r-1,xp_{s}}$ ${\displaystyle p_{r}}$ ${\displaystyle (xp_{r},p_{s})=(p_{r},xp_{s})=0}$ ${\displaystyle a_{r,s}=0}$ $${\displaystyle p_{r+1}(x)=(x-a_{r,r})p_{r}(x)-a_{r,r-1}p_{r-1}(x)}$$

или $${\displaystyle p_{r+1}(x)=(x-a_{r})p_{r}(x)-b_{r}p_{r-1}(x)}$$

(с соглашением ) где ${\displaystyle p_{-1}(x)\equiv 0}$ $${\displaystyle a_{r}:={\frac {(xp_{r},p_{r})}{(p_{r},p_{r})}},\qquad b_{r}:={\frac {(xp_{r},p_{r-1})}{(p_{r-1},p_{r-1})}}={\frac {(p_{r},p_{r})}{(p_{r-1},p_{r-1})}}}$$

(последнее из-за , так как отличается от на степень, меньшую r ). ${\displaystyle (xp_{r},p_{r-1})=(p_{r},xp_{r-1})=(p_{r},p_{r})}$ ${\displaystyle xp_{r-1}}$ ${\displaystyle p_{r}}$

Алгоритм Голуба-Вельша

Трехчленное рекуррентное соотношение можно записать в матричной форме , где , — -й стандартный базисный вектор, т.е. , а J — следующая трехдиагональная матрица , называемая матрицей Якоби: ${\displaystyle J{\tilde {P}}=x{\tilde {P}}-p_{n}(x)\mathbf {e} _{n}}$ ${\displaystyle {\tilde {P}}={\begin{bmatrix}p_{0}(x)&p_{1}(x)&\cdots &p_{n-1}(x)\end{bmatrix}}^{\mathsf {T}}}$ ${\displaystyle \mathbf {e} _{n}}$ ${\displaystyle n}$ ${\displaystyle \mathbf {e} _{n}={\begin{bmatrix}0&\cdots &0&1\end{bmatrix}}^{\mathsf {T}}}$ $${\displaystyle \mathbf {J} ={\begin{bmatrix}a_{0}&1&0&\cdots &0\\b_{1}&a_{1}&1&\ddots &\vdots \\0&b_{2}&\ddots &\ddots &0\\\vdots &\ddots &\ddots &a_{n-2}&1\\0&\cdots &0&b_{n-1}&a_{n-1}\end{bmatrix}}.}$$

Нули полиномов до степени n , которые используются в качестве узлов для квадратуры Гаусса, можно найти, вычислив собственные значения этой матрицы. Эта процедура известна как алгоритм Голуба–Велша . ${\displaystyle x_{j}}$

Для вычисления весов и узлов предпочтительнее рассматривать симметричную трехдиагональную матрицу с элементами ${\displaystyle {\mathcal {J}}}$ $${\displaystyle {\begin{aligned}{\mathcal {J}}_{k,i}=J_{k,i}&=a_{k-1}&k&=1,2,\ldots ,n\\[2.1ex]{\mathcal {J}}_{k-1,i}={\mathcal {J}}_{k,k-1}={\sqrt {J_{k,k-1}J_{k-1,k}}}&={\sqrt {b_{k-1}}}&k&={\hphantom {1,\,}}2,\ldots ,n.\end{aligned}}}$$

То есть,

$${\displaystyle {\mathcal {J}}={\begin{bmatrix}a_{0}&{\sqrt {b_{1}}}&0&\cdots &0\\{\sqrt {b_{1}}}&a_{1}&{\sqrt {b_{2}}}&\ddots &\vdots \\0&{\sqrt {b_{2}}}&\ddots &\ddots &0\\\vdots &\ddots &\ddots &a_{n-2}&{\sqrt {b_{n-1}}}\\0&\cdots &0&{\sqrt {b_{n-1}}}&a_{n-1}\end{bmatrix}}.}$$

J иявляются подобными матрицами и, следовательно, имеют одинаковые собственные значения (узлы). Веса могут быть вычислены из соответствующих собственных векторов: Если— нормализованный собственный вектор (т. е. собственный вектор с евклидовой нормой, равной единице), связанный с собственным значением x _j , соответствующий вес может быть вычислен из первого компонента этого собственного вектора, а именно: ${\displaystyle {\mathcal {J}}}$ ${\displaystyle \phi ^{(j)}}$ $${\displaystyle w_{j}=\mu _{0}\left(\phi _{1}^{(j)}\right)^{2}}$$

где - интеграл весовой функции ${\displaystyle \mu _{0}}$ $${\displaystyle \mu _{0}=\int _{a}^{b}\omega (x)dx.}$$

Более подробную информацию см., например, в (Gil, Segura & Temme 2007).

Оценка погрешности

Погрешность квадратурной формулы Гаусса можно сформулировать следующим образом. ^[5] Для подынтегрального выражения, имеющего 2 n непрерывных производных, для некоторого ξ из ( a , b ) , где p _n — монический (т.е. старший коэффициент равен 1 ) ортогональный многочлен степени n и где $${\displaystyle \int _{a}^{b}\omega (x)\,f(x)\,dx-\sum _{i=1}^{n}w_{i}\,f(x_{i})={\frac {f^{(2n)}(\xi )}{(2n)!}}\,(p_{n},p_{n})}$$ $${\displaystyle (f,g)=\int _{a}^{b}\omega (x)f(x)g(x)\,dx.}$$

В важном частном случае ω ( x ) = 1 мы имеем оценку погрешности ^[6] $${\displaystyle {\frac {\left(b-a\right)^{2n+1}\left(n!\right)^{4}}{(2n+1)\left[\left(2n\right)!\right]^{3}}}f^{(2n)}(\xi ),\qquad a<\xi <b.}$$

Стоер и Булирш отмечают, что эта оценка ошибки неудобна на практике, поскольку может быть трудно оценить производную порядка 2 n , и, кроме того, фактическая ошибка может быть намного меньше границы, установленной производной. Другой подход заключается в использовании двух квадратурных правил Гаусса разных порядков и оценке ошибки как разницы между двумя результатами. Для этой цели могут быть полезны квадратурные правила Гаусса–Кронрода.

Правила Гаусса-Кронрода

Если интервал [ a , b ] подразделяется, то точки оценки Гаусса новых подынтервалов никогда не совпадают с предыдущими точками оценки (за исключением нуля для нечетных чисел), и, таким образом, подынтегральное выражение должно оцениваться в каждой точке. Правила Гаусса–Кронрода являются расширениями квадратурных правил Гаусса, генерируемыми путем добавления n + 1 точек к правилу из n точек таким образом, что результирующее правило имеет порядок 2 n + 1 . Это позволяет вычислять оценки более высокого порядка, повторно используя значения функции оценки более низкого порядка. Разница между квадратурным правилом Гаусса и его расширением Кронрода часто используется в качестве оценки ошибки аппроксимации.

Правила Гаусса-Лобатто

Также известна как квадратура Лобатто , ^[7] названная в честь голландского математика Реуэля Лобатто . Она похожа на квадратуру Гаусса со следующими отличиями:

1. Точки интегрирования включают конечные точки интервала интегрирования.
2. Он точен для полиномов до степени 2n – 3 , где n – число точек интегрирования. ^[8]

Квадратура Лобатто функции f ( x ) на интервале [−1, 1] : $${\displaystyle \int _{-1}^{1}{f(x)\,dx}={\frac {2}{n(n-1)}}[f(1)+f(-1)]+\sum _{i=2}^{n-1}{w_{i}f(x_{i})}+R_{n}.}$$

Абсциссы: x _i - первый ноль , здесь обозначает стандартный полином Лежандра m -й степени, а черточка обозначает производную. ${\displaystyle (i-1)}$ ${\displaystyle P'_{n-1}(x)}$ ${\displaystyle P_{m}(x)}$

Вес: $${\displaystyle w_{i}={\frac {2}{n(n-1)\left[P_{n-1}\left(x_{i}\right)\right]^{2}}},\qquad x_{i}\neq \pm 1.}$$

Остаток: $${\displaystyle R_{n}={\frac {-n\left(n-1\right)^{3}2^{2n-1}\left[\left(n-2\right)!\right]^{4}}{(2n-1)\left[\left(2n-2\right)!\right]^{3}}}f^{(2n-2)}(\xi ),\qquad -1<\xi <1.}$$

Вот некоторые веса:

Адаптивный вариант этого алгоритма с 2 внутренними узлами ^[9] можно найти в GNU Octave и MATLAB как quadlи integrate. ^{[10] [11]}

Ссылки

Цитаты

1. Гаусс 1815
2. Якоби 1826
3. Абрамовиц и Стегун 1983, стр. 887
4. Стер и Булирш 2002, стр. 172–175.
5. Стер и Булирш 2002, Thm 3.6.24
6. Каханер, Молер и Нэш 1989, §5.2
7. Абрамовиц и Стегун 1983, стр. 888
8. Квартерони, Сакко и Салери 2000
9. Гандер и Гаучи 2000
10. MathWorks 2012
11. Итон и др. 2018

Библиография

• Abramowitz, Milton ; Stegun, Irene Ann , ред. (1983) [июнь 1964]. "Глава 25.4, Интеграция". Справочник по математическим функциям с формулами, графиками и математическими таблицами . Серия Applied Mathematics. Том 55 (Девятое переиздание с дополнительными исправлениями десятого оригинального издания с исправлениями (декабрь 1972 г.); первое изд.). Вашингтон, округ Колумбия; Нью-Йорк: Министерство торговли США, Национальное бюро стандартов; Dover Publications. ISBN 978-0-486-61272-0. LCCN  64-60036. MR  0167642. LCCN  65-12253.
• Андерсон, Дональд Г. (1965). "Гауссовы квадратурные формулы для ∫ 0 1 − ln ⁡ ( x ) f ( x ) d x {\displaystyle \int _{0}^{1}-\ln(x)f(x)dx} ". Math. Comp . 19 (91): 477–481. doi : 10.1090/s0025-5718-1965-0178569-1 .
• Данлой, Бернард (1973). «Численное построение квадратурных формул Гаусса для и ». Math. Comp . 27 (124): 861–869. doi :10.1090/S0025-5718-1973-0331730-X. MR  0331730. ${\displaystyle \int _{0}^{1}(-\log x)x^{\alpha }f(x)dx}$ ${\displaystyle \int _{0}^{\infty }E_{m}(x)f(x)dx}$
• Итон, Джон В.; Бейтман, Дэвид; Хауберг, Сорен; Вебринг, Рик (2018). «Функции одной переменной (GNU Octave)» . Проверено 28 сентября 2018 г.
• Гандер, Уолтер; Гаучи, Уолтер (2000). «Адаптивная квадратура — повторный взгляд». BIT Numerical Mathematics . 40 (1): 84–101. doi :10.1023/A:1022318402393.
• Гаусс, Карл Фридрих (1815). Methodus novaintegrium valores для аппроксимации изобретений. Комм. Соц. наук. Геттингенская математика. Том. 3. С. 29–76.datiert 1814, auch в Werke, Band 3, 1876, S. 163–196. Английский перевод Wikisource.
• Гаучи, Вальтер (1968). «Построение квадратурных формул Гаусса–Кристоффеля». Math. Comp . 22 (102): 251–270. doi :10.1090/S0025-5718-1968-0228171-0. MR  0228171.
• Гаучи, Вальтер (1970). «О построении квадратурных правил Гаусса из модифицированных моментов». Math. Comp . 24 (110): 245–260. doi :10.1090/S0025-5718-1970-0285117-6. MR  0285177.
• Гаучи, Вальтер (2020). Репозиторий программного обеспечения для квадратур Гаусса и функций Кристоффеля . SIAM. ISBN 978-1-611976-34-2.
• Гил, Ампаро; Сегура, Хавьер; Темме, Нико М. (2007), «§5.3: квадратура Гаусса», Численные методы для специальных функций , SIAM, ISBN 978-0-89871-634-4
• Голуб, Джин Х.; Уэлш, Джон Х. (1969). «Вычисление квадратурных правил Гаусса». Математика вычислений . 23 (106): 221–230. doi : 10.1090/S0025-5718-69-99647-1 . JSTOR  2004418.
• Якоби, CGJ (1826 г.). «Новый метод Ueber Gauß, die Werthe der Integrale näherungsweise zu finden». Журнал для королевы и математики . 1 . S. 301–308und Werke, Band 6.{{cite journal}}: CS1 maint: postscript (link)
• Кабир, Хоссейн; Матиколаи, Сайед Амир Хоссейн Хассанпур (2017). «Реализация точного обобщенного гауссовского квадратурного решения для поиска упругого поля в однородной анизотропной среде». Журнал Сербского общества вычислительной механики . 11 (1): 11–19. doi :10.24874/jsscm.2017.11.01.02.
• Каханер, Дэвид; Молер, Клив ; Нэш, Стивен (1989). Численные методы и программное обеспечение . Prentice-Hall . ISBN 978-0-13-627258-8.
• Лаудадио, Тереза; Мастронарди, Никола; Ван Доорен, Пол (2023). «Вычисление квадратурных правил Гаусса с высокой относительной точностью». Численные алгоритмы . 92 : 767–793. doi : 10.1007/s11075-022-01297-9 .
• Лори, Дирк П. (1999), «Точное восстановление коэффициентов рекурсии из квадратурных формул Гаусса», J. Comput. Appl. Math. , 112 (1–2): 165–180, doi : 10.1016/S0377-0427(99)00228-9
• Лори, Дирк П. (2001). «Вычисление квадратурных формул типа Гаусса». J. Comput. Appl. Math . 127 (1–2): 201–217. Bibcode : 2001JCoAM.127..201L. doi : 10.1016/S0377-0427(00)00506-9.
• MathWorks (2012). «Численное интегрирование — интеграл MATLAB».
• Piessens, R. (1971). "Гауссовы квадратурные формулы для численного интегрирования интеграла Бромвича и обращения преобразования Лапласа". J. Eng. Math . 5 (1): 1–9. Bibcode :1971JEnMa...5....1P. doi :10.1007/BF01535429.
• Press, WH ; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), «Раздел 4.6. Гауссовы квадратуры и ортогональные многочлены», Numerical Recipes: The Art of Scientific Computing (3-е изд.), Нью-Йорк: Cambridge University Press, ISBN 978-0-521-88068-8
• Квартерони, Альфио ; Сакко, Риккардо; Салери, Фаусто (2000). Численная математика. Нью-Йорк: Springer-Verlag . стр. 425–478. дои : 10.1007/978-3-540-49809-4_10. ISBN 0-387-98959-5.
• Ринер, Кордиан; Швайгхофер, Маркус (2018). «Оптимизационные подходы к квадратуре: новые характеристики гауссовой квадратуры на прямой и квадратуры с несколькими узлами на плоских алгебраических кривых, на плоскости и в более высоких измерениях». Журнал сложности . 45 : 22–54. arXiv : 1607.08404 . doi : 10.1016/j.jco.2017.10.002.
• Сагар, Робин П. (1991). «Квадратура Гаусса для вычисления обобщенных интегралов Ферми-Дирака». Comput. Phys. Commun . 66 (2–3): 271–275. Bibcode :1991CoPhC..66..271S. doi :10.1016/0010-4655(91)90076-W.
• Стоер, Йозеф; Булирш, Роланд (2002), Введение в численный анализ (3-е изд.), Springer , ISBN 978-0-387-95452-3
• Temme, Nico M. (2010), "§3.5(v): Квадратура Гаусса", в Olver, Frank WJ ; Lozier, Daniel M.; Boisvert, Ronald F.; Clark, Charles W. (ред.), Справочник по математическим функциям NIST , Cambridge University Press, ISBN 978-0-521-19225-5, г-н  2723248.
• Якимив, Э. (1996). «Точное вычисление весов в классических квадратурных правилах Гаусса–Кристоффеля». J. Comput. Phys . 129 (2): 406–430. Bibcode : 1996JCoPh.129..406Y. doi : 10.1006/jcph.1996.0258.

Внешние ссылки

• «Квадратурная формула Гаусса», Энциклопедия математики , EMS Press , 2001 [1994]
• ALGLIB содержит коллекцию алгоритмов численного интегрирования (на C# / C++ / Delphi / Visual Basic и т. д.)
• GNU Scientific Library — включает в себя C- версию алгоритмов QUADPACK (см. также GNU Scientific Library )
• От квадратуры Лобатто до постоянной Эйлера e
• Квадратурное правило Гаусса для интегрирования – заметки, PPT, Matlab, Mathematica, Maple, Mathcad в Институте целостных численных методов
• Вайсштейн, Эрик В. «Квадратура Лежандра-Гаусса». Математический мир .
• Гауссова квадратура Криса Маеса и Антона Антонова, Wolfram Demonstrations Project .
• Табличные веса и абсциссы с исходным кодом Mathematica, квадратурные веса и абсциссы Лежандра-Гаусса высокой точности (16 и 256 знаков после запятой) для n =2 – n =64 с исходным кодом Mathematica.
• Исходный код системы Mathematica, распространяемый под лицензией GNU LGPL, предназначен для генерации абсцисс и весов для произвольных весовых функций W(x), областей интегрирования и точности.
• Квадратура Гаусса в Boost.Math для произвольной точности и порядка аппроксимации
• Квадратура Гаусса–Кронрода в Boost.Math
• Узлы и веса квадратуры Гаусса Архивировано 2021-04-14 на Wayback Machine