Metody numeryczne fizyki/Całkowanie numeryczne funkcji interpolacyjnej

Całkowanie numeryczne funkcji interpolacyjnej

---

Licencja

Autor: Mirosław Makowiecki Absolwent UMCS Fizyki Komputerowej Uniwersytetu Marii Curie-Skłodowskiej w Lublinie Email: miroslaw(kropka)makowiecki(małpa)gmail(kropka)pl Dotyczy: książki, do której należy ta strona, oraz w niej zawartych stron i w nich podstron, a także w nich kolumn, wraz z zawartościami. Użytkownika książki, do której należy ta strona, oraz w niej zawartych stron i w nich podstron, a także w nich kolumn, wraz z zawartościami nie zwalnia z odpowiedzialności prawnoautorskiej nieprzeczytanie warunków licencjonowania. Umowa prawna: Creative Commons: uznanie autorstwa, na tych samych warunkach, z możliwością obowiązywania dodatkowych ograniczeń. Autor tej książki dołożył wszelką staranność, aby informacje zawarte w książce były poprawne i najwyższej jakości, jednakże nie udzielana jest żadna gwarancja, czy też rękojma. Autor nie jest odpowiedzialny za wykorzystanie informacji zawarte w książce, nawet jeśli wywołaby jakąś szkodę, straty w zyskach, zastoju w prowadzeniu firmy, przedsiębiorstwa lub spółki bądź utraty informacji, niezależnie czy autor (a nawet Wikibooks) został powiadomiony o możliwości wystąpienie szkód. Informacje zawarte w książce mogą być wykorzystane tylko na własną odpowiedzialność.

---

---

W tym rozdziale zapoznamy się z przybliżonym obliczanie całek znając tylko funkcje podcałkową. Często obliczenie funkcji pierwotnej jest uciążliwe, a nawet niemożliwe, bo podana jest za pomocą tablicy wartości funkcji podcałkowej, wtedy definicja funkcji pierwotnej traci sens. Mając tablicę funkcji można napisać wielomian interpolacyjny, który łączy te punkty, i napisać całkowanie tego wielomianu. Niech mamy wielomian interpolacyjny Lagrange'a określony wzorem (1.7) (gdzie y_k jest to właśnie funkcja F(x_k)) przy definicji Φ_j(x) (1.11), wtedy możemy napisać przybliżoną równość pomiędzy całką dokładną, a przybliżoną, którą piszemy w postaci:

$\int _{a}^{b}F(x)dx\simeq \int _{a}^{b}\varphi (x)dx=\int _{a}^{b}\sum _{k=1}^{n}\Phi _{k}(x)F(x_{k})=\sum _{k=1}^{n}A_{k}F(x_{k}){\mbox{ gdzie: }}A_{k}=\int _{a}^{b}\Phi _{k}(x)dx\;$

(4.1)

Bład bezwzględny pomiędzy wartością funkcji, a jej wartością przybliżoną, możemy przepisać przy pomocy schematu:

$|F(x)-\varphi (x)|<\epsilon {\mbox{, dla }}x\in \langle a,b\rangle \;$

(4.2)

to wtedy wykorzystując (4.2) i biorąc za funkcję φ(x) wielomian interpolacyjny, wtedy otrzymujemy otrzymujemy oszacowanie błędu popełnionego przy liczenia całki, który jest jakoby przybliżoną wartością funkcji interpolowanej:

$\left|\int _{a}^{b}\left[F(x)-\varphi (x)\right]dx\right|=\left|\int _{a}^{b}F(x)dx-\sum _{k=0}^{n}A_{k}F(x_{k})\right|\leq \epsilon (b-a)\;$

(4.3)

Za pomocą funkcji interpolowanej możemy policzyć całkę funkcji interpolowanej przy pomocy wielomianu interpolacyjnej Lagrange'a z dowolną dokładnością. Niech funkcja F(x) posiada pewne osobliwości, która utrudnia dobre przybliżenie jej przy pomocy funkcji interpolowanej, więc podzielmy ten wielomian na iloczyn dwóch funkcji, tzn. funkcji p(x) i f(x), gdzie f(x) jest funkcją, którą dobrze jest przybliżyć, a p(x) zawiera wszelkie osobliwości funkcji osobliwej F(x). Niech φ(x) będzie wielomianem interpolacyjnym funkcji f(x), wtedy przybliżenie funkcji interpolowanej f(x) funkcją interpolowaną φ(x) daje w wyniku czego:

$\int _{a}^{b}F(x)dx=\int _{a}^{b}p(x)f(x)dx\simeq \int _{a}^{b}p(x)\varphi (x)dx=\sum _{k=0}^{n}A_{k}f(x_{k})\;$

(4.4)

Współczynniki A_k określonego na przedziale <a,b> przy pomocy funkcji Φ_k(x) jako całka przy normie p(x) piszemy:

$A_{k}=\int _{a}^{b}p(x)\Phi _{k}(x)dx\;$

(4.5)

Ten sposób można wykorzystać do obliczania całek nieskończonych, wtedy współczynnik p(x) określa szybkość malenia funkcji F(x) dla x→±∞.

Wprowadzenie do wzorów całkowania numerycznego

W tym rozdziale będziemy przybliżać obliczanie całek, którego funkcją podcałkową jest f(x), a całka jest wykonywana przy wadze p(x), i ją przedstawiamy:

$I(f)=\int _{a}^{b}p(x)f(x)dx\;$

(4.6)

i będziemy je przybliżać w pewnym przedziale <a,b>, gdy obierzemy punkty x_k i będziemy je stosowali do funkcji f(x) przy stałych A_k podobnie zdefiniowanej jak w punkcie (4.5), wtedy mamy:

$S(f)=\sum _{k=1}^{n}A_{k}f(x_{k}){\mbox{ dla }}x_{k}\in \langle a,b\rangle \;$

(4.7)

Punkty x_K nazywamy węzłami kwadratury, a wzór (4.7) kwadraturą. Wzory te będziemy stosowali też do obliczania całek, które mają pewne osobliwości określone przez wagę p(x), zakładamy, że waga jest funkcją zależną od zmiennej "x" i przyjmuje wartości nieujemną. Błąd przybliżenia całki określonej przez funkcję interpolacyjną, a interpolowaną określamy przy pomocy wzorów (4.6) i (4.7) przy definicji błędu wzoru funkcji interpolującej (1.24), wtedy:

$E(f)=I(f)-S(f)={{f^{(n+1)}(\xi )} \over {(n+1)!}}\int _{a}^{b}p(x)\omega _{n}(x)dx\;$

(4.8)

Jako kryterium dokładności kwadratury będziemy pisali S(W) i I(W), gdy W jest pewnym wielomianem. Kwadratura (4.7) nazywamy kwadraturą rzędu "r", jeśli I(W)=S(W), gdy wielomian W(x) jest stopnia mniejszego niż "r", a jeśli I(W)≠S(W), gdy wielomian W(x) jest stopnia r≥1.

Twierdzenie
Kwadratura, która określimy przez (4.7) jest zbieżna do dowolnej funkcji f∈C(<a,b>), to wtedy warunek konieczny i dostateczny jest taki, że ona jest zbieżna do dowolnego wielomianu i istnieje liczba M dla współczynników $A_{k}^{n}\;$ niezależna od ilości węzłów "n":

$B_{n}=\sum _{k=0}^{n}\left|A_{k}^{n}\right|\leq M{\mbox{ dla }}n=1,2,...\;$

(4.9)

Całkowanie wielomianu interpolacyjnego z ustalonymi punktami

Będziemy tutaj rozpatrywali, że mamy ustalone węzły, tzn. x_k dla k=0,1,2,..,n. Poniżej będziemy liczyli współczynniki A_k, by wyznaczyć wartość całki wielomianu interpolującego jako przybliżenie całki funkcji interpolowanej w przedziale (a,b).

Metoda całkowa wielomianu interpolacyjnego Newtona-Cotesa

Będziemy się tutaj zajmowali kwadraturą Newtona-Cotesa określoną przez wzór początkowy (4.7), w którym zdefiniujemy A_K sposobem (4.1), a wielomian Ψ_k(x) jest zdefiniowany przez (1.11), w którym błąd kwadratury jest zdefiniowaliśmy poprzez wzór (1.24). Całkowanie wielomianu Φ_k(x) w przedziale (a,b) jest napisane jako definicję A_k:

$A_{k}=\int _{a}^{b}\Phi _{k}(x)dx=h{{(-1)^{n-k}} \over {k!(n-k)!}}\int _{0}^{n}{{t(t-1)...(t-n)} \over {t-k}}{\mbox{ gdzie: }}h={{b-a} \over {n}}{\mbox{, }}f_{k}=f(a+kh)\;$

(4.10)

Twierdzenie
Warunkiem dostatecznym i koniecznym dla kwadratury (4.7) jest jego rząd minimum n+1

Dowód powyższego twierdzenia
Niech wielomian f(x) będzie stopnia "n" i mamy "n" węzłów, wtedy f(x)=L_n(x) i I(f)=I(L_n)=S(f), wtedy rząd kwadratury wynosi n+1. Załóżmy, że rząd kwadratury wynosi co najmniej n+1, wtedy jest ona dokładna dla wielomianów stopnia mniejszego niż n+1, co stąd wynika, że jest dokładna dla wielomianów f(x), czyli kwadraturę określamy poprzez wzór (4.7).

Aby policzyć różne współczynniki A_k należy skorzystać ze wzoru (4.10), wtedy poszczególne wspomniane współczynniki dla n=1 liczymy jako:

$A_{0}=-h\int _{0}^{1}(t-1)={{1} \over {2}}h\;$

(4.11)

$A_{1}=h\int _{0}^{1}tdt={{1} \over {2}}h\;$

(4.12)

Wtedy kwadraturę (4.6) przy pomocy współczynników (4.11) i (4.12) piszemy przy pomocy wartości funkcji f₁ i f₂ wyznaczonych w punktach "a" i "b":

$S(f)={{1} \over {2}}(f_{0}+f_{1})\;$

(4.13)

Błąd kwadratury (4.9) liczymy przy pomocy wzoru napisanego w punkcie (1.13), dla n=1 możemy napisać wzór poniżej, do którego będziemy wykorzystali metodę całkowania przez części przy liczeniu błędu kwadratury (4.13) i twierdzenie o wartości średniej przy wykorzystaniu definicji odstępu h=b-a:

$E(f)={{1} \over {2}}\int _{a}^{b}(x-a)(x-b)f^{(2)}({\tilde {\xi }})dx={{1} \over {2}}f^{(2)}(\xi )\int _{a}^{b}(x-a)(x-b)dx=\;$
$={{1} \over {2}}f^{(2)}(\xi )\left[{{1} \over {2}}(x-a)^{2}(x-b)\right]_{a}^{b}-{{1} \over {2}}\int _{a}^{b}(x-a)^{2}dx=\;$
$=-{{1} \over {4}}{{1} \over {3}}f^{(2)}(\xi )\left[(x-a)^{3}\right]_{a}^{b}=-{{1} \over {12}}h^{3}f^{(2)}(\xi )\;$

(4.14)

Zajmiemy się teraz wielomianami stopnia drugiego n=2, wtedy możemy zapisać (4.10) licząc współczynniki A_k po kolei:

$A_{0}=h{{1} \over {2!}}\int _{0}^{2}{{t(t-1)(t-2)} \over {t}}dt={{1} \over {2}}h\int _{0}^{2}(t-1)(t-2)dt={{1} \over {2}}h\int _{0}^{2}(t^{2}-2t-t+2)dt=\;$
$={{1} \over {2}}h\left({{1} \over {3}}t^{3}-{{3} \over {2}}t^{2}+2t\right)_{0}^{2}={{1} \over {2}}h\left({{8} \over {3}}-6+4\right)=h\left({{4} \over {3}}-1\right)={{1} \over {3}}h\;$

(4.15)

$A_{1}=-h\int _{0}^{2}{{t(t-1)(t-2)} \over {t-1}}dt=-h\int _{0}^{2}t(t-2)dt=-h\left(\left[{{1} \over {2}}t^{2}(t-2)\right]_{0}^{2}-\int _{0}^{2}{{1} \over {2}}t^{2}dt\right)=\;$
$=h{{1} \over {2\cdot 3}}8={{4} \over {3}}h\;$

(4.16)

$A_{2}=h{{1} \over {2}}\int _{0}^{2}{{t(t-1)(t-2)} \over {t-2}}dt=h{{1} \over {2}}\int _{0}^{2}t(t-1)dt={{1} \over {2}}h\int _{0}^{2}(t^{2}-t)={{1} \over {2}}h\left({{1} \over {3}}t^{3}-{{1} \over {2}}t^{2}\right)_{0}^{2}=\;$
$={{1} \over {2}}h\left({{8} \over {3}}-2\right)=h\left({{4} \over {3}}-1\right)={{1} \over {3}}h\;$

(4.17)

Kwadratura (4.7) przy obliczeniach na A₀ (4.14), A₁ (4.15) i A₂ (4.16) można zapisać wzorem Simpsona przedstawioną:

$S(f)={{1} \over {3}}h(f_{0}+4f_{1}+f_{2})\;$

(4.18)

Obierzmy sobie teraz punkty a, (a+b)/2,b, x₄, takie, że x₄→(a+b)/2, wtedy okazuje się słuszny wzór Simsona (4.16) dla wielomianów trzeciego stopnia, wtedy błąd kwadratury po zastosowaniu twierdzenia o wartości średniej przy założeniu $_{h={{a+b} \over {2}}}\;$ rysujemy:

$E(f)={{1} \over {4!}}\int _{a}^{b}(x-a)\left(x-{{a+b} \over {2}}\right)^{2}(x-b)f^{(4)}({\tilde {\xi }})dx=\;$
$={{1} \over {24}}f^{(4)}(\xi )\left[{{1} \over {3}}\left(x-{{a+b} \over {2}}\right)^{3}(x-a)(x-b){\Bigg |}_{a}^{b}-{{1} \over {3}}\int _{a}^{b}\left(x-{{a+b} \over {2}}\right)^{3}(x-a+x-b)dx\right]=\;$
$=-{{1} \over {36}}f^{(4)}(\xi )\int _{a}^{b}\left(x-{{a+b} \over {2}}\right)^{4}dx=-{{1} \over {90}}f^{(4)}(\xi )h^{5}\;$

(4.19)

Złożone metody całkowania wielomianu interpolującego Newtona-Cotesa

Do całkowania zwykle nie używa się w praktyce wielomianów wyższych niż stopnia dwa, dlatego kwadratury używa się w taki sposób, że przedział (a,b) dzieli się na m równych części, dla każdej z tych przedziałów stosuje sie metodę trapezów (4.13) lub metodę Simpsona (4.18). Weźmy sobie kwadratrurę:

$E(f)=C_{r}f^{(r)}(\xi )\;$

(4.20)

to przy zmniejszaniu przedziału na p równych cząści (a-b) C_r zmniejszy się p^r razy, a błąd zmniejszy się p^r-1 razy. W metodze trapezów mamy h=(b-a)/m, która jest długością każdego z tychże przedziałów, w ten sposób piszemy kwadraturę przy zastosowaniu wzorów trapezów:

$S(f)=\sum _{k=0}^{m-1}{{1} \over {2}}h(f_{k}+f_{k+1})=h\left({{1} \over {2}}f_{0}+f_{1}+..+f_{m-1}+{{1} \over {2}}f_{m}\right){\mbox{ gdzie: }}f_{k}=f(a+kh)\;$

(4.21)

Weźmy funkcję f określonej na przedziale <a,b> określonej na zbiorze liczb zespolonych, wtedy błąd popełniony przy liczeniu metodą trapezów po podziale naszego przedziału na równe cząści piszemy wedle wzoru:

$E(f)=-{{h^{3}} \over {12}}\sum _{k=0}^{m-1}f^{(2)}(\xi _{k})=-{{(b-a)^{3}} \over {12m^{2}}}{{1} \over {m}}\sum _{k=0}^{m-1}f^{(2)}(\xi _{k})=-{{(b-a)^{3}} \over {12m^{2}}}f^{(2)}(\xi )\;$

(4.22)

Powyżej zastosowaliśmy, że liczby ξ_k leżą w przedziale <a,b>, ale oczywiste jest, że wartość średnia jest powyżej minimum funkcji i poniżej maksimum funkcji, więc ξ mieści się pomiędzy nimi, zatem ξ∈<a,b>.

Zastosujmy teraz metodę Simsona przy podziale na m/2 części przedział <a,b>, wtedy długość każdego przedziału jest h liczoną wedle h=(b-a)/m, wtedy każde z tych przedziałów ma się jako <a,a+2h>,...,<a+m-2)h,b>, wtedy każde z tych przedziałów ma długość 2h, co po zastosowaniu wzoru (4.18) dla każdego z tychże przedziałów możemy powiedzieć:

$S(f)={{h} \over {3}}\sum _{k=1}^{m/2}(f_{2k-2}+4f_{2k-1}+f_{2k})=\;$
$={{h} \over {3}}\left[f_{0}+f_{m}+2(f_{2}+f_{4}+..+f_{m-2})+4(f_{1}+f_{3}+...+f_{m-1})\right]\;$

(4.23)

Błąd złożony z parabol, możemy przepisać jako suma kwadrapul złożonych z m/2 części, których każda ma długość 2h, piszemy:

$E(f)=-{{h^{5}} \over {90m^{5}}}\sum _{k=1}^{m/2}f^{(4)}(\xi _{k})={{(b-a)^{5}} \over {180m^{4}}}{{2} \over {m}}\sum _{k=0}^{m/2}f^{(4)}(\xi _{k})={{(b-a)^{5}} \over {180m^{4}}}f^{(4)}(\xi )\;$

(4.24)

Metoda całkowania wielomianu interpolacyjnego Romberga

Niech mamy przedział <a,b> i podzielmy go na 2ⁱ równych części, wtedy długość każdego z przedziału, a także wzór na węzeł zależnego od "i" i "k" i wzór na wartość funkcji przedstawiamy je w jednej linijce jako:

$h_{i}={{b-a} \over {2^{i}}}\;$

(4.25)

$x_{i,k}=a+kh_{i}\;$

(4.26)

$f_{i,k}=f(x_{i,k})\;$

(4.27)

Wykorzystamy tutaj kwadraturę złożoną Newtona-Cotesa (4.21), który przepiszemy w postaci zależnego od "i" i od punktów "a" i "b", w ten sposób:

$T_{0,i}=h_{i}\left(\sum _{k=0}^{2^{i}}f_{i,k}-{{1} \over {2}}\left(f(a)+f(b)\right)\right)\;$

(4.28)

Napiszmy teraz błąd kwadratury (4.28) zakładając, że wzór na ten błąd jest funkcją parzystą, gdy rozszerzyć go na przedział nieskończony, i pamiętając, że dla h_i→0 bład kwadratury powninien dążyć do zera:

$I(f)-T_{0,i}=\sum _{j=1}^{\infty }C_{0,j}h_{i}^{2j}\;$

(4.29)

Określmy teraz T_1,i w zależności od T_0,i i T_1,i+1 w sposób zależny od wskaźnika i-tego:

$T_{1,i}=T_{0,j+1}+{{T_{0,j+1}-T_{0,j}} \over {2^{2}-1}}\;$

(4.30)

Mając wzór na T_1,i (4.30) a także błąd bezwzględny (4.29) na kwadraturę (4.28), to błąd T_1,i względem I(f), wykorzystując, że C_1,1=0, możemy zapisać:

$I(f)-T_{1,i}=I(f)-T_{0,j+1}+{{T_{0,j+1}-T_{0,j}} \over {2^{2}-1}}=I(f)-{{3T(0,j+1)+T_{0,i+1}-T_{0,j}} \over {3}}=\;$
$=I(f)-{{4T_{0,j+1}-T_{0,j}} \over {3}}={{1} \over {3}}\left(3I(f)-4T_{0,j+1}+T_{0,j}\right)=\;$
$={{1} \over {3}}\left[4\left(I(f)-T_{0,j+1}\right)-\left(I(f)-T_{0,j}\right)\right]={{1} \over {3}}\sum _{j=0}^{\infty }C_{0,j}\left(2^{2}h_{i+1}^{2j}-h_{i}^{2j}\right)=\;$
$={{1} \over {3}}\sum _{j=1}^{\infty }C_{0,j}h_{i}^{2j}\left({{2^{2}h_{i+1}^{2j}} \over {h_{i}^{2j}}}-1\right)=\sum _{j=1}^{\infty }C_{0,j}h_{i}^{2j}\underbrace {{{1} \over {3}}\left({{2^{2}} \over {2^{2j}}}-1\right)} _{C_{1,j}}=\sum _{j=2}^{\infty }C_{1,j}h_{i}^{2j}$

(4.31)

Główny czynnik w wyrażeniu (4.31) jest 2⁴. Napiszemy teraz wzór na T_1,i, i dowiemy się, że jest on wzorem złożony z parabol określonym w przedziale <a,b> podzielonego na 2ⁱ równych części i określonych dla m=2ⁱ⁺¹ wzorem (4.23):

$T_{1,i}={{1} \over {3}}\left(4T_{0,i+1}-T_{0,i}\right)=\;$
$={{1} \over {3}}\left[4h_{i+1}\left(\sum _{k=0}^{2^{i+1}}f_{i+1,k}-{{1} \over {2}}\left(f(a)+f(b)\right)\right)-h_{i}\left(\sum _{k=0}^{2^{i}}f_{i,k}-{{1} \over {2}}\left(f(a)+f(b)\right)\right)\right]=$
$={{1} \over {3}}h_{i+1}\left[-(f(a)+f(b))+4\left(\sum _{k=0}^{2^{i}}f_{i+1,2k}+\sum _{k=1}^{2^{i}}f_{i+1,2k-1}\right)-2\sum _{k=0}^{2^{i}}f_{i+1,2k}\right]=\;$
$={{1} \over {3}}h_{i+1}\left[-(f(a)+f(b))+4\sum _{k=1}^{2^{i}}f_{i+1,2k-1}+2\sum _{k=0}^{2^{i}}f_{i+1,2k}\right]=\;$
$={{1} \over {3}}h_{i+1}\left(f(a)+f(b)+2\sum _{k=1}^{2^{i}-1}f_{i+1,2k}+4\sum _{k=1}^{2^{i}}f_{i+1,2k-1}\right)$

(4.32)

Uogólnimy wzór (4.30) dla dowolnego "m" i "i", w ten sposób:

$T_{m,i}=T_{m-1,i+1}+{{T_{m-1,i+1}-T_{m-1,i}} \over {2^{2m}-1}}\;$

(4.33)

wtedy napiszmy błąd kwadratury T_2,i znając błąd kwadratury (4.31), w ten sposób:

$I(f)-T_{2,i}=I(f)=I(f)-\left(T_{1,i+1}+{{T_{1,i+1}-T_{1,i}} \over {2^{4}-1}}\right)={{1} \over {15}}\left(15I(f)-16T_{1,i+1}+T_{1,i}\right)=\;$
$={{1} \over {15}}\left[16\left(I(f)-T_{1,i+1}\right)-\left(I(f)-T_{1,i}\right)\right]={{1} \over {15}}\left[2^{4}\sum _{j=2}^{\infty }C_{1,j}h_{i+1}^{2j}-\sum _{j=2}^{\infty }C_{1,j}h_{i}^{2j}\right]=\;$
$={{1} \over {15}}\left[\sum _{j=2}^{\infty }h_{i}^{2j}\left({{2^{4}h_{i+1}^{2j}} \over {h_{i}^{2j}}}-1\right)\right]=\sum _{j=2}^{\infty }h_{i}^{2j}\underbrace {{{1} \over {15}}\left({{2^{4}} \over {2^{2j}}}-1\right)} _{C_{2,j}}=\sum _{j=3}^{\infty }C_{2,i}h_{i}^{2j}$

(4.34)

Główny czynnik występujący w (4.34) jest h⁶. Widzimy, że we wzorach (4.28) i (4.32) czynniki stojące przy funkcjach f(a+kh_I) są współczynnikami dodatnimi, zatem napiszmy kwadraturę T_m,i w postaci wzoru:

$T_{m,i}=h\sum _{j=0}^{2^{m+i}}d_{m,j}f(a+jh_{i})\;$

(4.35)

i udowodnimy, że dla (4.35) współczynniki d_m,j są współczynnikami dodatnimi metodą indukcji matematycznej. Dla m=0,1 zostało to udowodnione, więc udowodnimy to, że jeśli dla m współczynniki d_m,j są dodatnie to dla m+1 współczynniki d_m+1,j są też współczynnikami dodatnimi, zatem przejdźmy do sedna dowodu:

$T_{m+1,i}=T_{m,i+1}+{{T_{m,i+1}-T_{m,i}} \over {2^{2m}-1}}={{1} \over {2^{2m}-1}}\left(2^{2m}T_{m,i+1}-T_{m,i}\right)=\;$
$={{1} \over {2^{2m}-1}}\left[2^{2m}\sum _{j=0}^{2^{m+i+1}}d_{m,j}f(a+jh_{i+1})-\sum _{j=0}^{2^{m+i}}d_{m,j}f(a+jh_{i})\right]=\;$
$={{1} \over {2^{2m}-1}}{\Bigg [}2^{2m}\left(\sum _{j=0}^{2^{m+i}}d_{m,2j}f(a+2jh_{i+1})+\sum _{j=1}^{2^{m+i}}d_{m,2j-1}f(a+(2j-1)h_{i+1})\right)-\;$
$-\sum _{j=0}^{2^{m+i}}d_{m,2j}f(a+2jh_{i+1}){\Bigg ]}=\;$
$={{1} \over {2^{2m}-1}}\left[(2^{2m}-1)\sum _{j=0}^{2^{m+i}}d_{m,2j}f(a+2jh_{i+1})+2^{2m}\sum _{j=1}^{2^{m+i}}d_{m,2j-1}f(a+(2j-1)h_{i+1})\right]=\;$
$=\sum _{j=0}^{2^{m+i}}\underbrace {d_{m,2j}} _{d_{m+1,2j}}f(a+2jh_{i+1})+\sum _{j=1}^{2^{m+i}}\underbrace {{2^{2m}d_{m,2j-1}} \over {2^{2m}-1}} _{d_{m+1,2j-1}}f(a+(2j-1)h_{i+1})$

(4.36)

Zatem na podstawie dowodu (4.36) wszystkie współczynniki d_m,j są współczynnikami nieujemnymi. Błąd kwadratury (4.35) możemy przedstawić w zależności od "m" i "h" w postaci:

$I(f)-T_{m,i}=Ch_{m+i}^{2m+2}f^{(2m+2)}(\xi )\;$

(4.37)

Dla m=0 kwadratura (4.35) przechodzi w kwadraturę (4.28) i wzór (4.37) dla tej kwadratury przechodzi w (4.22) dla zdefiniowanego h_i=(b-a)/2ⁱ dla m=2ⁱ, zatem dla tego m powyższy wzór jest poprawny. Udowodnijmy metodą indukcji matematycznej twierdzenie (4.37), gdy ona jest spełnione dla m, to powinno być spełnione dla m+1 przy wykorzystaniu wzoru (4.33) i definicji $_{h_{m+i}={{b-a} \over {2^{m+i}}}}\;$, wtedy:

$I(f)-T_{m+1,i}=I(f)-T_{m,i+1}-{{T_{m,i+1}-T_{m,i}} \over {2^{2(m+1)}-1}}=\;$
$={{1} \over {2^{2m+2}-1}}\left[2^{2m+2}\left(I(f)-T_{m,i+1}\right)-(I(f)-T_{m,i})\right]=\;$
$={{1} \over {2^{2m+2}-1}}\left[2^{2m+2}Ch_{m+(i+1)}^{2m+2}f^{(2m+2)}(\xi _{2})-C2^{2m+2}h_{(m+1)+i}^{2m+2}f^{2m+2}(\xi _{1})\right]=\;$
$={{2^{2m+2}} \over {2^{2m+2}-1}}h_{(m+1)+i}^{2m+2}C\left[f^{(2m+2)}(\xi _{2})-f^{(2m+2)}(\xi _{2})\right]=C_{2}h_{(m+1)+i}^{2m+2}h_{(m+1)+i}^{2}f^{2(m+1)+2}(\xi )=\;$
$=C_{2}h_{(m+1)+i}^{2(m+1)+2}f^{(2(m+1)+2)}(\xi )$

(4.38)

W powyższym dowodzie przyjęto, że pochodna $f^{(2m+3)}(\xi )={{f^{(2m+2)}(\xi _{2})-f^{(2m+2)}(\xi _{1})} \over {(\xi _{2}-\xi _{1})}}\;$ jest równa zero, bo wartości pochodnych f^(2m+2)(ξ) w punktach ξ_i i ξ₂ są sobie równe, bo błąd kwadratury (4.35) jest zależny od głównego czynnika $h_{i}^{2(m+1)+2}\;$, zatem na podstawie tego i twierdzenia Lagrange'a pochodna f^(2(m+1)+1) jest równa zero, zatem tylko pochodna f^(2(m+1)+2)(ξ) jest nie równa zero w tym naszym przypadku. Co kończy dowód twierdzenia (4.37) na podstawie dowodu (4.38).

Całkowanie wielomianu interpolacyjnego Gaussa

Dotychczas rozpatrywaliśmy wielomiany interpolacyjnego, których kwadratury są określone przez wzór (4.7) przy definicji A_k (4.4). Były to kwadratury rzędu co najmniej N+1. Rozpatrzmy przy liczbie N, aby rząd kwadratury był najwyższy i wynosił 2(N+1). Taką kwadraturę nazywamy kwadraturą Gaussa, wtedy mamy do wyznaczenia 2(N+1) stałych do wyznaczenia. Tak więc problem sprowadza się do odpowiedniego wyboru węzłów w kwadraturze.

Rozpatrzmy teraz wielomiany ortogonalne, czyli rozpatrzmy ciąg wielomianów {P_n(x)}=P₀(x), P₁(x),...,P_n(x), w którym n oznacza najwyższy numer wielomanu ortogonalnego. Wielomiany ortogonalne spełniają warunek ortogonalności przy wadze p(x), który piszemy:

$(P_{r},P_{s})=\int _{a}^{b}p(x)P_{r}(x)P_{s}(x)dx{\mbox{ dla }}r\neq s\;$

(4.39)

Rozpatrzmy twierdzenie mówiące o liczbie pierwiastków wielomianów ortogonalnych:

Twierdzenie
Wielomiany ortogonalne P_s(x) mają jednokrotne pierwiastki rzeczywiste należące do przedziału (a,b)

Dowód
Niech pierwiastkami wielomianu P_s(x) są pierwiastkami nieparzystej niejednokrotności, i oznaczmy te pierwiastki symbolami ξ₁, ξ₂,..,ξ_n, i niech m<k, wtedy możemy zbudować wielomian W(x)=(x-ξ₁)...(x-ξ_n), który możemy rozłożyć w bazie funkcji ortogonalnych P₀(x),.., P_m(x) wedle sposobu:

$W(x)=a_{0}P_{0}(x)+...+a_{m}P_{m}(x)\;$

(4.40)

Wtedy możemy napisać iloczyn skalarny, który z definicji wielomianów ortogonalnych (4.39) jest równy zero:

$(W,P_{k})=\sum _{i=0}^{m}a_{i}(P_{i},P_{k})=0\;$

(4.41)

Z drugiej jednak strony wyrazenie (4.41) jest większe od zera, bo P_k posiada pierwiastki o nieparzystej krotności, także i W(x), zatem wyrażenie W(x)P_k(x) posiada potęgi parzyste przy (x-ξ_i), wtedy powinno być (W,P_k)>0, zatem dochodzimy do sprzeczności, stąd wynika, że wielomian P_k(x) posiada k pierwiastków jednokrotnych.

Twierdzenie
 

• 1) Nie istnieje wielomian rzędu wyższego niż 2(N+1).
• 2) Kwadratura określona wzorem (4.7) poprzez współczynniki A_k jest rzędu 2(N+1), gdy x_k są pierwiastkami wielomianu P_N+1 dla ciągu wielomianów ortogonalnych {P_n(x)}.

Dowód
1)Weźmy sobie wielomian $W(x)=\omega _{N+1}^{2}(x)\;$, który jest stopnia 2(N+1), który jest zdefiniowany o wielomian ω_N+1=(x-₀)..(x-x_N), wtedy możemy powiedzieć:

$I(W)=\int _{a}^{b}p(x)W(x)dx>0\;$

(4.42)

$S(W)=\sum _{k=0}^{N}A_{k}W(x_{k})=0\;$

(4.43)

(4.43) jest równe zero dla ω_n+1 liczona w punktach x_k, które są jego pierwiastkami tego wielomianu, dla którego kwadratura jest równa zero. Zatem mamy wielomian rzędu 2(N+1), dla których kwadratura (4.7) jest niedokładna, stąd z definicji rzędu kwadratury jest ona rzędu 2(N+1).

2)Weźmy sobie wielomian W(x) stopnia niższego niż N+1, dalej rozłóżmy wielomian ω_N+1 w bazie funkcji ortogonalnych {P_n(x)} dla a_N+1≠0, wtedy mamy:

$\omega _{N+1}(x)=a_{0}P_{0}(x)+...+a_{N+1}P_{N+1}(x)\;$

(4.44)

A wielomian W(x) możemy rozłożyć w bazie funkcji ortogonalnych dla wielomianów co najwyżej rzędu N+1:

$W(x)=a_{0}P_{0}(x)+..+a_{N}P_{N}(x)\;$

(4.45)

Wtedy możemy napisać kwadraturę w postaci wzoru, dla której kwadratura jest dokładnie równa zero w punktach pierwiastków ω_N+1:

$\int _{a}^{b}p(x)W(x)\omega _{N+1}dx=\sum _{k=0}^{N}A_{k}R(x_{k})\;$

(4.46)

Wtedy iloczyn skalarny wielomianu W(x) (4.45) z wielomianem ω_N+1 (4.44) w punktach pierwiastków rzeczywistych jednokrotnych wielomianu P_k(x) jest równy zera, stąd mamy:

$(W,\omega _{N+1})=(a_{0}P_{0}(x)+..+a_{N}P_{N}(x),a_{0}P_{0}(x)+...+a_{N+1}P_{N+1}(x))=\;$
$=a_{0}(P_{0},P_{0})+...+a_{N}(P_{N},P_{N})=0\;$

(4.47)

Aby powyższe wyrażenie było zawsze równe zero, to współczynniki a_k powinny być równe zero dla i=0,1,..,N, stąd wniosek, że funkcja (4.44) powinna być równa:

$\omega _{N+1}(x)=a_{N+1}P_{N+1}(x)\;$

(4.48)

Twierdzenie
Współczynniki A_k w kwadraturach są zawsze dodatnie.

Dowód
Weźmy sobie wielomian stopnia 2N, który zawsze przyjmuje wartości nieujemne, a jego definicja jest dla i=0,1,..,N:

$R(x)=\left[(x-x_{0})^{2}...(x-x_{i-1})(x-x_{i+1})...(x-x_{N})\right]^{2}\;$

(4.49)

Na podstawie definicji kwadratury i nieujemności (4.49) współczynniki A_k są zawsze dodatnie, co pokazuje dowód:

$0<\int _{a}^{b}p(x)R_{i}(x)dx=\sum _{k=0}^{N}A_{k}R_{i}(x_{k})=A_{i}R_{i}(x_{i})\;$

(4.50)

Ponieważ R_i(x_i) jest zawsze nieujemne na mocy (4.49), to A_i wedle (4.50) jest dodatnie.

Weźmy sobie punkty $x_{0}\;$, $x_{1}\;$,..,$x_{N}\;$, ${\tilde {x}}_{0}\;$, ${\tilde {x}}_{1}\;$,...,${\tilde {x}}_{N}\;$, takie że: $x_{k}\neq {\tilde {x}}_{k}\;$, ale by zachodziło $x_{k}\rightarrow {\tilde {x}}_{k}\;$, wtedy błąd interpolacyjny funkcji interpolującej względem funkcji interpolowanej wyrażamy przy pomocy wzoru (1.24), zatem piszemy ją:

$f(x)-W_{2N+1}={{1} \over {(2N+2)!}}f^{(2N+2)}({\tilde {\xi }})\omega _{N+1}^{2}(x)\;$

(4.51)

Błąd kwadratury (4.7), przy pomocy wagi p(x), piszemy podobnym wzorem do (4.8) dla wielomianu $\omega _{N+1}^{2}(x)\;$, który jest rzędu 2N+2:

$E(f)={{1} \over {(2N+2)!}}f^{(2N+2)}(\xi )\int _{a}^{b}p(x)\omega _{N+1}^{2}(x)dx\;$

(4.52)

Wprowadzenie do kwadratur określonych w przedziale skończonym, a kwadratura Gaussa-Legendre'a

Będziemy tutaj określać kwadraturę przy wadze p(x)=1, którą nazywamy kwadraturą Gaussa-Legendre'a. Określamy teraz ciąg wielomianów Legendre'a jako:

$P_{n}(x)={{1} \over {2^{n}n!}}{{d^{n}} \over {dx^{n}}}(x^{2}-1)^{n}\;$

(4.53)

wtedy współczynniki A_n (4.4) i błąd kwadratury Gaussa-Legendre'a określać będziemy wzorami dla -1<ξ<1:

$A_{k}={{2} \over {(N+1)P_{N+2}(x_{k})P_{N+1}^{'}(x_{k})}}\;$

(4.54)

$E(f)={{2^{2N+3}((N+1)!)^{4}} \over {(2N+3)((2N+2)!)^{3}}}f^{(2N+2)}(\xi )\;$

(4.55)

Określmy teraz nową zmienną t, którą określamy względem przedziału <a,b>, którą określamy względem punktów "a" i "b" i zmiennej x określonej na przedziale <-1,1>:

$t={{a+b} \over {2}}+{{b-a} \over {2}}x\;$

(4.56)

wtedy korzystając z definicji ω_N+1, a także definicji całki na funkcję błędu E(f) określonych na przedziale <-1,1>, możemy poszerzyć funkcję błędu E(f) (4.55) na przedziale <a,b> wychodząc od wzoru (4.55), co jest błędem kwadratury dla tego przedziału, wykorzystując przy tym fakt $t_{k}={{a+b} \over {2}}+{{b-a} \over {2}}x_{k}$ wynikającą z (4.56):

$E(f)={{1} \over {(2N+2)!}}f^{(2N+2)}(\xi )\int _{a}^{b}\omega _{N+1}^{2}(x)dx={{2^{2N+3}((N+1)!)^{4}} \over {(2N+3)((2N+2)!)^{3}}}f^{(2N+2)}(\xi )=\;$
$={{1} \over {(2N+2)!}}f^{(2N+2)}(\xi )\int _{-1}^{1}\left[\left({{b-a} \over {2}}\right)^{N+1}\omega _{N+1}(t)\right]^{2}d{{b-a} \over {2}}dt=\;$
$=\left({{b-a} \over {2}}\right)^{2N+3}{{1} \over {(2N+2)!}}f^{(2N+2)}(\xi )\int _{-1}^{1}\omega _{N+1}^{2}(x)dx=\;$
$=\left({{b-a} \over {2}}\right)^{2N+3}{{2^{2N+3}((N+1)!)^{4}} \over {(2N+3)((2N+2)!)^{3}}}f^{(2N+2)}(\xi )=\;$
$={{(b-a)^{2N+3}((N+1)!)^{4}} \over {(2N+3)((2N+2)!)^{3}}}f^{(2N+2)}(\xi ){\mbox{ dla }}a<\xi <b\;$

(4.57)

Całkowanie interpolacyjne dla przedziału skończonego, a kwadratura Gaussa-Jacobiego

Rozpatrzmy teraz kwadraturę Gauussa-Jacobiego przy normie z definiowanej p(x)=(1-x)^α(1+x)^β, w którem parametry α i β są większe niż minus jeden, wtedy możemy zdefiniować funkcję (wielomian) Gaussa-Jacobiego, które są wielomianami ortogonalnymi z powyższą wagą:

$J_{n}(x,\alpha ,\beta )={{(-1)^{n}} \over {2^{n}n!}}(1-x)^{-\alpha }(1+x)^{-\beta }{{d^{n}} \over {dx^{n}}}\left[(1-x)^{\alpha +n}(1+x)^{\beta +n}\right]\;$

(4.58)

wtedy możemy zdefiniować współczynniki A_k (4.5) kwadratury Gaussa-Jacobiego zależnego od N, α,β:

$A_{k}=-{{(2N+\alpha +\beta +4)\Gamma (N+\alpha +2)\Gamma (N+\beta +2)\cdot 2^{\alpha +\beta }} \over {(N+\alpha +\beta +2)\Gamma (N+\alpha +\beta +2)J_{N+1}(x_{k},\alpha ,\beta )J_{N+2}(x_{k},\alpha ,\beta )(N+2)}}\;$

(4.59)

Funkcja błędu dla kwadratury Gaussa-Jacobiego przy normie p(x) wyżej podanej przedstawiamy w zależności od funkcji interpolowanej f(x):

$E(f)={{\Gamma (N+\alpha +2)\Gamma (N+\beta +2)\Gamma (N+\alpha +\beta +2)(N+1)!2^{2N+\alpha +\beta +3}} \over {(2N+\alpha +\beta +3)[\Gamma (2N+\alpha +\beta +3)]^{2}(2N+2)!}}f^{(2N+2)}(\xi )\;$

(4.60)

Powyższe wzory są słuszne dla α,β>-1, ale dla α=β=0 przechodzą one we wzory Gaussa-Legendre'a, bo wtedy norma p(x) jest równa jedności, przyjmują zaś one konkretne wartości uwzględniające osobliwości normy p(x) podanej powyżej.

Całkowanie interpolacyjne dla przedziału skończonego, a kwadratura Gaussa-Czybyszewa

W powyższych wzorach przyjmijmy α=β=-1/2, wtedy powyżej zdefiniowana norma przechodzi w p(x)=(1-x²)^-1/2. Kwadratury z tak określoną wagą nazywamy kwadraturami Gaussa-Czebyszewa. Ciąg wielomianów ortogonalnych określamy przez (1.27), wtedy węzły x_k i współczynniki A_k określamy przez:

$x_{k}=\cos {{(2k+1)\pi } \over {2N+2}}\;$

(4.61)

$A_{k}={{\pi } \over {N+1}}\;$

(4.62)

Całkowanie interpolacyjne dla przedziału nieskończonego

Omówimy tutaj zastosowanie całek Gaussa dla przedziału nieskończonego, czyli jednostronnie nieskończonego lub obustronnie. Rozpatrzmy tutaj normę p(x)=e^-x, wtedy możemy zdefiniować wzory na wielomiany Laguerre'a, które określamy przy pomocy współczynnika n dla pierwszego rozważanego przypadku, czyli dla przedziału od zera do nieskończoności:

$L_{n}(x)=(-1)^{n}e^{x}{{d^{n}} \over {dx^{n}}}\left(x^{n}e^{-x}\right)\;$

(4.63)

wtedy powyższa waga przy zbieżności całki $\int _{0}^{\infty }p(x)W(x)dx\;$, gdzie W(x) jest wielomianem dowolnego stopnia, wynika że całka określająca współczynniki A_k jest zbieżna, wtedy możemy określić współczynniki A_k i funkcję błędu zależnego od N i funkcji f:

$A_{k}={{((N+1)!)^{2}} \over {L_{N+1}(x_{k})L_{N+1}(x_{k})}}\;$

(4.64)

$E(f)={{(N+1)!)^{2}} \over {(2N+2)!}}f^{(2N+2)}(\eta )\;$

(4.65)

Omówimy tutaj kwadraturę dla przedziału obustronnie nieskończonego przy normie p(x)=e^-x2, dla której zdefiniujemy ciąg wielomianów ortogonalnych zwanych wzorami Hermite'a:

$H_{n}(x)=(-1)^{n}e^{x^{2}}{{d^{n}} \over {dx^{n}}}e^{-x^{2}}\;$

(4.66)

i przy których zdefiniujemy współczynniki A_k i błąd dla kwadratury określonych przy pomocy wielomianów Hermite'a:

$A_{k}={{2^{N+2}(N+1)!} \over {H_{N+1}^{'}(x_{k})H_{N+2}(x_{k})}}\;$

(4.67)

$E(f)={{(N+1)!{\sqrt {\pi }}} \over {2^{n+1}(2n+2)!}}f^{2N+2}(\eta ){\mbox{ dla }}\eta \in (-\infty ,\infty )\;$

(4.68)

Złożone metody całkowania interpolacyjnego Gaussa

Będziemy tutaj stosowali kwadratury Gaussa-Legendre'a, którego błąd kwadratury dla przedziału $h={{b-a} \over {n}}\;$, który określamy na podstawie wzoru (4.57), pisaną:

$E(f)=n\left({{b-a} \over {n}}\right)^{2N+3}{{((N+1)!)^{4}} \over {(2N+3)((2N+2)!)^{3}}}{{1} \over {n}}\sum _{k=1}^{n}f^{(2N+2)}(\xi _{k})=\;$
$=n\left({{b-a} \over {n}}\right)^{2N+3}{{((N+1)!)^{4}} \over {(2N+3)((2N+2)!)^{3}}}f^{(2N+2)}(\xi )\;$

(4.69)

We wzorze (4.68) przyjmijmy, że N=1 i n=m/2, gdy m parzyste, wtedy całkę funkcji f(x) w przedziale <a,b> po podziale na n części określamy:

$\int _{a}^{b}f(x)dx=\sum _{j=0}^{{{m} \over {2}}-1}\int _{a+2hj}^{a+2h(j+1)}f(x)dx\;$

(4.70)

Wyznaczmy teraz pierwiastki wielomianu P_N+1(x), który określamy poprzez wzór (4.53), wtedy jego pierwiastki wynoszą według obliczeń:

$0=P_{2}(x)={{1} \over {4\cdot 2}}{{d^{2}} \over {dx^{2}}}(x^{2}-1)^{2}={{1} \over {4\cdot 2}}{{d^{2}} \over {dx^{2}}}(x^{4}-2x^{2}+1)={{1} \over {8}}{{d} \over {dx}}\left(4x^{3}-4x\right)=\;$
$={{1} \over {2}}\left(3x^{2}-1\right)\Rightarrow x=\pm {{1} \over {\sqrt {3}}}\;$

(4.71)

Wyznaczmy czemu odpowiadają węzły w przedziale <a+2hj,a+2h(j+1)> na podstawie (4.56) wiedząc, że węzły w przedziale <-1,1> są określone przez (4.71):

$t_{k}={{a+2h(j+1)+a+2hj} \over {2}}+{{a+2h(j+1)-a-2hj} \over {2}}x_{k}={{2a+4hj+2h} \over {2}}+{{2h} \over {2}}x_{k}=\;$
$=a+2hj+h+hx_{k}=a+h(1+x_{k}+2j)\;$

(4.72)

Wtedy kwadraturę złozoną Gaussa, gdy węzły w przedziale wyżej określone są według wzoru (4.72), piszemy na podstawie (4.70):

$S(f)=h\sum _{j=0}^{{{m} \over {2}}-1}\left[f\left(a+h\left(1-1/{\sqrt {3}}+2j\right)\right)+f\left(a+h\left(1+1/{\sqrt {3}}+2j\right)\right)\right]\;$

(4.73)

Wyznaczmy teraz błąd złożonej kwadratury Gaussa wykorzystując przy tym wzór (4.69) dla N=1 i n=m/2, wtedy:

$E(f)=m/2\left({{b-a} \over {m/2}}\right)^{5}{{2^{4}} \over {5\cdot 24^{3}}}f^{(4)}(\xi )={{(b-a)^{5}} \over {270m^{4}}}f^{(4)}(\xi ){\mbox{ dla }}\xi \in (a,b)\;$

(4.74)