Квадратура Кленшоу – Кертиса

редактировать

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

Кленшоу – Кертис квадрат и квадратурные выражения Фейера - это методы численного интегрирования, или "квадратур", которые основаны на разложении подынтегрального выражения в терминах полиномов Чебышева. Эквивалентно, они используют замену переменных $x = cos ⁡ θ {\ displaystyle x = \ cos \ theta}$ $x = \ cos \ theta$ и используют дискретное косинусное преобразование (DCT) приближение для серии косинусов . Помимо высокой точности сходимости, сравнимой с квадратурными правилами Гаусса, квадратура Кленшоу – Кертиса естественным образом приводит к (где разные порядки точности разделяют точки), что важно как для адаптивной квадратурной, так и для многомерной квадратурной (кубатура ).

Вкратце, функция $f (x) {\ displaystyle f (x)}$ $f (x)$ , подлежащая интеграции, оценивается в $N {\ displaystyle N}$ $N$ экстремумов или корней полинома Чебышева, и эти значения используются для построения полиномиального приближения для функции. Затем этот многочлен точно интегрируется. На практике веса интегрирования для значения функции в каждом узле предварительно вычисляются, и это вычисление может быть выполнено в $O (N log ⁡ N) {\ displaystyle O (N \ log N)}$ $O(N\log N)$ времени с помощью алгоритмов быстрого преобразования Фурье для DCT.

Содержание

1 Общий метод
2 Связь с полиномами Чебышева
3 Квадратура Фейера
4 Сравнение с квадратурой Гаусса
5 Интегрирование с весовыми функциями
6 Интегрирование по бесконечным и полубесконечным интервалам
7 Предварительное вычисление квадратурных весов
8 См. Также
9 Ссылки

Общий метод

Простой способ понять алгоритм - понять, что квадратура Кленшоу – Кертиса (предложенная этими авторами в 1960 г.) сводится к интегрированию посредством изменения переменной x = cos (θ). Алгоритм обычно выражается для интегрирования функции f (x) в интервале [-1,1] (любой другой интервал может быть получен путем соответствующего изменения масштаба). Для этого интеграла можно записать:

∫ - 1 1 f (x) d x = ∫ 0 π f (cos ⁡ θ) sin ⁡ (θ) d θ. {\ Displaystyle \ int _ {- 1} ^ {1} е (х) \, dx = \ int _ {0} ^ {\ pi} е (\ соз \ тета) \ грех (\ тета) \, д \ theta.}

{\ displaystyle \ int _ {- 1} ^ {1} f (x) \, dx = \ int _ {0} ^ {\ pi} f (\ cos \ theta) \ sin (\ theta) \, d \ theta.}

То есть мы преобразовали задачу от интегрирования $f (x) {\ displaystyle f (x)}$ $f (x)$ к одной из задач интегрирования $f (cos ⁡ θ) грех ⁡ θ {\ Displaystyle е (\ соз \ тета) \ грех \ тета}$ $f (\ cos \ theta) \ sin \ theta$ . Это можно сделать, если мы знаем ряд косинусов для $f (cos ⁡ θ) {\ displaystyle f (\ cos \ theta)}$ $f (\ cos \ theta)$ :

f (cos ⁡ θ) = a 0 2 + ∑ К знак равно 1 ∞ ак соз ⁡ (к θ) {\ Displaystyle F (\ соз \ тета) = {\ гидроразрыва {a_ {0}} {2}} + \ sum _ {k = 1} ^ {\ infty } a_ {k} \ cos (k \ theta)}

f (\ cos \ theta) = {\ frac {a_ {0}} {2}} + \ sum _ {k = 1} ^ {\ infty} a_ {k} \ cos (k \ theta)

в этом случае интеграл становится следующим:

∫ 0 π f (cos ⁡ θ) sin ⁡ (θ) d θ = a 0 + ∑ k = 1 ∞ 2 а 2 к 1 - (2 к) 2. {\ Displaystyle \ int _ {0} ^ {\ pi} е (\ соз \ theta) \ sin (\ theta) \, d \ theta = a_ {0} + \ sum _ {k = 1} ^ {\ infty } {\ frac {2a_ {2k}} {1- (2k) ^ {2}}}.}

\int _{0}^{\pi }f(\cos \theta)\sin(\theta)\,d\theta =a_{0}+\sum _{k=1}^{\infty }{\frac {2a_{2k}}{1-(2k)^{2}}}.

Конечно, для вычисления коэффициентов косинусного ряда

ak = 2 π ∫ 0 π f ( соз ⁡ θ) соз ⁡ (к θ) d θ, к знак равно 0, 1, 2,…, {\ displaystyle a_ {k} = {\ frac {2} {\ pi}} \ int _ {0} ^ { \ pi} f (\ cos \ theta) \ cos (k \ theta) \, d \ theta, \ quad k = 0,1,2, \ dots,}

{\ displaystyle a_ {k} = {\ frac {2} {\ pi}} \ int _ {0 } ^ {\ pi} f (\ cos \ theta) \ cos (k \ theta) \, d \ theta, \ quad k = 0,1,2, \ dots,}

нужно снова выполнить числовое интегрирование, поэтому при Во-первых, может показаться, что это не упростило проблему. Однако, в отличие от вычисления произвольных интегралов, интегрирования ряда Фурье для периодических функций (например, $f (cos ⁡ θ) {\ displaystyle f (\ cos \ theta)}$ $f (\ cos \ theta)$ , по построению), вплоть до частоты Найквиста $k = N {\ displaystyle k = N}$ $k = N$ , точно вычисляются с помощью $N + 1 {\ displaystyle N +1}$ $N+1$ равноотстоящие и одинаково взвешенные точки $θ n = n π / N {\ displaystyle \ theta _ {n} = n \ pi / N}$ $\theta _{n}=n\pi /N$ для $n = 0,…, N {\ displaystyle n = 0, \ ldots, N}$ $n = 0, \ ldots, N$ (за исключением того, что конечные точки взвешиваются на 1/2, чтобы избежать двойного счета, что эквивалентно трапециевидной трапеции правило или формула Эйлера – Маклорена ). То есть мы аппроксимируем интеграл косинусной серии дискретным косинусным преобразованием типа I (DCT):

ak ≈ 2 N [f (1) 2 + f (- 1) 2 (- 1) К + ∑ N = 1 N - 1 е (соз ⁡ [N π / N]) соз ⁡ (NK π / N)] {\ Displaystyle а_ {к} \ приблизительно {\ гидроразрыва {2} {N}} \ left [{\ frac {f (1)} {2}} + {\ frac {f (-1)} {2}} (- 1) ^ {k} + \ sum _ {n = 1} ^ { N-1} е (\ соз [п \ пи / N]) \ соз (nk \ пи / N) \ справа]}

a_{k}\approx {\frac {2}{N}}\left[{\frac {f(1)}{2}}+{\frac {f(-1)}{2}}(-1)^{k}+\sum _{n=1}^{N-1}f(\cos[n\pi /N])\cos(nk\pi /N)\right]

для $к = 0,…, N {\ displaystyle k = 0, \ ldots, N}$ $к = 0, \ ldots, N$ , а затем используйте приведенную выше формулу для интеграла в терминах этих $ak {\ displaystyle a_ {k}}$ $a_ {k}$ . Поскольку требуется только $a 2 k {\ displaystyle a_ {2k}}$ $a _ {{2k}}$ , формула дополнительно упрощается до DCT типа I порядка N / 2, предполагая, что N - четное число . :

a 2 k ≈ 2 N [f (1) + f (- 1) 2 + f (0) (- 1) k + ∑ n = 1 N / 2 - 1 {f (cos ⁡ [n π / N]) + е (- соз ⁡ [n π / N])} соз ⁡ (nk π N / 2)] {\ displaystyle a_ {2k} \ приблизительно {\ frac {2} {N}} \ left [{ \ frac {f (1) + f (-1)} {2}} + f (0) (- 1) ^ {k} + \ sum _ {n = 1} ^ {N / 2-1} \ left \ {f (\ cos [n \ pi / N]) + f (- \ cos [n \ pi / N]) \ right \} \ cos \ left ({\ frac {nk \ pi} {N / 2} } \ right) \ right]}

a_ {2k} \ приблизительно {\ frac {2} {N}} \ left [{\ frac {f (1) + f (-1)} {2}} + f (0) (- 1) ^ {k} + \ sum _ {n = 1} ^ {N / 2-1} \ left \ {f (\ cos [n \ pi / N]) + f (- \ cos [n \ pi / N]) \ right \} \ cos \ left ({\ frac {nk \ pi} {N / 2}} \ right) \ right]

Из этой формулы ясно, что квадратурное правило Кленшоу – Кертиса является симметричным в том смысле, что оно взвешивает f (x) и f (−x) одинаково.

Из-за псевдонима вычисляются только коэффициенты $a 2 k {\ displaystyle a_ {2k}}$ $a _ {{2k}}$ до k = N / 2, поскольку дискретная выборка функции делает частоту 2k неотличимой от частоты N – 2k. Эквивалентно, $a 2 k {\ displaystyle a_ {2k}}$ $a _ {{2k}}$ - это амплитуды уникального ограниченного диапазона полинома тригонометрической интерполяции, проходящего через N + 1 точки, где вычисляется f (cos θ), и мы аппроксимируем интеграл интегралом этого интерполяционного полинома. Есть некоторая тонкость в том, как трактовать коэффициент $a N {\ displaystyle a_ {N}}$ $a _ {{N}}$ в интеграле, однако, чтобы избежать двойного счета с его псевдонимом, он включается с весом 1 / 2 в окончательном приближенном интеграле (что также можно увидеть, исследуя интерполирующий полином):

∫ 0 π f (cos ⁡ θ) sin ⁡ (θ) d θ ≈ a 0 + ∑ k = 1 N / 2 - 1 2 а 2 к 1 - (2 к) 2 + а N 1 - N 2. {\ Displaystyle \ int _ {0} ^ {\ pi} е (\ соз \ тета) \ грех (\ тета) \, д \ тета \ приблизительно а_ {0} + \ сумма _ {к = 1} ^ {N / 2-1} {\ frac {2a_ {2k}} {1- (2k) ^ {2}}} + {\ frac {a_ {N}} {1-N ^ {2}}}.}

\ int _ {0} ^ {\ pi} f (\ cos \ theta) \ sin (\ theta) \, d \ theta \ приблизительно a_ {0} + \ sum _ {k = 1} ^ {N / 2-1} {\ frac {2a_ {2k}} {1- (2k) ^ {2}}} + {\ frac {a_ {N}} {1-N ^ {2}}}.

Связь с многочленами Чебышева

Причина, по которой это связано с многочленами Чебышева $T k (x) {\ displaystyle T_ {k} (x)}$ $T_ {k} (x)$ , заключается в том, что по по определению, $T K (соз ⁡ θ) = cos ⁡ (k θ) {\ displaystyle T_ {k} (\ cos \ theta) = \ cos (k \ theta)}$ $T_ {k} (\ cos \ theta) = \ cos (k \ theta)$ , и так приведенный выше ряд косинусов на самом деле является приближением $f (x) {\ displaystyle f (x)}$ $f (x)$ полиномами Чебышева:

f (x) = a 0 2 T 0 (x) + ∑ К знак равно 1 ∞ ак T К (Икс), {\ Displaystyle F (х) = {\ гидроразрыва {а_ {0}} {2}} T_ {0} (х) + \ сумма _ {к = 1} ^ {\ infty} a_ {k} T_ {k} (x),}

f (x) = {\ frac {a_ {0}} {2}} T_ {0} (x) + \ sum _ {k = 1} ^ {\ infty} a_ {k} T_ {k} (x),

и, таким образом, мы «действительно» интегрируем $f (x) {\ displaystyle f (x)}$ $f (x)$ интегрируя его приближенное разложение по многочленам Чебышева. Точки оценки $xn = cos ⁡ (n π / N) {\ displaystyle x_ {n} = \ cos (n \ pi / N)}$ $x_ {n} = \ cos (n \ pi / N)$ соответствуют экстремумам полинома Чебышева $TN (x) {\ displaystyle T_ {N} (x)}$ $T_ {N } (x)$ .

Тот факт, что такое приближение Чебышева представляет собой просто ряд косинусов при замене переменных, отвечает за быстрая сходимость аппроксимации по мере включения большего количества членов $T k (x) {\ displaystyle T_ {k} (x)}$ ${\ displaystyle T_ {k} (x)}$ . Ряд косинусов сходится очень быстро для функций, которые являются даже, периодическими и достаточно гладкими. Здесь это верно, поскольку $f (cos ⁡ θ) {\ displaystyle f (\ cos \ theta)}$ $f (\ cos \ theta)$ четный и периодический в $θ {\ displaystyle \ theta}$ $\theta$ по построению, и дифференцируем в k раз везде, если $f (x) {\ displaystyle f (x)}$ $f (x)$ дифференцируем в k раз на $[- 1, 1] {\ displaystyle [-1,1]}$ $[- 1,1]$ . (Напротив, прямое применение разложения ряда косинусов к $f (x) {\ displaystyle f (x)}$ $f (x)$ вместо $f (cos ⁡ θ) {\ displaystyle f (\ cos \ theta)}$ $f (\ cos \ theta)$ обычно не сходится быстро, потому что наклон четно-периодического расширения обычно будет прерывистым.)

Квадратура Фейера

Фейер предложил два квадратурных правила очень похоже на квадратуру Кленшоу – Кертиса, но намного раньше (в 1933 г.).

Из этих двух «второе» квадратурное правило Фейера почти идентично правилу Кленшоу – Кертиса. Единственное отличие состоит в том, что конечные точки $f (- 1) {\ displaystyle f (-1)}$ $f (-1)$ и $f (1) {\ displaystyle f (1)}$ $f (1)$ установлены на ноль. То есть Фейер использовал только внутренние экстремумы полиномов Чебышева, то есть истинные стационарные точки.

"первое" квадратурное правило Фейера оценивает $ak {\ displaystyle a_ {k}}$ $a_ {k}$ путем вычисления $f (cos ⁡ θ) {\ displaystyle f (\ cos \ theta)}$ $f (\ cos \ theta)$ в другом наборе равноотстоящих точек, на полпути между экстремумами: $θ n = (n + 0,5) π / N {\ displaystyle \ theta _ {n} = (n +0,5) \ pi / N}$ $\ theta _ {n} = (n + 0,5) \ pi / N$ для $0 ≤ n < N {\displaystyle 0\leq n$ $0\leq n<N$ . Это корни $T N (cos ⁡ θ) {\ displaystyle T_ {N} (\ cos \ theta)}$ $T_ {N} (\ cos \ theta)$ и известны как узлы Чебышева. (Эти равноотстоящие средние точки являются единственным другим выбором квадратурных точек, которые сохраняют как четную симметрию косинусного преобразования, так и трансляционную симметрию периодического ряда Фурье.) Это приводит к формуле:

ak ≈ 2 N ∑ N знак равно 0 N - 1 е (соз ⁡ [(n + 0,5) π / N]) соз ⁡ [(n + 0,5) к π / N] {\ displaystyle a_ {k} \ приблизительно {\ frac {2} {N}} \ sum _ {n = 0} ^ {N-1} f (\ cos [(n + 0,5) \ pi /N])\cos[(n+0,5)k\pi / N ]}

a_ {k} \ приблизительно {\ frac {2} {N}} \ sum _ {n = 0} ^ {N-1} f ( \cos[(n+0.5)\pi /N])\cos[(n+0.5)k\pi /N]

что и есть DCT типа II. Однако первое квадратурное правило Фейера не является вложенным: точки оценки для 2N не совпадают ни с одной из точек оценки для N, в отличие от квадратур Кленшоу-Кертиса или второго правила Фейера.

Несмотря на то, что Фейер открыл эти техники до Кленшоу и Кертиса, название «квадратура Кленшоу – Кертиса» стало общепринятым.

Сравнение с квадратурой Гаусса

Классический метод квадратур Гаусса вычисляет подынтегральное выражение при $N + 1 {\ displaystyle N + 1}$ $N+1$ точек и построен для точного интегрирования многочленов до степени $2 N + 1 {\ displaystyle 2N + 1}$ $2N+1$ . Напротив, квадратура Кленшоу – Кертиса, приведенная выше, вычисляет подынтегральное выражение в $N + 1 {\ displaystyle N + 1}$ $N+1$ точек и точно интегрирует многочлены только до степени $N {\ displaystyle N }$ $N$ . Таким образом, может показаться, что Кленшоу – Кертис по своей природе хуже, чем гауссова квадратура, но на самом деле это не так.

На практике несколько авторов заметили, что точность Кленшоу – Кертиса может быть сопоставима с точностью квадратуры Гаусса для того же количества точек. Это возможно, потому что большинство числовых подынтегральных выражений не являются полиномами (особенно потому, что полиномы можно интегрировать аналитически), а приближение многих функций с помощью полиномов Чебышева быстро сходится (см. приближение Чебышева ). Фактически, недавние теоретические результаты утверждают, что и квадратура Гаусса, и квадратура Кленшоу – Кертиса имеют ошибку, ограниченную $O ([2 N] - k / k) {\ displaystyle O ([2N] ^ {- k} / k)}$ $O([2N]^{{-k}}/k)$ для k-кратно дифференцируемого подынтегрального выражения.

Одним из часто упоминаемых преимуществ квадратур Кленшоу – Кертиса является то, что квадратурные веса могут быть вычислены в $O (N log ⁡ N) {\ displaystyle O (N \ log N)}$ $O(N\log N)$ времени с помощью алгоритмов быстрого преобразования Фурье (или их аналогов для DCT), тогда как для большинства алгоритмов квадратурных весов Гаусса требуется $O (N 2) {\ displaystyle O (N ^ {2})}$ $O(N^{2})$ время вычислить. Однако последние алгоритмы достигли $O (N) {\ displaystyle O (N)}$ $O (N)$ сложности для квадратур Гаусса – Лежандра. На практике численное интегрирование высокого порядка редко выполняется простым вычислением квадратурной формулы для очень большого $N {\ displaystyle N}$ $N$ . Вместо этого обычно используется схема адаптивной квадратурной, которая сначала оценивает интеграл до низкого порядка, а затем последовательно улучшает точность, увеличивая количество точек выборки, возможно, только в областях, где интеграл неточен. Чтобы оценить точность квадратур, сравнивают ответ с ответом квадратурного правила даже более низкого порядка. В идеале это правило квадратуры низшего порядка оценивает подынтегральное выражение в подмножестве исходных N точек, чтобы минимизировать оценки подынтегрального выражения. Это называется a, и здесь Кленшоу – Кертис имеет то преимущество, что правило для порядка N использует подмножество точек из порядка 2N. Напротив, квадратурные правила Гаусса не являются естественно вложенными, и поэтому необходимо использовать квадратурные формулы Гаусса – Кронрода или аналогичные методы. Вложенные правила также важны для разреженных сеток в многомерной квадратуре, и квадратура Кленшоу – Кертиса является популярным методом в этом контексте.

Интеграция с весовыми функциями

В более общем смысле, можно поставить задачу интегрирования произвольного $f (x) {\ displaystyle f (x)}$ $f (x)$ против функции фиксированного веса $w (x) {\ displaystyle w (x)}$ $w (x)$ , который известен заранее:

∫ - 1 1 f (x) w (x) dx = ∫ 0 π f (cos ⁡ θ) w (cos ⁡ θ) sin ⁡ (θ) d θ. {\ Displaystyle \ int _ {- 1} ^ {1} е (х) вес (х) \, dx = \ int _ {0} ^ {\ pi} е (\ соз \ тета) вес (\ соз \ тета) \ sin (\ theta) \, d \ theta.}

{\ Displaystyle \ int _ {- 1} ^ {1} е (х) вес (х) \, dx = \ int _ {0} ^ {\ pi} е (\ соз \ тета) ш (\ соз \ тета) \ грех (\ тета) \, д \ тета.}

Наиболее частый случай: $w (x) = 1 {\ displaystyle w (x) = 1}$ $w(x)=1$ , как указано выше, но в некоторых приложениях желательна другая весовая функция. Основная причина заключается в том, что, поскольку $w (x) {\ displaystyle w (x)}$ $w (x)$ может приниматься во внимание априори, ошибка интегрирования может зависеть только от точности аппроксимации $f (x) {\ displaystyle f (x)}$ $f (x)$ , независимо от того, насколько плохо может себя вести весовая функция.

Квадратуру Кленшоу – Кертиса можно обобщить на этот случай следующим образом. Как и прежде, он работает, находя разложение в ряд косинусов $f (cos ⁡ θ) {\ displaystyle f (\ cos \ theta)}$ $f (\ cos \ theta)$ через DCT, а затем интегрируя каждый член в косинусный ряд. Однако теперь эти интегралы имеют вид

W k = ∫ 0 π w (cos ⁡ θ) cos ⁡ (k θ) sin ⁡ (θ) d θ. {\ displaystyle W_ {k} = \ int _ {0} ^ {\ pi} w (\ cos \ theta) \ cos (k \ theta) \ sin (\ theta) \, d \ theta.}

W_ {k} = \ int _ {0} ^ {\ pi} w (\ cos \ theta) \ cos (k \ theta) \ sin (\ theta) \, d \ theta.

для most $w (x) {\ displaystyle w (x)}$ $w (x)$ , этот интеграл нельзя вычислить аналитически, в отличие от предыдущего. Поскольку одна и та же весовая функция обычно используется для многих подынтегральных выражений $f (x) {\ displaystyle f (x)}$ $f (x)$ , однако, можно позволить вычислить эти $W k {\ displaystyle W_ {k}}$ $W_ {k}$ заранее численно с высокой точностью. Более того, поскольку $w (x) {\ displaystyle w (x)}$ $w (x)$ обычно задается аналитически, иногда можно использовать специальные методы для вычисления $W k {\ displaystyle W_ {k}}$ $W_ {k}$ .

Например, были разработаны специальные методы для применения квадратуры Кленшоу – Кертиса к подынтегральным выражениям вида $f (x) w (x) {\ displaystyle f (x) w (x)}$ $f (x) w (x)$ с весовой функцией $w (x) {\ displaystyle w (x)}$ $w (x)$ , которая сильно колеблется, например синусоидой или функцией Бесселя (см., например, Evans Webster, 1999). Это полезно для высокоточных вычислений рядов Фурье и рядов Фурье – Бесселя, где простой $w (x) = 1 {\ displaystyle w (x) = 1} <Квадратурные методы 282>$ $w(x)=1$ проблематичны из-за высокой точности, необходимой для определения вклада быстрых колебаний. Здесь быстроосциллирующая часть подынтегрального выражения учитывается с помощью специальных методов для $W k {\ displaystyle W_ {k}}$ $W_ {k}$ , тогда как неизвестная функция $f (x) { \ displaystyle f (x)}$ $f (x)$ обычно ведет себя лучше.

Другой случай, когда весовые функции особенно полезны, - это если подынтегральное выражение неизвестно, но имеет известную особенность некоторой формы, например известный разрыв или интегрируемая расходимость (например, 1 / √x) в некоторой точке. В этом случае сингулярность может быть преобразована в весовую функцию $w (x) {\ displaystyle w (x)}$ $w (x)$ , и ее аналитические свойства могут использоваться для вычисления $W k {\ displaystyle W_ {k}}$ $W_ {k}$ точно заранее.

Обратите внимание, что квадратура Гаусса также может быть адаптирована для различных весовых функций, но методика несколько отличается. В квадратуре Кленшоу – Кертиса подынтегральное выражение всегда вычисляется в одном и том же наборе точек независимо от $w (x) {\ displaystyle w (x)}$ $w (x)$ , соответствующего экстремумам или корням чебышевского полином. В квадратуре Гаусса разные весовые функции приводят к разным ортогональным многочленам и, следовательно, к разным корням, по которым вычисляется подынтегральное выражение.

Интегрирование на бесконечных и полубесконечных интервалах

Также можно использовать квадратуру Кленшоу – Кертиса для вычисления интегралов вида $∫ 0 ∞ f (x) dx {\ displaystyle \ int _ {0} ^ {\ infty} f (x) \, dx}$ ${\ displaystyle \ int _ {0} ^ {\ infty} f (x) \, dx}$ и $∫ - ∞ ∞ f (x) dx {\ displaystyle \ int _ {- \ infty} ^ {\ infty} f (x) \, dx}$ $\ int _ {- \ infty} ^ {\ infty} f (x) \, dx$ , используя технику переназначения координат. Высокая точность, даже экспоненциальная сходимость для гладких подынтегральных выражений, может сохраняться до тех пор, пока $f (x) {\ displaystyle f (x)}$ $f (x)$ затухает достаточно быстро, как | x | приближается к бесконечности.

Одна из возможностей - использовать типичное преобразование координат, такое как x = t / (1 − t)

∫ - ∞ ∞ f (x) dx = ∫ - 1 1 f (t 1 - t 2) 1 + T 2 (1 - T 2) 2 dt, {\ displaystyle \ int _ {- \ infty} ^ {\ infty} f (x) \, dx = \ int _ {- 1} ^ {1} f \ left ({\ frac {t} {1-t ^ {2}}} \ right) {\ frac {1 + t ^ {2}} {(1-t ^ {2}) ^ {2}}} \, dt,}

\int _{-\infty }^{\infty }f(x)\,dx=\int _{-1}^{1}f\left({\frac {t}{1-t^{2}}}\right){\frac {1+t^{2}}{(1-t^{2})^{2}}}\,dt,

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

Например, можно использовать переназначение координат $x = L cot 2 ⁡ (θ / 2) {\ displaystyle x = L \ cot ^ {2} (\ theta / 2)}$ $x = L \ cot ^ {2} (\ theta / 2)$ , где L - заданная пользователем константа (можно просто использовать L = 1; оптимальный выбор L может ускорить сходимость, но зависит от проблемы), чтобы преобразовать полубесконечный интеграл в:

∫ 0 ∞ f (x) dx = 2 L ∫ 0 π f [L cot 2 ⁡ (θ / 2)] [1 - cos ⁡ (θ)] 2 sin ⁡ (θ) d θ. {\ displaystyle \ int _ {0} ^ {\ infty} f (x) \, dx = 2L \ int _ {0} ^ {\ pi} {\ frac {f [L \ cot ^ {2} (\ theta / 2)]} {[1- \ cos (\ theta)] ^ {2}}} \ sin (\ theta) \, d \ theta.}

\int _{0}^{\infty } f(x)\,dx=2L\int _{0}^{\pi }{\frac {f[L\cot ^{2}(\theta /2)]}{[1-\cos(\theta)]^{2}}}\sin(\theta)\,d\theta.

Коэффициент, умножающий sin (θ), f (...) / (...), затем можно разложить в косинусный ряд (приблизительно, используя дискретное косинусное преобразование) и интегрировать почленно, точно так же, как это было сделано для f (cos θ) выше. Чтобы устранить сингулярность при θ = 0 в этом подынтегральном выражении, достаточно просто, чтобы f (x) стремилась к нулю достаточно быстро, когда x приближается к бесконечности, и, в частности, f (x) должна убывать по крайней мере так быстро, как 1 / x.>

Для дважды бесконечного интервала интегрирования можно использовать перераспределение координат $x = L cot ⁡ (θ) {\ displaystyle x = L \ cot (\ theta)}$ $x = L \ cot (\ theta)$ (где L - заданная пользователем константа, как указано выше), чтобы преобразовать интеграл в:

∫ - ∞ ∞ f (x) dx = L ∫ 0 π f [L cot ⁡ (θ)] sin 2 ⁡ (θ) d θ ≈ L π N ∑ n = 1 N - 1 f [L cot ⁡ (n π / N)] sin 2 ⁡ (n π / N). {\ displaystyle \ int _ {- \ infty} ^ {\ infty} f (x) \, dx = L \ int _ {0} ^ {\ pi} {\ frac {f [L \ cot (\ theta)] } {\ sin ^ {2} (\ theta)}} \, d \ theta \ приблизительно {\ frac {L \ pi} {N}} \ sum _ {n = 1} ^ {N-1} {\ frac {f [L \ cot (n \ pi / N)]} {\ sin ^ {2} (n \ pi / N)}}.}

{\ displaystyle \ int _ {- \ infty} ^ {\ infty} f (x) \, dx = L \ int _ {0} ^ {\ pi} {\ frac {f [L \ cot (\ theta)]} {\ sin ^ {2} (\ theta)}} \, d \ theta \ приблизительно { \ frac {L \ pi} {N}} \ sum _ {n = 1} ^ {N-1} {\ frac {f [L \ cot (n \ pi / N)]} {\ sin ^ {2} (n \ pi /N)}}.}

В этом случае мы использовали тот факт, что подынтегральное выражение f (L cotθ) / sin (θ) уже является периодическим, поэтому его можно напрямую интегрировать с высокой (даже экспоненциальной) точностью, используя правило трапеций (при условии, что f достаточно гладкая и быстро затухает); нет необходимости вычислять ряд косинусов в качестве промежуточного шага. Обратите внимание, что квадратурное правило не включает конечные точки, где мы предположили, что подынтегральная функция стремится к нулю. Приведенная выше формула требует, чтобы f (x) затухал быстрее, чем 1 / x, когда x стремится к ± ∞. (Если f убывает точно как 1 / x, тогда подынтегральное выражение переходит к конечному значению в конечных точках, и эти пределы должны быть включены в качестве конечных членов в правило трапеций.) Однако, если f убывает только полиномиально быстро, то может оказаться необходимым использовать следующий шаг квадратур Кленшоу – Кертиса для получения экспоненциальной точности переназначенного интеграла вместо правила трапеций, в зависимости от более подробной информации об ограничивающих свойствах f: Проблема в том, что, хотя f (L cotθ) / sin (θ) действительно периодичен с периодом π, он не обязательно будет гладким на концах, если все производные там не обращаются в нуль [например, функция f (x) = tanh (x) / x убывает как 1 / x, но имеет скачкообразный разрыв в наклоне преобразованной функции при θ = 0 и π].

Другой подход с перераспределением координат был предложен для интегралов вида $∫ 0 ∞ e - xg (x) dx {\ displaystyle \ int _ {0} ^ {\ infty} e ^ {- x } g (x) \, dx}$ ${\ displaystyle \ int _ {0} ^ {\ infty} e ^ {- x} g (x) \, dx}$ , и в этом случае можно использовать преобразование $x = - ln ⁡ [(1 + cos ⁡ θ) / 2] {\ displaystyle x = - \ ln [(1+ \ cos \ theta) / 2]}$ $x=-\ln[(1+\cos \theta)/2]$ для преобразования интеграла в форму $∫ 0 π f (cos cos θ) sin ⁡ θ d θ {\ displaystyle \ int _ {0} ^ {\ pi} f (\ cos \ theta) \ sin \ theta \, d \ theta}$ $\ int _ {0} ^ {\ pi} f (\ cos \ theta) \ sin \ theta \, d \ theta$ где $f (u) = g (- ln ⁡ [(1 + u) / 2]) / 2 {\ displaystyle f (u) = g (- \ ln [(1 + u) / 2]) / 2}$ $f (u) = g (- \ ln [(1 + u) / 2]) / 2$ , после чего можно перейти точно к Кленшоу– Квадратура Кертиса для f, как указано выше. Однако из-за особенностей концевых точек в этом переотображении координат используется первое квадратурное правило Фейера [которое не вычисляет f (−1)], если g (∞) не является конечным.

Предварительное вычисление квадратурных весов

На практике неудобно выполнять DCT значений выборочной функции f (cosθ) для каждого нового подынтегрального выражения. Вместо этого обычно предварительно вычисляются квадратурные веса $wn {\ displaystyle w_ {n}}$ $w _ {n}$ (для n от 0 до N / 2, предполагая, что N четно), так что

∫ - 1 1 f (x) dx ≈ ∑ n = 0 N / 2 wn {f (cos ⁡ [n π / N]) + f (- cos ⁡ [n π / N])}. {\ displaystyle \ int _ {- 1} ^ {1} е (х) \, dx \ приблизительно \ сумма _ {n = 0} ^ {N / 2} w_ {n} \ left \ {f (\ cos [ n \ pi / N]) + f (- \ cos [n \ pi / N]) \ right \}.}

\int _{-1}^{1}f(x)\,dx\approx \sum _{n=0}^{N/2}w_{n}\left\{f(\cos [n\pi /N])+f(- \cos[n\pi /N])\right\}.

Эти веса $wn {\ displaystyle w_ {n}}$ $w _ {n}$ также вычисляются с помощью DCT, что легко увидеть, выразив вычисление в терминах матрицы алгебры. В частности, мы вычислили коэффициенты ряда косинусов $a 2 k {\ displaystyle a_ {2k}}$ $a _ {{2k}}$ с помощью выражения вида:

c = (a 0 a 2 a 4 ⋮ a N) знак равно D (Y 0 Y 1 Y 2 ⋮ Y N / 2) = D Y, {\ displaystyle c = {\ begin {pmatrix} a_ {0} \\ a_ {2} \\ a_ {4} \\ \ vdots \\ a_ {N} \ end {pmatrix}} = D {\ begin {pmatrix} y_ {0} \\ y_ {1} \\ y_ {2} \\\ vdots \\ y_ {N / 2} \ end {pmatrix}} = Dy,}

c = {\ begin {pmatrix} a_ {0} \\ a_ {2} \\ a_ {4} \\\ vdots \\ a_ {N} \ end {pmatrix}} = D {\ begin {pmatrix} y_ {0} \\ y_ {1} \\ y_ {2} \\\ vdots \\ y_ {N / 2} \ end {pmatrix}} = Dy,

где D - матричная форма (N / 2 + 1) -точки типа I DCT сверху, с записями (для индексы с отсчетом от нуля):

D kn = 2 N cos ⁡ (nk π N / 2) × {1/2 n = 0, N / 2 1 в противном случае {\ displaystyle D_ {kn} = {\ frac {2} {N}} \ cos \ left ({\ frac {nk \ pi} {N / 2}} \ right) \ times {\ begin {cases} 1/2 n = 0, N / 2 \\ 1 \ mathrm {иначе} \ end {case}}}

D_ {kn} = {\ frac {2} {N}} \ cos \ left ({\ frac {nk \ pi} { N / 2}} \ right) \ times {\ begin {cases} 1/2 n = 0, N / 2 \\ 1 \ mathrm {в противном случае} \ end {ases}}

и $yn {\ displaystyle y_ {n}}$ $y_ {n}$ is

yn = f (cos ⁡ [n π / N]) + f (- cos ⁡ [n π / N]). {\ displaystyle y_ {n} = f (\ cos [n \ pi / N]) + f (- \ cos [n \ pi / N]). \!}

y_ {n} = f (\ cos [n \ pi / N]) + е (- \ соз [п \ пи /N]).\!

Как обсуждалось выше, из-за псевдонима, нет смысла вычислять коэффициенты за пределами $a N {\ displaystyle a_ {N}}$ $a_{N}$ , поэтому D является $(N / 2 + 1) × (N / 2 + 1) {\ displaystyle (N / 2 + 1) \ times (N / 2 + 1)}$ $(N/2+1)\times (N/2+1)$ матрица. В терминах этих коэффициентов c интеграл приблизительно равен:

∫ - 1 1 f (x) dx ≈ a 0 + ∑ k = 1 N / 2 - 1 2 a 2 k 1 - (2 k) 2 + a N 1 - N 2 знак равно d T c, {\ displaystyle \ int _ {- 1} ^ {1} f (x) \, dx \ приблизительно a_ {0} + \ sum _ {k = 1} ^ {N / 2-1} {\ frac {2a_ {2k}} {1- (2k) ^ {2}}} + {\ frac {a_ {N}} {1-N ^ {2}}} = d ^ {T } c,}

\int _{-1}^{1}f(x)\,dx\approx a_{0}+\sum _{k=1}^{N/2-1}{\frac {2a_{2k}}{1-(2k)^{2}}}+{\frac {a_{N}}{1-N^{2}}}=d^{T}c,

сверху, где c - вектор коэффициентов $a 2 k {\ displaystyle a_ {2k}}$ $a _ {{2k}}$ выше, а d - вектор интегралов для каждого коэффициента Фурье:

d = (1 2 / (1 - 4) 2 / (1 - 16) ⋮ 2 / (1 - [N - 2] 2) 1 / (1 - N 2)). {\ Displaystyle d = {\ begin {pmatrix} 1 \\ 2 / (1-4) \\ 2 / (1-16) \\\ vdots \\ 2 / (1- [N-2] ^ {2}) \\ 1 / (1-N ^ {2}) \ end {pmatrix}}.}

d={\begin{pmatrix}1\\2/(1-4)\\2/( 1-16)\\\vdots \\2/(1-[N-2]^{2})\\1/(1-N^{2})\end{pmatrix}}.

(Обратите внимание, однако, что эти весовые коэффициенты изменяются, если изменить DCT-матрицу D для использования другого соглашения о нормализации. Например, обычно DCT типа I определяют с дополнительными множителями 2 или √2 множителей в первой и последней строках или столбцах, что приводит к соответствующим изменениям в записях d.) $d T c { \ displaystyle d ^ {T} c}$ $d ^ {T} c$ суммирование можно преобразовать в:

∫ - 1 1 f (x) dx ≈ d T c = d TD y = (DT d) T y знак равно вес T Y {\ Displaystyle \ int _ {- 1} ^ {1} f (x) \, dx \ приблизительно d ^ {T} c = d ^ {T} Dy = (D ^ {T} d) ^ {T} y = w ^ {T} y}

{\ displaystyle \ int _ {- 1} ^ {1} f (x) \, dx \ приблизительно d ^ {T} c = d ^ {T} Dy = (D ^ {T} d) ^ {T} y = w ^ { T} y}

где w - вектор желаемых весов $wn {\ displaystyle w_ {n}}$ $w _ {n}$ выше, с:

w = DT d. {\ displaystyle w = D ^ {T} d. \!}

w = D ^ {T} d. \!

Поскольку транспонированная матрица $DT {\ displaystyle D ^ {T}}$ $D ^ {T}$ также является DCT (например, транспонирование DCT типа I является DCT типа I, возможно, с немного другой нормализацией в зависимости от используемых соглашений), квадратурные веса w могут быть предварительно вычислены за время O (N log N) для заданное N с использованием быстрых алгоритмов DCT.

Веса $wn {\ displaystyle w_ {n}}$ $w _ {n}$ положительны и их сумма равна единице.

См. Также

Литература

^W. Морвен Джентльмен, "Реализация квадратуры Кленшоу-Кертиса I: Методология и опыт", Коммуникации ACM 15 (5), p. 337-342 (1972).
^Йорг Вальдфогель, "Быстрое построение квадратурных правил Фейера и Кленшоу-Кертиса ", BIT Numerical Mathematics 46 (1), p. 195-202 (2006).
^С. W. Clenshaw and AR Curtis "Метод численного интегрирования на автоматическом компьютере Numerische Mathematik 2, 197 (1960).
^JP Boyd, Chebychev and Fourier Spectral Methods, 2nd ред. (Довер, Нью-Йорк, 2001).
^См., например, С.Г. Джонсон, «Заметки о конвергенции квадратуры трапецеидального правила », онлайн-курсы MIT (2008).
^Леопольд Фейер, "О бесконечных последовательностях, возникающих в теориях гармонического анализа, интерполяции и механических квадратур ", Бюллетень Американского математического общества 39 (1933), стр. 521–534. Леопольд Фейер, "Mechanische Quadraturen mit Positiven Cotesschen Zahlen, Mathematische Zeitschrift 37, 287 (1933)".
^Трефетен, Ллойд Н. (2008). «Квадратура Гаусса лучше, чем Кленшоу-Кертис?». SIAM Обзор. 50 (1): 67–87. CiteSeerX 10.1.1.157.4174. doi : 10.1137 / 060659831.
^Игнас Богерт, Безоперационное вычисление квадратурных узлов и весов Гаусса - Лежандра, SIAM Journal on Scientific Computing vol. 36, стр. A1008 – A1026 (2014)
^Эрих Новак и Клаус Риттер, «Интегрирование гладких функций с высокой размерностью над кубами», Numerische Mathematik vol. 75, стр. 79–97 (1996).
^Г. А. Эванс и Дж. Р. Вебстер, "Сравнение некоторых методов для вычисления сильно колеблющихся интегралов", Журнал вычислительной и прикладной математики, вып. 112, с. 55-69 (1999).
^ Джон П. Бойд, "Экспоненциально сходящиеся квадратурные схемы Фурье – Чебшева [sic] на ограниченных и бесконечных интервалах", J. Scientific Computing 2 (2), с. 99-109 (1987).
^Нирмал Кумар Басу и Мадхав Чандра Кунду, «Некоторые методы численного интегрирования на полубесконечном интервале», Приложения математики 22 (4), с. 237-243 (1977).
^Дж. П. Имхоф, "О методе численного интегрирования Кленшоу и Кертиса", Numerische Mathematik 5, p. 138-141 (1963).