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

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

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

В численном анализе n - точечное квадратурное правило Гаусса , названное в честь Карла Фридриха Гаусса ^[1], представляет собой квадратурное правило , построенное для получения точного результата для многочленов степени 2 n - 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}}}}$

Можно показать (см. Пресс и др. или Стоер и Булирш), что квадратурные узлы 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 x=\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] . То есть задача состоит в том, чтобы вычислить

$${\displaystyle \int _{a}^{b}\omega (x)\,f(x)\,dx}$$

abωa = −1b = 1ω ( x ) = 1Абрамовица и Стегуна

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

Пусть 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 , чтобы получить

$${\displaystyle h(x)=p_{n}(x)\,q(x)+r(x).}$$

q ( x )n − 1p _nr ( x )n − 1pnn , он должен быть _{ортогонален}q ( 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 , поэтому приведенная выше формула деления говорит нам, что

$${\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 , поскольку все его корни теперь имеют четную кратность. Итак, интеграл

$${\displaystyle \int _{a}^{b}p_{n}(x)\,\left(\prod _{i}(x-x_{i})\right)\,\omega (x)\,dx\neq 0,}$$

ω ( x )pn _{ортогонален}n -1

$${\displaystyle \prod _{i}(x-x_{i})}$$

npnn различных корней, все _{действительные}ab

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

Веса могут быть выражены как

где коэффициент при . Чтобы доказать это, заметим, что, используя интерполяцию Лагранжа , можно выразить r ( x ) как ${\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}}}}$$

r ( x )nnω ( x )ab

$${\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}$$

предоставлено , потому что ${\displaystyle k\leq n}$

$${\displaystyle {\frac {1-\left({\frac {x}{x_{i}}}\right)^{k}}{x-x_{i}}}}$$

k - 1q ( x ) ${\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 , мы имеем ${\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)}$$

s ( x )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)}$

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

Рассмотрим следующий полином степени ${\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}}}}$$

x _jji ${\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 ( n2 ) ^{операций} , метод Ньютона для решения с использованием трехчленной рекуррентности для оценки, требующий O ( n2 ) операций, и асимптотические формулы для больших n , требующие O ( n ) операций . ${\displaystyle p_{n}(x)=0}$

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

Ортогональные полиномы с for для скалярного произведения , степени и старшего коэффициента один (т.е. монические ортогональные полиномы) удовлетворяют рекуррентному соотношению ${\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.}$$

Однако, если скалярное произведение удовлетворяет требованиям (что имеет место для квадратуры Гаусса), рекуррентное соотношение сводится к трехчленному рекуррентному соотношению: For — многочлен степени, меньшей или равной 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 непрерывных производных,

$${\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})}$$

ξ( a , b )pn1 ) _{ортогональный}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. Это верно для полиномов до степени 2 n – 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.}$$

Некоторые из весов:

Адаптивный вариант этого алгоритма с двумя внутренними узлами ^[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 год

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

• Абрамовиц, Милтон ; Стегун, Ирен Энн , ред. (1983) [июнь 1964 г.]. «Глава 25.4, Интеграция». Справочник по математическим функциям с формулами, графиками и математическими таблицами . Серия «Прикладная математика». Том. 55 (Девятое переиздание с дополнительными исправлениями десятого оригинального издания с исправлениями (декабрь 1972 г.); первое изд.). Вашингтон, округ Колумбия; Нью-Йорк: Министерство торговли США, Национальное бюро стандартов; Дуврские публикации. ISBN 978-0-486-61272-0. LCCN  64-60036. МР  0167642. LCCN  65-12253.
• Андерсон, Дональд Г. (1965). «Квадратурные формулы Гаусса для ∫ 0 1 - ln ⁡ ( Икс ) ж ( Икс ) d Икс {\ displaystyle \ int _ {0} ^ {1} - \ ln (x) f (x) dx} ". Математика. Комп . 19 (91): 477–481. дои : 10.1090/s0025-5718-1965-0178569-1 .
• Дэнлой, Бернард (1973). «Численное построение квадратурных формул Гаусса для и ». Математика. Комп . Том. 27, нет. 124. стр. 861–869. doi : 10.1090/S0025-5718-1973-0331730-X. МР  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). «Адаптивная квадратура - еще раз». БИТ Численная математика . 40 (1): 84–101. дои : 10.1023/А: 1022318402393.
• Гаусс, Карл Фридрих (1815). Methodus novaintegrium valores для аппроксимации изобретений. Комм. Соц. наук. Геттингенская математика. Том. 3. С. 29–76.datiert 1814, auch в Werke, Band 3, 1876, S. 163–196. Английский перевод Wikisource.
• Гаучи, Уолтер (1968). «Построение квадратурных формул Гаусса – Кристоффеля». Математика. Комп . Том. 22, нет. 102. С. 251–270. doi : 10.1090/S0025-5718-1968-0228171-0. МР  0228171.
• Гаучи, Уолтер (1970). «О построении правил квадратур Гаусса из модифицированных моментов». Математика. Комп . Том. 24. С. 245–260. дои : 10.1090/S0025-5718-1970-0285117-6. МР  0285177.
• Гаучи, Вальтер (2020). Репозиторий программного обеспечения для гауссовских квадратур и функций Кристоффеля . СИАМ. ISBN 978-1-611976-34-2.
• Гил, Ампаро; Сегура, Хавьер; Темме, Нико М. (2007), «§5.3: квадратура Гаусса», Численные методы для специальных функций , SIAM, ISBN 978-0-89871-634-4
• Голуб, Джин Х .; Уэлш, Джон Х. (1969). «Расчет по правилам квадратур Гаусса». Математика вычислений . 23 (106): 221–230. дои : 10.1090/S0025-5718-69-99647-1 . JSTOR  2004418.
• Якоби, CGJ (1826 г.). «Новый метод Ueber Gauß, die Werthe der Integrale näherungsweise zu finden». Журнал für Reine und Angewandte Mathematik . 1 . S. 301–308und Werke, Band 6.{{cite journal}}: CS1 maint: postscript (link)
• Кабир, Хосейн; Матиколаи, Сайед Амир Хосейн Хасанпур (2017). «Реализация точного обобщенного гауссова квадратурного решения для поиска упругого поля в однородной анизотропной среде» (PDF) . Журнал Сербского общества вычислительной механики . 11 (1): 11–19. doi : 10.24874/jsscm.2017.11.01.02.
• Каханер, Дэвид; Молер, Клив ; Нэш, Стивен (1989). Численные методы и программное обеспечение . Прентис-Холл . ISBN 978-0-13-627258-8.
• Лаудадио, Тереза; Мастронарди, Никола; Ван Доорен, Пол (2023). «Вычисление правил квадратур Гаусса с высокой относительной точностью». Численные алгоритмы . 92 : 767–793.
• Лори, Дирк П. (1999), «Точное восстановление коэффициентов рекурсии из квадратурных формул Гаусса», J. Comput. Прил. Математика. , 112 (1–2): 165–180, doi : 10.1016/S0377-0427(99)00228-9
• Лори, Дирк П. (2001). «Вычисление квадратурных формул типа Гаусса». Дж. Компьютер. Прил. Математика . 127 (1–2): 201–217. Бибкод : 2001JCoAM.127..201L. дои : 10.1016/S0377-0427(00)00506-9.
• Математика (2012). «Численное интегрирование - Интеграл MATLAB».
• Писсенс, Р. (1971). «Квадратурные формулы Гаусса для численного интегрирования интеграла Бромвича и обращения преобразования Лапласа». Дж. Инж. Математика . Том. 5, нет. 1. С. 1–9. Бибкод : 1971JEnMa...5....1P. дои : 10.1007/BF01535429.
• Пресс, WH ; Теукольский, С.А.; Феттерлинг, WT; Фланнери, BP (2007), «Раздел 4.6. Гауссовы квадратуры и ортогональные полиномы», Численные рецепты: искусство научных вычислений (3-е изд.), Нью-Йорк: Cambridge University Press, ISBN 978-0-521-88068-8
• Квартерони, Альфио ; Сакко, Риккардо; Салери, Фаусто (2000). Численная математика. Нью-Йорк: Springer-Verlag . стр. 422, 425. ISBN. 0-387-98959-5.
• Ринер, Кордиан; Швайгофер, Маркус (2018). «Подходы к оптимизации квадратуры: новые характеристики гауссовой квадратуры на прямой и квадратуры с небольшим количеством узлов на плоских алгебраических кривых, на плоскости и в более высоких измерениях». Журнал сложности . 45 : 22–54. arXiv : 1607.08404 . дои : 10.1016/j.jco.2017.10.002.
• Сагар, Робин П. (1991). «Гауссова квадратура для расчета обобщенных интегралов Ферми-Дирака». Вычислить. Физ. Коммун . 66 (2–3): 271–275. Бибкод : 1991CoPhC..66..271S. дои : 10.1016/0010-4655(91)90076-W.
• Стер, Йозеф; Булирш, Роланд (2002), Введение в численный анализ (3-е изд.), Springer , ISBN 978-0-387-95452-3
• Темме, Нико М. (2010), «§3.5(v): Квадратура Гаусса», в Олвере, Фрэнке В.Дж .; Лозье, Дэниел М.; Буасверт, Рональд Ф.; Кларк, Чарльз В. (ред.), Справочник NIST по математическим функциям , издательство Кембриджского университета, ISBN 978-0-521-19225-5, МР  2723248.
• Якимив, Э. (1996). «Точное вычисление весов в классических квадратурных правилах Гаусса – Кристоффеля». Дж. Компьютер. Физ . 129 (2): 406–430. Бибкод : 1996JCoPh.129..406Y. дои : 10.1006/jcph.1996.0258.

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

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