Gauss-kvadratúra

A Wikipédiából, a szabad enciklopédiából

A numerikus analízisben használt kvadratúra fogalma egy adott függvény integráljának a közelítése, általában az integrálási tartományból választott függvényértékek súlyozott átlagával. A Carl Friedrich Gaussról elnevezett Gauss-kvadratúra olyan kvadratúra, amely pontos értéket ad 2n − 1 vagy alacsonyabb fokú polinomok esetén az xi pontok és wi súlyok megfelelő megválasztása esetén (i = 1,...,n).

Az integrálási tartományt általában [−1, 1]-nek véve, a kvadratúra

\int_{-1}^1 f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i).

A fenti Gauss-kvadratúra csak akkor szolgáltat pontos értékeket, ha az f(x) függvény polinommal jól közelíthető a [-1,1] intervallumban.

Általában használt súlyok például: W(x)=(1-x^2)^{-1/2}\, (Gauss-Csebisev) és W(x)=e^{-x^2} (Gauss-Hermite).

Gauss-kvadratúrák és ortogonális polinomok[szerkesztés | forrásszöveg szerkesztése]

A Gauss-kvadratúrák lényege, hogy mi magunk válasszuk meg nem csak a súlyokat, hanem az abszcisszákat is, ahol a függvényt meg szeretném közelíteni. A Gauss-kvadratúrák foka kétszer magasabb lesz mint a Newton-Cotes formulák foka, ugyanannyi függvényértékelés esetén. Azonban magasabb fok nem mindig jelent nagyobb pontosságot is. Nagyobb fok akkor jelent nagyobb pontosságot amikor az integrálandó függvény sima és jól meg lehet közelíteni egy polinommal. A Gauss kvadratúrák előnye, hogy adott W(x) és N egész szám esetén találunk egy rend wj súlyt és xj abszcisszát úgy, hogy a megközelítés

 \int_a^b f(x)dx\approx\sum_{i=0}^N w_i f_i (1)

,egzakt ha f(x) polinomiális.

például \int_a^b \frac{exp(-\cos^2{x})}{\sqrt{1-x^2}}

Ebben az esetben

W(x)=\frac{1}{\sqrt{1-x^2}}

Ez a partikuláris megválasztása súlyfüggvénynek a Gauss-Csebisev integrálás. Az integrálandó kifejezést úgy is fel lehet írni, hogy a súlyfüggvény W(x) nem látható egyértelműen: Legyen g(x)=W(x)f(x) és vj=wj / w(xj) és akkor (1)-es így alakul:

 \int_a^b g(x)dx\approx\sum_{j=1}^N v_j g(x_j) (2)

Amikor a súlyokat és az abszcisszákat keressük egy adott W(x) függvény esetén jól meg kell gondolnunk, hogy az (1)-es a vagy a (2)-es képletet szeretnénk alkalmazni.

A Gauss-kvadratúrák elmélete visszanyúlik egészen Gaussig 1814-re, aki folytonos törteket használt, hogy kifejlessze a témát. 1826-ban Jacobi újraértékelte Gauss eredményeit az ortogonális polinomok segítségével. Az ortogonális polinomok bevezetéséhez lerögzítünk egy (a,b) intervallumot, és bevezetjük egy f és g függvény skaláris szorzatát egy súlyfüggvényen W-n:

\langle f|g \rangle=\int_a^b W(x)f(x)g(x)dx

A skaláris szorzat egy szám és nem x-nek egy függvénye. Két függvény ortogonális ha a skaláris szorzatuk nulla. Egy függvény normalizált ha az önmagával való skaláris szorzata egységnyi. Egy olyan polinomrendszert, amely mind ortogonális és egyenként normalizált polinomokat tartalmaz ortonormált polinomrendszernek nevezzük.

Ortogonális polinomok generálása[szerkesztés | forrásszöveg szerkesztése]

Inicializálunk majd generálunk:

p_{-1}(x)=0
p_0(x)=1
p_{j+1}(x)=(x-a_j)p_j(x)-b_jp_{j-1}(x), j=1,2,… (3)

,ahol

a_j=\frac{\langle xp_j|p_j\rangle}{\langle p_j|p_j\rangle} , j=0,1,..
b_j=\frac{\langle p_j|p_j\rangle}{\langle p_{j-1}|p_{j-1}\rangle}, j=1,2,.. (4)
b0-át vehetjük 0 nak

A (3)-as képletben definiált polinomok főpolinomok, vagyis a főegyütthatója [p1(x)-nek az x^j] egységnyi. Ha elosztunk minden pj(x) polinomot \mid\langle p_j|p_j \rangle|^{1/2} skalárral akkor egy ortonormális polinomrendszert kapunk. Kimutatható, hogy minden pj polinomnak pontosan j különböző gyöke van az (a,b) intervallumban. Továbbá kimutatható, hogy pj(x) polinom gyökei közreékelődnek a pj-1(x) polinom j-1 gyöke közé úgy, hogy két szomszédos pj-1(x) gyök között van egy pj(x) gyök. Ez akkor nagyon fontos amikor meg akarjuk találni a polinom összes gyökét. Az összes gyököt azért szeretnénk megtalálni mert az N pontos Gauss kvadratúra formula abszcisszái pontosan a PN(x) polinom gyökei. Természetesen ugyanazon intervallum és súlyfüggvény esetén. Miután ismerjük x1,…,xN abszcisszákat meg kell keresnünk a wj súlyfüggvényt j=1,2,…,N esetén. Az egyik módszer a súlyok kiszámolására:

\begin{pmatrix}
  p_0(x_1) & \dots & p_0(x_N \\
  p_1(x_1) & \dots & p_1(x_N)\\
  \vdots & \vdots & \vdots\\
  p_{N-1}(x_1) & \dots & p_{N-1}(x_N)
\end{pmatrix} \begin{pmatrix}
  w_1\\
  w_2\\
  \vdots\\
  w_n
\end{pmatrix}=\begin{pmatrix}
  \int_a^b W(x)p_0(x)dx\\
  0\\
  \vdots\\
  0
\end{pmatrix}

Megoldva az egyenletrendszert megkapjuk a wj-ket. Ezen wj-k megfelelők az első N ortogonális polinom esetén. Ki lehet mutatni, hogy az integrálja következő N-1 polinomnak is egzakt, tehát a kvadratúra egzakt minden olyan polinomra amelynek foka 2N-1 -nél kisebb. Egy másik módszer a súlyok kiszámítására:

w_j=\frac{\langle p_{N-1}|p_{N-1}\rangle}{p_{N-1}(x_j)p'_N(x_j)}(5),ahol

a p'_N{x_j} az ortogonális polinom deriváltja.

Összefoglalva a Gauss kvadratúrák kiszámítása 2 különböző fázisból áll.

-a p0,…,pN ortogonális polinomok generálása, aj, bj együtthatók kiszámolása
-meghatározni a pN(X) gyökeit és kiszámolni a nekik megfelelő súlyokat
Klasszikus ortogonális polinomok esetében aj, bj együtthatók ismertek így az első lépés kihagyható. Viszont ha nem egy klasszikus súlyfüggvényünk van akkor az ortogonális polinomok felépítése nem olyan egyszerű.

A súlyok és az abszcisszák kiszámolása[szerkesztés | forrásszöveg szerkesztése]

Ez a művelet változhat könnyűtől egészen nehézig annak függvényében, hogy mennyit tudunk súlyfüggvényekről és a hozzuk rendelt polinomokról. Klasszikus esetben, tanulmányozott ortogonális polinomok esetében gyakorlatilag minden ismert még a polinomok zérus helyei is jól meg vannak közelítve. Ezt használjuk mint kezdeti megközelítést, lehetővé téve, hogy a Newton módszer gyorsan konvergáljon. A Newton módszer megköveteli a P'N(x) derivált ismeretét, amit a PN PN-1 függvényében fejezünk ki.
A következőkben itt vannak azok a súlyfüggvények, intervallumok, rekurencia képletek/relációk amelyek generálják a leggyakrabban használt ortogonális polinomokat és a nekik megfelelő Gauss kvadratúra formulákat

Gauss-Legendre:
 W(x)=1 ,-1<x<1
 (j+1)P_{j+1}=(2j+1)xP_j-jP_{j-1}
Gauss-Chebyshev:
 W(x)=(1-x^2)^{-1/2},-1<x<1
 T_{j+1}=2xT_j-T_{j-1}
Gauss-Laguerre:
 W(x)=x^{\alpha}e^{-x},0<x<\infin
(j+1)L_{j+1}^{\alpha}=(-x+2j+{\alpha}+1)L_j^{\alpha}-(j+\alpha)L_{j-1}^{\alpha}
Gauss-Hermite:
 W(x)=e^{-x^2},-\infin<x<\infin(5)
 H_{j+1}=2xH_j-2jH_{j-1}

Most egyedi módszereket adunk arra, hogy hogyan számoljuk az abszcisszákat és súlyokat ezekben az esetekben. Először a leggyakoribb abszcisszák és súlyok halmaza, a Gauss-Legendre.

 w_j=\frac{2}{(1-x_j^2)[P'_N(x_j)]^2}

A módszer az (x1,x2) integrálási határt átviszi (-1,1)-be és xj-ket és wj-ket nyújt a Gauss formuláknak.

\int_{x_1}^{x_2}f(x)dx=\sum_{j=1}^{N}w_jf(x_j)

A következő 3 módszer kezdetben megközelíti a gyököket. Az első a Gauss-Laguerre abszcisszák és súlyok, amelyeket a következő integrálási formulával használunk:

 \int_{0}^{\infin}x^{\alpha}e^{-x}f(x)dx=\sum_{j=1}^{N}w_jf(x_j)

A következő módszer a Gauss-Hermite abszcisszák és súlyok esetén használható. Ha ezen függvények standard normalizálását használjuk, ahogy az meg van adva az (5)-ös képletben, nagy N-ek esetén a számolás túlcsordul a különböző faktorok megjelenése miatt. Ezt elkerülhetjük ha egy ortogonális polinomrendszert használunk \overline{H}_j
. A rekurencia a következő lesz:

\overline{H}_{-1}=0
\overline{H}_0=\frac{1}{\pi^{1/4}}
\overline{H}_{j+1}=x\sqrt{\frac{2}{j+1}}\overline{H}_j-\sqrt{\frac{j}{j+1}}\overline{H}_{j-1}

A súly képlete így alakul

w_j=\frac{2}{(\overline{H}'_j)^2}

a derivált képlete pedig

\overline{H}'_j=\sqrt{2j}\overline{H}_{j-1}.

És végül az integrálandó formula

\int_{-\infin}^{\infin} e^{-x^2}f(x)dx=\sum_{j=1}^N w_jf(x_j)

Ismerve a rekurenciákat[szerkesztés | forrásszöveg szerkesztése]

A továbbiakban azt az esetet tárgyaljuk, amikor nem tudunk jó megközelítést adni az ortogonális polinomok zérus helyeire,de ismerjük az aj, bj együtthatókat amik generálják őket. Amint láttuk az előbbiekben a pN(x) gyökei az N pontos Gauss kvadratúra abszcisszái. A leghasznosabb összefüggés a súlyok kiszámolására az (5)-ös. A p'N deriváltat ki lehet számolni a (3)-as képlet lederiválásával. A (3)-as összefüggés viszont csak monikus polinomok esetében igaz, más normalizálások esetén bejön még egy plusz faktor \lambda_{N}/\lambda_{N-1}, ahol \lambda_{N} az x^N együtthatója a pN polinomban.
Kivéve azokat a speciális eseteket amelyeket már letárgyaltunk, a Newton módszer nem a legjobb gyökkeresési módszer. Sokkal gyorsabb ha a Galub-Welsch algoritmust használjuk, amely Wilf eredményeire épül. Az algoritmus szerint ha a (3)-as összefüggésben az xpj tagokat bal oldalra visszük és a p{j+1}-es tagokat a jobb oldalra akkor a rekurencia reláció felírható mátrix alakban:

 x\begin{pmatrix}
  p_0\\
  p_1\\
  \vdots\\
  p_{N-2}\\
  p_{N-1}
\end{pmatrix}=\begin{pmatrix}
  a_0 & 1 & \\
  b_1 & a_1 & 1 & \\
   & \vdots & \vdots\\
   & & b_{N-2} & a_{N-2} & 1\\
   & & & b_{N-1} & a_{N-1}
\end{pmatrix} \begin{pmatrix}
  p_0\\
  p_1\\
  \vdots\\
  p_{N-2}\\
  p_{N-1}
\end{pmatrix}+\begin{pmatrix}
  0\\
  0\\
  \vdots\\
  0 \\
  p_N
\end{pmatrix}

vagy

x p = T p + pN eN-1
(6)

Itt a T tridiagonális mátrix, p a p0,p1,…,pN-1 oszlop vektor és eN-1 egy egységvektor egy 1-sel az (N-1) pozícióban és zérusok minden más pozícióban. A T mátrix szimmetrizálható egy diagonális hasonlósági transzformációval D:

J = D T T^{-1}= \begin{pmatrix}
  a_0 & \sqrt{b_1} & \\
  \sqrt{b_1} & a_1 & \sqrt{b_2} & \\
   & \vdots & \vdots\\
   & & \sqrt{b_{N-2}} & a_{N-2} & \sqrt{b_{N-1}}\\
   & & & \sqrt{b_{N-1}} & a_{N-1}
\end{pmatrix}

A J mátrixot Jacobi-mátrixnak nevezzük. A (6)-ból láthatjuk, hogy p_N(x_j)=0 ekvivalens azzal, hogy xj sajátértéke a T -nek. Amíg a sajátértékeket egy hasonlósági transzformáció szolgáltatja addig xj sajátértéke a J szimmetrikus tridiagonális mátrixnak. Továbbá Wilf rámutatott arra, hogy ha vj megfelelő sajétvektora az xj sajátértéknek v v=1 és akkor

w_j=\mu_0v_{j,1}

ahol

\mu_0=\int_a^b W(x)dx

ahol vj,1 a v első komponense.

Ortogonális polinomok nem klasszikus súlyokkal[szerkesztés | forrásszöveg szerkesztése]

A Stieltjes eljárás kiszámolja az a0-át és p1(x)-et a (3)(4)-ből
Ismerve p0 és p1-et ki tudjuk számolni a1-et és b1-et a (4)-ből. Most ki kell számolnunk a skaláris szorzatot! p_j(x) polinomiális x-ben, számolom a skaláris szorzatot de csak ha ismerem az első 2N momentumát a súlyfüggvénynek :

\mu_j=\int_a^b x^jW(x)\,dx , j=0,1,...,2N-1

Ezt nem használjuk mert nagyon pontatlan.
Suck és Donován rájöttek, hogy a numerikus stabilitás nagyon megjavul, ha az x^n alapú polinomrendszert helyett használjunk egy ismert ortogonális polinomrendszert a \pi_j(x)-et
-A változtatott momentum

\nu_j=\int_a^b \pi_jW(x)dx , j=0,1,…,2N-1

ahol \pi_j kielégíti analóg módon a (3)(4)-et

\pi_{-1}(x)=0
\pi_0(x)=1
\pi_{j+1}(x)=(x-\alpha_j)\pi_j(x)-\beta_j\pi_{j-1}(x) ,j=1,2,...

ahol \alpha_j, \beta_j természetesen ismertek.
Majd Wheeler adott egy hatékony algoritmust O(N^2):

\sigma_{k,l}=\langle p_k|\pi_l\rangle , k,l>-1

Inicializálunk:

\sigma_{-1,l}=0 ,l=1,...,2N-2
\sigma_{0,l}=\nu_,l=0,1,...,2N-1
a_0=\alpha_0+\frac{\nu_1}{\nu_0}
b_0=0

Majd k=1,2,..,N-1 számoljuk:

\sigma_{k,l}=\sigma_{k-1,l+1}-(a_{k-1}-\alpha_l)\sigma_{k-1,l}-b_{k-1}\sigma_{k-2,l}+\beta_l\sigma_{k-1,l-1},l=k,k+1..,2N-1
a_k=\alpha_k-\frac{\sigma_{k-1,k}}{\sigma_{k-1,k-1}}+\frac{\sigma_{k,k+1}}{\sigma_{k,k}}
b_k=\frac{\sigma_{k,k}}{\sigma_{k-1,k-1}}

A Wheeler-algoritmus megköveteli, hogy a megváltoztatott momentum legyen pontosan kiszámítva. Az algoritmus nagyon hatásos véges (a,b) intervallum esetén, végtelen intervallum esetén az algoritmus nem távolítja el a hibák nagy részét. Ebben az esetben Gautschi azt javasolja, hogy csökkentsük le az intervallumot egy véges intervallumra változócsere segítségével.

Lásd még[szerkesztés | forrásszöveg szerkesztése]

Források[szerkesztés | forrásszöveg szerkesztése]

  • Gisbert Stoyan - Takó Galina: Numerikus módszerek / Gisbert Stoyan ; Takó Galina. 2. átdolg. kiad. Budapest : Typotex Elektronikus Kiadó Kft., 2002. 439 p. : ill. ISBN 963-9326-20-8

Külső hivatkozások[szerkesztés | forrásszöveg szerkesztése]