Квадратура Кленшоу – Кертиса - Clenshaw–Curtis quadrature

Clenshaw-Curtis квадратурные и Фейеровские квадратурные являются методами численного интегрирования , или «квадратуры», которые основаны на разложении подынтегрального в терминах полиномов Чебышева . Эквивалентно, они используют замену переменных и приближение дискретного косинусного преобразования (DCT) для ряда косинусов . Помимо наличия быстрой сходящейся точности, сравнимой с квадратурными правилами Гаусса, квадратурная формула Кленшоу – Кертиса естественным образом приводит к вложенным квадратурным правилам (где разные порядки точности разделяют точки), что важно как для адаптивной квадратурной, так и для многомерной квадратурной ( кубатуры ).

Вкратце, интегрируемая функция оценивается в экстремумах или корнях полинома Чебышева, и эти значения используются для построения полиномиального приближения для функции. Затем этот многочлен точно интегрируется. На практике веса интегрирования для значения функции в каждом узле предварительно вычисляются, и это вычисление может быть выполнено во времени с помощью алгоритмов быстрого преобразования Фурье для DCT.

Общий метод

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

То есть мы превратили задачу из интеграции в задачу интеграции . Это можно сделать, если мы знаем ряд косинусов для :

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

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

нужно снова выполнить числовое интегрирование, так что сначала может показаться, что это не упростило задачу. Однако, в отличие от вычисления произвольных интегралов, интегрирования ряда Фурье для периодических функций (например , по построению), вплоть до частоты Найквиста , точно вычисляются по равноотстоящим и одинаково взвешенным точкам для (за исключением того, что конечные точки взвешиваются на 1/2 , чтобы избежать двойного счета, эквивалентного правилу трапеций или формуле Эйлера – Маклорена ). То есть, мы аппроксимируем интеграл ряда косинусов с помощью дискретного косинусного преобразования типа I (DCT):

для, а затем используйте приведенную выше формулу для интеграла через них . Поскольку требуется только , формула упрощается до DCT типа I порядка N / 2, предполагая, что N - четное число :

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

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

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

Причина, по которой это связано с многочленами Чебышева, заключается в том, что по определению , и поэтому косинусный ряд, приведенный выше, на самом деле является приближением многочленами Чебышева:

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

Тот факт, что такое приближение Чебышева представляет собой просто ряд косинусов при замене переменных, отвечает за быструю сходимость приближения по мере включения большего количества членов . Ряд косинусов сходится очень быстро для функций, которые являются четными , периодическими и достаточно гладкими. Это верно здесь, так как даже и периодические по по построению, и к шрифт Times дифференцируема всюду , если есть к шрифт Times дифференцируема на . (Напротив, прямое применение разложения в ряд косинусов вместо , как правило, не будет быстро сходиться, потому что наклон четно-периодического расширения обычно будет прерывистым.)

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

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

Из этих двух «второе» квадратурное правило Фейера почти идентично правилу Кленшоу – Кертиса. Единственная разница в том, что конечные точки и установлены в ноль. То есть Фейер использовал только внутренние экстремумы полиномов Чебышева, то есть истинные стационарные точки.

«Первое» квадратурное правило Фейера оценивает , оценивая в другом наборе равноотстоящих точек, на полпути между экстремумами: для . Таковы корни из , и известны как узлы Чебышева . (Эти равноотстоящие средние точки являются единственным другим выбором квадратурных точек, которые сохраняют как четную симметрию косинусного преобразования, так и трансляционную симметрию периодического ряда Фурье.) Это приводит к формуле:

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

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

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

Классический метод квадратуры Гаусса вычисляет подынтегральное выражение в точках и построен для точного интегрирования многочленов с точностью до степени . Напротив, квадратура Кленшоу – Кертиса, приведенная выше, вычисляет подынтегральное выражение в точках и точно интегрирует многочлены только до степени . Таким образом, может показаться, что Кленшоу – Кертис по своей природе хуже, чем гауссовская квадратура, но на самом деле это не так.

На практике несколько авторов заметили, что точность Кленшоу – Кертиса может быть сопоставима с точностью квадратуры Гаусса для того же количества точек. Это возможно, потому что большинство числовых подынтегральных выражений не являются многочленами (особенно потому, что многочлены могут быть проинтегрированы аналитически), а приближение многих функций с помощью многочленов Чебышева быстро сходится (см. Приближение Чебышева ). Фактически, недавние теоретические результаты доказывают, что и квадратура Гаусса, и квадратура Кленшоу – Кертиса имеют ошибку, ограниченную для k- кратного дифференцируемого подынтегрального выражения.

Одним из часто упоминаемых преимуществ квадратур Кленшоу – Кертиса является то, что квадратурные веса могут быть вычислены во времени с помощью алгоритмов быстрого преобразования Фурье (или их аналогов для DCT), тогда как большинство алгоритмов для квадратурных весов Гаусса требовали времени для вычисления. Однако недавние алгоритмы достигли сложности для квадратур Гаусса – Лежандра. На практике числовое интегрирование высокого порядка редко выполняется простым вычислением квадратурной формулы для очень больших значений . Вместо этого обычно используется адаптивная квадратурная схема, которая сначала оценивает интеграл до низкого порядка, а затем последовательно повышает точность, увеличивая количество точек выборки, возможно, только в тех областях, где интеграл неточен. Чтобы оценить точность квадратуры, нужно сравнить ответ с ответом квадратурного правила еще более низкого порядка. В идеале это правило квадратуры низшего порядка оценивает подынтегральное выражение в подмножестве исходных N точек, чтобы минимизировать оценки подынтегрального выражения. Это называется вложенное правило квадратурного , и здесь Clenshaw-Curtis имеет то преимущество , что правило для порядка N использует подмножество точек из порядка 2 N . Напротив, квадратурные правила Гаусса не являются естественными вложенными, поэтому необходимо использовать квадратурные формулы Гаусса – Кронрода или аналогичные методы. Вложенные правила также важны для разреженных сеток в многомерной квадратуре, и квадратура Кленшоу – Кертиса является популярным методом в этом контексте.

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

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

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

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

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

Например, были разработаны специальные методы для применения квадратуры Кленшоу – Кертиса к подынтегральным выражениям формы с сильно колеблющейся весовой функцией , например синусоидой или функцией Бесселя (см., Например, Evans & Webster, 1999). Это полезно для высокоточных вычислений рядов Фурье и Фурье – Бесселя , где простые квадратурные методы проблематичны из-за высокой точности, необходимой для определения вклада быстрых колебаний. Здесь быстроосциллирующая часть подынтегрального выражения учитывается с помощью специализированных методов для , тогда как неизвестная функция обычно ведет себя лучше.

Другой случай, когда весовые функции особенно полезны, - это если подынтегральное выражение неизвестно, но имеет известную особенность некоторой формы, например, известный разрыв или интегрируемую расходимость (например, 1 / x ) в некоторой точке. В этом случае особенность может быть включена в весовую функцию, и ее аналитические свойства могут быть использованы для точного вычисления заранее.

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

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

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

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

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

Например, можно использовать перераспределение координат , где L - заданная пользователем константа (можно просто использовать L = 1; оптимальный выбор L может ускорить сходимость, но зависит от проблемы), чтобы преобразовать полубесконечный интеграл в:

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

Для вдвойне бесконечного интервала интегрирования можно использовать перераспределение координат (где L - заданная пользователем константа, как указано выше), чтобы преобразовать интеграл в:

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

Другой подход перераспределения координат был предложен для интегралов формы , и в этом случае можно использовать преобразование, чтобы преобразовать интеграл в форму, где , в этой точке можно перейти идентично квадратуре Кленшоу – Кертиса для f, как указано выше. Однако из-за особенностей концевых точек в этом преобразовании координат используется первое квадратурное правило Фейера [которое не вычисляет f (−1)], если g (∞) не является конечным.

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

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

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

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

и это

Как обсуждалось выше, из-за наложения спектров нет смысла вычислять дополнительные коэффициенты , поэтому D является матрицей. В терминах этих коэффициентов c интеграл приблизительно равен:

сверху, где c - вектор вышеуказанных коэффициентов, а d - вектор интегралов для каждого коэффициента Фурье:

(Обратите внимание, однако, что эти весовые коэффициенты изменяются, если изменить матрицу DCT D, чтобы использовать другое соглашение о нормализации. Например, обычно DCT типа I определяют с дополнительными коэффициентами 2 или 2 в первом и последние строки или столбцы, что приводит к соответствующим изменениям в элементах d .) Суммирование может быть преобразовано в:

где w - вектор желаемых весов, указанных выше, с:

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

Веса положительны и их сумма равна единице.

Смотрите также

Рекомендации

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