Kwadratura Clenshawa-Curtisa

Kwadratura Clenshawa – Curtisa i kwadratura Fejéra to metody całkowania numerycznego lub „kwadratury”, które opierają się na rozwinięciu całki pod względem wielomianów Czebyszewa . Równoważnie wykorzystują zmianę zmiennych $($ używają przybliżenia dyskretnej transformacji kosinusowej DCT dla Oprócz szybko zbieżnej dokładności porównywalnej z regułami kwadratury Gaussa , kwadratura Clenshawa-Curtisa w naturalny sposób prowadzi do zagnieżdżonych reguł kwadraturowych (gdzie różne rzędy dokładności mają wspólne punkty), co jest ważne zarówno w przypadku kwadratury adaptacyjnej , jak i kwadratury wielowymiarowej ( kubatura ).

W skrócie, funkcja $, która ma być$$zintegrowana,$ jest oceniana w ekstremach lub pierwiastkach używane do skonstruowania wielomianowego przybliżenia funkcji Ten wielomian jest następnie całkowany dokładnie. W praktyce wagi integracji dla wartości funkcji w każdym węźle są wstępnie obliczane, a obliczenie to $w$ szybkiej Fouriera -powiązane algorytmy dla DCT.

Metoda ogólna

Prostym sposobem zrozumienia algorytmu jest uświadomienie sobie, że kwadratura Clenshawa-Curtisa (zaproponowana przez tych autorów w 1960 r.) sprowadza się do całkowania poprzez zmianę zmiennej x = cos( θ ) . Algorytm jest zwykle wyrażany dla całkowania funkcji f ( x ) w przedziale [−1,1] (dowolny inny przedział można uzyskać przez odpowiednie przeskalowanie). Dla tej całki możemy napisać:

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

$\ Displaystyle f (\ cos \ theta) \ sin$ przekształciliśmy problem z całkowania na problem całkowania $θ$ . Można to wykonać, jeśli znamy szereg kosinusowy dla ${\ Displaystyle f (\ cos \ theta)}$ :

$${\ Displaystyle f (\ cos \ teta) = {\ Frac {a_ {0}} {2}} + \ suma _{k=1}^{\infty}a_{k}\cos(k\theta)}$$

w takim przypadku całka staje się:

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

Oczywiście w celu obliczenia współczynników serii cosinus

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

trzeba ponownie wykonać całkowanie numeryczne, więc na początku może się wydawać, że nie uprościło to problemu. Jednak w przeciwieństwie do obliczania dowolnych całek, całki w szereg Fouriera dla $displaystyle k = N} , są$$konstrukcję)$ jak przez aż do Nyquista $dokładnie$ obliczane przez punkty dla $\ Displaystyle \ theta _ {n} = n \ pi / N$ } ${\ Displaystyle n = 0, \ ldots, N}$ (z wyjątkiem tego, że punkty końcowe są ważone przez 1/2, aby uniknąć podwójnego liczenia, co jest równoważne z regułą trapezów lub wzorem Eulera – Maclaurina ). Oznacza to, że aproksymujemy całkę cosinusową za pomocą dyskretnej transformaty kosinusowej typu I (DCT):

$${\ Displaystyle a_ {k} \ około {\ Frac {2} {N}} \ lewo [{\ Frac {f (1)} {2}} + {\ Frac {f (-1)} {2 }}(-1)^{k}+\suma _{n=1}^{N-1}f(\cos[n\pi /N])\cos(nk\pi /N)\prawo]}$$

dla $}$$displaystyle a_ {k}$ następnie użyj powyższego wzoru dla całki w kategoriach tych za . Ponieważ potrzebny jest tylko ${\ displaystyle a_ {2k}}$ , wzór upraszcza się dalej do DCT typu I rzędu N / 2 , zakładając, że N jest parzystą :

$${\ Displaystyle a_ {2k} \ około {\ Frac {2} {N}} \ lewo [{\ Frac { f(1)+f(-1)}{2}}+f(0)(-1)^{k}+\suma _{n=1}^{N/2-1}\lewo\{f (\cos[n\pi /N])+f(-\cos[n\pi /N])\right\}\cos \left({\frac {nk\pi }{N/2}}\right )\Prawidłowy]}$$

Z tego wzoru jasno wynika, że ​​​​reguła kwadratur Clenshawa-Curtisa jest symetryczna, ponieważ równo waży f ( x ) i f ( - x ).

Z powodu aliasingu współczynniki oblicza się tylko $do$ k N / 2 , ponieważ dyskretne próbkowanie funkcji powoduje, że częstotliwość 2 k nie do odróżnienia od częstotliwości N k . Równoważnie, to amplitudy unikalnego wielomianu interpolacji trygonometrycznej o ograniczonym paśmie przechodzącym przez punkty N + 1 , w których oblicza się f ( cos θ ) i przybliżamy całkę przez całkę tego za 2 k {\ displaystyle a_ $}$ wielomian interpolacyjny. Istnieje jednak pewna subtelność w sposobie traktowania $jak$ w całce - aby uniknąć podwójnego liczenia z jego aliasem, jest on uwzględniany z wagą 1/2 w końcowej przybliżonej całce ( można również zobaczyć, badając wielomian interpolujący):

$${\ Displaystyle \ int _ {0} ^ {\ pi} f (\ cos \ teta) \ sin (\ teta) \, d \ teta \ około a_ {0} + \ suma _ {k = 1} ^ {N /2-1}{\frac {2a_{2k}}{1-(2k)^{2}}}+{\frac {a_{N}}{1-N^{2}}}.}$$

Połączenie z wielomianami Czebyszewa

Powodem, dla którego jest to połączone z wielomianami Czebyszewa, $Displaystyle T_ {k} (\ cos \ theta ) = \ cos (k \ theta )}$$θ$ definicji , a więc powyższy szereg cosinusów jest w rzeczywistości przybliżeniem wielomianów Czebyszewa : ${\ Displaystyle f (x)}$

$${\ Displaystyle f (x) = {\ Frac {a_ {0}} {2}} T_ {0} (x)+\suma _{k=1}^{\infty}a_{k}T_{k}(x),}$$

a zatem „naprawdę” integrujemy, $rozwinięcie$ pod względem wielomianów Czebyszewa Punkty $_$ wielomianu $_ N}(x)}$ T_ _ .

Fakt, że takie $, jest odpowiedzialny za szybką$ Czebyszewa przybliżenia, ponieważ uwzględniono więcej . Szereg kosinusowy zbiega się bardzo szybko dla funkcji, które są parzyste , okresowe i wystarczająco gładkie. Jest to prawdą tutaj, ponieważ $,$ parzysty i okresowy ${\ Displaystyle f (x)}$$jest$ i - różniczkowalny jest k - razy różniczkowalna na ${\ Displaystyle [-1,1]}$ . (W przeciwieństwie do tego, bezpośrednie zastosowanie rozwinięcia $\ Displaystyle f (\ cos \ theta)}$ szereg kosinusowy do zamiast zamiast $)$ zbieżne szybko , ponieważ nachylenie parzystookresowego wydłużenia byłoby na ogół nieciągłe).

Kwadratura Fejéra

Fejér zaproponował dwie zasady kwadraturowe bardzo podobne do kwadratury Clenshawa-Curtisa, ale znacznie wcześniej (w 1933 r.).

Z tych dwóch „druga” reguła kwadraturowa Fejéra jest prawie identyczna z regułą Clenshawa – Curtisa. Jedyna różnica $są$$że$ punkty i Oznacza to, że Fejér użył tylko ekstremów wewnętrznych wielomianów Czebyszewa, tj. rzeczywistych punktów stacjonarnych.

$punktów , w połowie$$kwadratury$ Fejéra ocenia wartość w między ekstremami : ${\ Displaystyle \ theta _ {n} = (n + 0,5) \ pi / N}$ dla ${\ Displaystyle 0 \ równoważnik n <N}$ . To są korzenie $_$ _ _ _ (Te równomiernie rozmieszczone punkty środkowe są jedynym innym wyborem punktów kwadraturowych, które zachowują zarówno parzystą symetrię transformaty kosinusowej, jak i translacyjną symetrię okresowego szeregu Fouriera). Prowadzi to do wzoru:

$${\ Displaystyle a_ {k} \ około {\ frac {2}{N}}\suma _{n=0}^{N-1}f(\cos[(n+0,5)\pi /N])\cos[(n+0,5)k\pi / N]}$$

który jest właśnie DCT typu II. Jednak pierwsza reguła kwadraturowa Fejéra nie jest zagnieżdżona: punkty oceny dla 2 N nie pokrywają się z żadnym z punktów oceny dla N , w przeciwieństwie do kwadratury Clenshawa-Curtisa lub drugiej reguły Fejéra.

Pomimo faktu, że Fejér odkrył te techniki przed Clenshawem i Curtisem, nazwa „kwadratura Clenshawa – Curtisa” stała się standardem.

Porównanie do kwadratury Gaussa

Klasyczna metoda kwadratury Gaussa ocenia całkę w punktach i jest skonstruowana tak, aby $1 {$$displaystyle N + 1$ wielomiany do stopnia . $stopnia$$.$ tego, powyższa kwadratura Clenshawa – Curtisa ocenia całkę w punktach i integruje wielomiany tylko do Może się zatem wydawać, że Clenshaw-Curtis jest z natury gorszy niż kwadratura Gaussa, ale w rzeczywistości tak nie jest.

W praktyce kilku autorów zauważyło, że Clenshaw-Curtis może mieć dokładność porównywalną z kwadraturą Gaussa dla tej samej liczby punktów. Jest to możliwe, ponieważ większość całek liczbowych nie jest wielomianami (zwłaszcza, że ​​wielomiany można całkować analitycznie), a przybliżenie wielu funkcji za pomocą wielomianów Czebyszewa szybko się zbiega (patrz przybliżenie Czebyszewa ). $, że zarówno kwadratura$ Gaussa przez k dla -czasy całki różniczkowalnej.

$z często$ wagi kwadraturowe można oszacować w czasie za pomocą szybkich algorytmów Fouriera (lub dla DCT) $gdy$ algorytmów wag kwadraturowych Gaussa wymagała obliczenia. $Jednak$ algorytmy osiągnęły złożoność . Ze względów praktycznych całkowanie numeryczne wysokiego rzędu jest $rzadko$ wykonywane przez proste obliczenie wzoru kwadraturowego dla bardzo . Zamiast tego zwykle stosuje się adaptacyjny schemat kwadraturowy, który najpierw ocenia całkę do niskiego rzędu, a następnie sukcesywnie poprawia dokładność, zwiększając liczbę punktów próbkowania, prawdopodobnie tylko w obszarach, w których całka jest niedokładna. Aby ocenić dokładność kwadratury, porównuje się odpowiedź z regułą kwadratury jeszcze niższego rzędu. W idealnym przypadku ta reguła kwadraturowa niższego rzędu ocenia całkę w podzbiorze oryginalnych N punktów, aby zminimalizować oceny całki. Nazywa się to zagnieżdżoną regułą kwadraturową i tutaj Clenshaw-Curtis ma tę zaletę, że reguła dla rzędu N wykorzystuje podzbiór punktów z rzędu 2 N . W przeciwieństwie do tego, reguły kwadraturowe Gaussa nie są naturalnie zagnieżdżone, dlatego należy zastosować formuły kwadraturowe Gaussa-Kronroda lub podobne metody. Reguły zagnieżdżone są również ważne dla rzadkich siatek w wielowymiarowej kwadraturze, a kwadratura Clenshawa-Curtisa jest popularną metodą w tym kontekście.

Integracja z funkcjami wagowymi

$displaystyle f (x)}$ ogólnie, można postawić problem zintegrowania dowolnej ustaloną wagą, ${$ jest znana z wyprzedzeniem: fa

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

Najczęstszym przypadkiem jest , $.$ jest inna funkcja wagi Podstawowym powodem jest to, że skoro $a$ wziąć pod uwagę priori , błąd całkowania zależeć tylko od dokładności przybliżenia $displaystyle f ( x)}$ , niezależnie od tego, jak źle zachowuje się funkcja wagi.

Kwadraturę Clenshawa – Curtisa można uogólnić na ten przypadek w następujący sposób. Tak jak poprzednio, działa to poprzez $znalezienie$ rozwinięcia w szereg cosinusowy DCT, a następnie całkowanie każdego wyrazu Teraz jednak całki te mają postać

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

Dla większości $całki$ nie można obliczyć analitycznie, w przeciwieństwie $,$$sama funkcja wagi$ zwykle używana dla wielu całek sobie pozwolić na wcześniejsze obliczenie tych z dużą dokładnością $,$ ponieważ $określony$ analitycznie , czasami można zastosować wyspecjalizowane

Na przykład opracowano specjalne metody stosowania kwadratury Clenshawa-Curtisa do całek postaci z funkcją wagi $x) w (x$$Displaystyle$ , która jest wysoce oscylacyjna, np. sinusoida lub funkcja Bessela (patrz np. Evans & Webster, 1999). Jest to przydatne w obliczeniach szeregów Fouriera i szeregów Fouriera-Bessela z dużą dokładnością , gdzie proste $\ displaystyle w (x) = 1}$ kwadraturowe są problematyczne ze względu na dużą dokładność wymaganą do rozwiązania wkładu szybkiego w ( oscylacje. Tutaj $podczas$ całki o szybkich oscylacjach jest brana pod uwagę za pomocą specjalistycznych metod dla $jest$ zwykle lepiej zachowywana

Innym przypadkiem, w którym funkcje wagowe są szczególnie przydatne, jest sytuacja, gdy całka jest nieznana, ale ma znaną osobliwość jakiejś postaci, np. znaną nieciągłość lub całkowalną rozbieżność (taką jak 1/ √ x ) w pewnym punkcie. $właściwości$ można wciągnąć do funkcji wagi, $jej$ wykorzystać do dokładnego obliczenia.

Należy zauważyć, że kwadraturę Gaussa można również dostosować do różnych funkcji wagowych, ale technika jest nieco inna. $jest$ zawsze oceniana w tym samym zestawie punktów, niezależnie od tego, co wielomianu Czebyszewa. W kwadraturze Gaussa różne funkcje wagi prowadzą do różnych wielomianów ortogonalnych , a tym samym do różnych pierwiastków, w których obliczana jest całka.

Całkowanie na przedziałach nieskończonych i półnieskończonych

Możliwe jest również użycie kwadratury Clenshawa – Curtisa do obliczenia całek postaci ${\ textstyle \ int _ {0} ^ {\ infty} f (x) \, dx}$ i ${\ textstyle \ int _ {- \ infty} ^ {\ infty} f (x) \, dx}$ , przy użyciu techniki ponownego mapowania współrzędnych. Wysoką dokładność, nawet zbieżność wykładniczą dla gładkich całek, można zachować tak długo, jak rozpada się wystarczająco szybko, ponieważ ${\ displaystyle f (x)}$ x | zbliża się do nieskończoności.

Jedną z możliwości jest użycie ogólnej transformacji współrzędnych, takiej jak x = t /(1− t ² )

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

przekształcić nieskończony lub półnieskończony interwał w skończony, jak opisano w Całkowanie numeryczne . Istnieją również dodatkowe techniki, które zostały opracowane specjalnie dla kwadratury Clenshawa-Curtisa.

Na $współrzędnych$ można użyć określoną jeden mógłby po prostu użyć L = 1; optymalny wybór L może przyspieszyć zbieżność, ale jest zależny od problemu), aby przekształcić całkę półnieskończoną w:

$${\ Displaystyle \ int _ {0} ^ {\ infty} f (x) \, dx = 2L \ int _ {0} ^ {\ pi} {\ Frac {f [L \ łóżeczko ^ {2} (\ teta /2)]}{[1-\cos(\theta)]^{2}}}\sin(\theta)\,d\theta.}$$

Czynnik mnożący sin( θ ), f (...)/(...) ² , można następnie rozwinąć w szereg kosinusowy (w przybliżeniu, stosując dyskretną transformatę kosinusową) i scałkować wyraz po wyrazie, dokładnie tak, jak to było zrobione dla f (cos θ ) powyżej. Aby wyeliminować osobliwość w θ = 0 w tej całce, wystarczy, że f ( x ) dąży do zera wystarczająco szybko, gdy x zbliża się do nieskończoności, aw szczególności f ( x ) musi zanikać co najmniej tak szybko, jak 1/ x ^3/2 .

$aby$ można użyć ponownego odwzorowania współrzędnych L jest przez użytkownika, jak powyżej integralna w:

$${\ Displaystyle \ int _ {- \ infty} ^ {\ infty} f (x) \, dx = L \ int _ {0} ^ {\ pi} {\ Frac {f [L \ łóżeczko (\ teta)] }{\sin ^{2}(\theta)}}\,d\theta \about {\frac {L\pi}}{N}}\sum _{n=1}^{N-1}{\frac {f[L\cot(n\pi /N)]}{\sin ^{2}(n\pi /N)}}.}$$

W tym przypadku wykorzystaliśmy fakt, że przemapowana całka f ( L cot θ )/sin ² ( θ ) jest już okresowa, więc może być bezpośrednio całkowana z dużą (nawet wykładniczą) dokładnością przy użyciu reguły trapezów (zakładając, że f jest wystarczająco gładkie i szybko psujące się); nie ma potrzeby obliczania szeregu cosinusowego jako kroku pośredniego. Zauważ, że reguła kwadraturowa nie obejmuje punktów końcowych, gdzie założyliśmy, że całka dąży do zera. Powyższy wzór wymaga, aby f ( x ) rozpadało się szybciej niż 1/ x ² , gdy x dąży do ±∞. (Jeżeli f rozpada się dokładnie jak 1/ x ² , to całka osiąga wartość skończoną w punktach końcowych i te granice muszą być uwzględnione jako wyrazy końcowe w regule trapezów.). Jeśli jednak f rozpada się szybko tylko wielomianowo, może być konieczne zastosowanie kolejnego kroku kwadratury Clenshawa-Curtisa, aby uzyskać dokładność wykładniczą ponownie odwzorowanej całki zamiast reguły trapezów, w zależności od większej liczby szczegółów ograniczających właściwości f : the problem polega na tym, że chociaż f ( L cot θ )/sin ² ( θ ) jest rzeczywiście okresowe z okresem π, niekoniecznie jest gładkie w punktach końcowych, jeśli wszystkie pochodne tam nie znikają [np. funkcja f ( x ) = tanh ( x ³ )/ x ³ rozpada się jako 1/ x ³ , ale ma skokową nieciągłość nachylenia przemapowanej funkcji przy θ=0 i π].

Zasugerowano inne podejście do ponownego mapowania współrzędnych dla całek postaci ${\ textstyle \ int _ {0} ^ {\ infty} e ^ {- x} g (x) \, dx }$ , w takim przypadku można użyć transformacji x $\ Displaystyle x = - \ ln [(1 + \ cos \ theta)/2]}$$\ textstyle \ int _ {0} ^ {\ pi} f (\ cos \ theta) \ sin \ theta \, d \ theta}$ grzech gdzie ${\ Displaystyle f (u) = g (- \ ln [(1 + u) / 2]) / 2} ,$ w którym punkt jeden można postępować identycznie z kwadraturą Clenshawa-Curtisa dla f jak powyżej. Jednak ze względu na osobliwości punktów końcowych w tym ponownym odwzorowaniu współrzędnych, stosuje się pierwszą regułę kwadraturową Fejéra [która nie ocenia f (−1)], chyba że g (∞) jest skończone.

Wstępne obliczanie wag kwadraturowych

W praktyce niewygodne jest wykonywanie DCT próbkowanych wartości funkcji f (cos θ) dla każdej nowej całki. Zamiast tego zwykle wstępnie oblicza się wagi kwadraturowe (dla $w_ {n}}$ od 0 do N / 2, zakładając, że N jest parzyste), więc w n displaystyle

$${\ Displaystyle \ int _ {-1} ^ {1} f (x) \, dx \ w przybliżeniu \ suma _ {n = 0} ^ {N/2} w_ {n} \ lewo \ {f (\ cos [ n\pi /N])+f(-\cos[n\pi /N])\prawo\}.}$$

$kategoriach$ te są również obliczane przez DCT, co łatwo zauważyć, wyrażając obliczenia w macierzowej . W szczególności obliczyliśmy współczynniki szeregu cosinusowego za pomocą wyrażenia w postaci: ${\ displaystyle a_ {2k}}$

$${\ Displaystyle c = {\ rozpocząć {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,}$$

gdzie D jest formą macierzową ( N /2+1)-punktowego DCT typu I z góry, z wpisami (dla indeksów od zera ):

$${\ Displaystyle D_ {kn} = {\ Frac {2 }{N}}\cos \left({\frac {nk\pi}}{N/2}}\right)\times {\begin{przypadki}1/2&n=0,N/2\\1&\mathrm { inaczej} \end{przypadki}}}$$

i jest ${\ displaystyle y_ {n}}$

$${\ Displaystyle y_ {n} = f (\ cos [n \ pi / N]) + f (- \ cos [n \ pi / N]).}$$

Jak omówiono powyżej, z powodu aliasingu nie ma sensu obliczanie współczynników poza $+$ więc D jest za $Displaystyle ( ( N/2+1)\razy (N/2+1)}$ macierz. Pod względem tych współczynników c całka wynosi w przybliżeniu:

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

z góry, gdzie c jest wektorem współczynników ${\ displaystyle a_ {2k}}$ , a d jest wektorem całek dla każdego współczynnika Fouriera: za

$${\ Displaystyle d = {\ rozpocząć {pmatrix} 1 \\ 2 / (1-4) \\ 2/ (1-16) \\\ vdots \\ 2/ (1- [N-2] ^ {2} )\\1/(1-N^{2})\end{pmacierz}}.}$$

(Należy jednak zauważyć, że te współczynniki wagi ulegają zmianie, jeśli zmieni się macierz DCT D , aby zastosować inną konwencję normalizacji. Na przykład często definiuje się DCT typu I z dodatkowymi czynnikami 2 lub √ 2 czynniki w pierwszym i ostatnie wiersze lub $wpisach$ we d .) można przeorganizować w następujący sposób:

$${\ Displaystyle \ int _ {-1} ^ {1} f (x) \ ,dx\około d^{T}c=d^{T}Dy=(D^{T}d)^{T}y=w^{T}y}$$

gdzie w jest wektorem pożądanych wag powyżej, przy czym ${\ displaystyle w_ {n}}$

$${\ Displaystyle w = D ^ {T} re.}$$

Ponieważ transponowana macierz jest również DCT (np. transponowana DCT typu I jest DCT typu I, prawdopodobnie z nieco inną normalizacją w zależności od zastosowanych konwencji D T {\ Displaystyle D ^ { $}$ ), wagi kwadraturowe w można wstępnie obliczyć w czasie O ( N log N ) dla danego N przy użyciu szybkich algorytmów DCT.

Wagi $jest$ równa jeden.

Zobacz też

• Formuła Eulera-Maclaurina
• Formuła kwadraturowa Gaussa-Kronroda