Metody numeryczne fizyki/Całkowanie numeryczne funkcji interpolacyjnej

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 Szablon:LinkWzór (gdzie ySzablon:Sub jest to właśnie funkcja F(xSzablon:Sub)) przy definicji ΦSzablon:Sub(x) Szablon:LinkWzór, wtedy możemy napisać przybliżoną równość pomiędzy całką dokładną, a przybliżoną, którą piszemy w postaci: Szablon:CentrujWzór Bład bezwzględny pomiędzy wartością funkcji, a jej wartością przybliżoną, możemy przepisać przy pomocy schematu: Szablon:CentrujWzór to wtedy wykorzystując Szablon:LinkWzór 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: Szablon:CentrujWzór 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: Szablon:CentrujWzór Współczynniki ASzablon:Sub określonego na przedziale <a,b> przy pomocy funkcji ΦSzablon:Sub(x) jako całka przy normie p(x) piszemy: Szablon:CentrujWzór 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: Szablon:CentrujWzór i będziemy je przybliżać w pewnym przedziale <a,b>, gdy obierzemy punkty xSzablon:Sub i będziemy je stosowali do funkcji f(x) przy stałych ASzablon:Sub podobnie zdefiniowanej jak w punkcie Szablon:LinkWzór, wtedy mamy: Szablon:CentrujWzór Punkty xSzablon:Sub nazywamy węzłami kwadratury, a wzór Szablon:LinkWzór 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 Szablon:LinkWzór i Szablon:LinkWzór przy definicji błędu wzoru funkcji interpolującej Szablon:LinkWzór, wtedy: Szablon:CentrujWzór Jako kryterium dokładności kwadratury będziemy pisali S(W) i I(W), gdy W jest pewnym wielomianem. Kwadratura Szablon:LinkWzór 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. Szablon:Twierdzenie

Całkowanie wielomianu interpolacyjnego z ustalonymi punktami

Będziemy tutaj rozpatrywali, że mamy ustalone węzły, tzn. xSzablon:Sub dla k=0,1,2,..,n. Poniżej będziemy liczyli współczynniki ASzablon:Sub, 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 Szablon:LinkWzór, w którym zdefiniujemy ASzablon:Sub sposobem Szablon:LinkWzór, a wielomian ΨSzablon:Sub(x) jest zdefiniowany przez Szablon:LinkWzór, w którym błąd kwadratury jest zdefiniowaliśmy poprzez wzór Szablon:LinkWzór. Całkowanie wielomianu ΦSzablon:Sub(x) w przedziale (a,b) jest napisane jako definicję ASzablon:Sub: Szablon:CentrujWzór Szablon:Twierdzenie Szablon:Dowód Aby policzyć różne współczynniki ASzablon:Sub należy skorzystać ze wzoru Szablon:LinkWzór, wtedy poszczególne wspomniane współczynniki dla n=1 liczymy jako: Szablon:ElastycznyWiersz Wtedy kwadraturę Szablon:LinkWzór przy pomocy współczynników Szablon:LinkWzór i Szablon:LinkWzór piszemy przy pomocy wartości funkcji fSzablon:Sub i fSzablon:Sub wyznaczonych w punktach "a" i "b": Szablon:CentrujWzór Błąd kwadratury Szablon:LinkWzór liczymy przy pomocy wzoru napisanego w punkcie Szablon:LinkWzór, 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 Szablon:LinkWzór i twierdzenie o wartości średniej przy wykorzystaniu definicji odstępu h=b-a: Szablon:CentrujWzór Zajmiemy się teraz wielomianami stopnia drugiego n=2, wtedy możemy zapisać Szablon:LinkWzór licząc współczynniki ASzablon:Sub po kolei: Szablon:CentrujWzór Szablon:CentrujWzór Szablon:CentrujWzór Kwadratura Szablon:LinkWzór przy obliczeniach na ASzablon:Sub Szablon:LinkWzór, ASzablon:Sub Szablon:LinkWzór i ASzablon:Sub Szablon:LinkWzór można zapisać wzorem Simpsona przedstawioną: Szablon:CentrujWzór Obierzmy sobie teraz punkty a, (a+b)/2,b, xSzablon:Sub, takie, że xSzablon:Sub→(a+b)/2, wtedy okazuje się słuszny wzór Simsona Szablon:LinkWzór dla wielomianów trzeciego stopnia, wtedy błąd kwadratury po zastosowaniu twierdzenia o wartości średniej przy założeniu ${}_{h=\frac{a+b}{2}}$ rysujemy: Szablon:CentrujWzór

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 Szablon:LinkWzór lub metodę Simpsona Szablon:LinkWzór. Weźmy sobie kwadratrurę: Szablon:CentrujWzór to przy zmniejszaniu przedziału na p równych cząści (a-b) CSzablon:Sub zmniejszy się pSzablon:Sup razy, a błąd zmniejszy się pSzablon:Sup 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: Szablon:CentrujWzór 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: Szablon:CentrujWzór Powyżej zastosowaliśmy, że liczby ξSzablon:Sub 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 Szablon:LinkWzór dla każdego z tychże przedziałów możemy powiedzieć: Szablon:CentrujWzór 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: Szablon:CentrujWzór

Metoda całkowania wielomianu interpolacyjnego Romberga

Niech mamy przedział <a,b> i podzielmy go na 2Szablon:Sup 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: Szablon:ElastycznyWiersz Wykorzystamy tutaj kwadraturę złożoną Newtona-Cotesa Szablon:LinkWzór, który przepiszemy w postaci zależnego od "i" i od punktów "a" i "b", w ten sposób: Szablon:CentrujWzór Napiszmy teraz błąd kwadratury Szablon:LinkWzór 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 hSzablon:Sub→0 bład kwadratury powninien dążyć do zera: Szablon:CentrujWzór Określmy teraz TSzablon:Sub w zależności od TSzablon:Sub i TSzablon:Sub w sposób zależny od wskaźnika i-tego: Szablon:CentrujWzór Mając wzór na TSzablon:Sub Szablon:LinkWzór a także błąd bezwzględny Szablon:LinkWzór na kwadraturę Szablon:LinkWzór, to błąd TSzablon:Sub względem I(f), wykorzystując, że CSzablon:Sub=0, możemy zapisać: Szablon:CentrujWzór Główny czynnik w wyrażeniu Szablon:LinkWzór jest 2Szablon:Sup. Napiszemy teraz wzór na TSzablon:Sub, i dowiemy się, że jest on wzorem złożony z parabol określonym w przedziale <a,b> podzielonego na 2Szablon:Sup równych części i określonych dla m=2Szablon:Sup wzorem Szablon:LinkWzór: Szablon:CentrujWzór Uogólnimy wzór Szablon:LinkWzór dla dowolnego "m" i "i", w ten sposób: Szablon:CentrujWzór wtedy napiszmy błąd kwadratury TSzablon:Sub znając błąd kwadratury Szablon:LinkWzór, w ten sposób: Szablon:CentrujWzór Główny czynnik występujący w Szablon:LinkWzór jest hSzablon:Sup. Widzimy, że we wzorach Szablon:LinkWzór i Szablon:LinkWzór czynniki stojące przy funkcjach f(a+khSzablon:Sub) są współczynnikami dodatnimi, zatem napiszmy kwadraturę TSzablon:Sub w postaci wzoru: Szablon:CentrujWzór i udowodnimy, że dla Szablon:LinkWzór współczynniki dSzablon:Sub 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 dSzablon:Sub są dodatnie to dla m+1 współczynniki dSzablon:Sub są też współczynnikami dodatnimi, zatem przejdźmy do sedna dowodu: Szablon:CentrujWzór Zatem na podstawie dowodu Szablon:LinkWzór wszystkie współczynniki dSzablon:Sub są współczynnikami nieujemnymi. Błąd kwadratury Szablon:LinkWzór możemy przedstawić w zależności od "m" i "h" w postaci: Szablon:CentrujWzór Dla m=0 kwadratura Szablon:LinkWzór przechodzi w kwadraturę Szablon:LinkWzór i wzór Szablon:LinkWzór dla tej kwadratury przechodzi w Szablon:LinkWzór dla zdefiniowanego hSzablon:Sub=(b-a)/2Szablon:Sup dla m=2Szablon:Sup, zatem dla tego m powyższy wzór jest poprawny. Udowodnijmy metodą indukcji matematycznej twierdzenie Szablon:LinkWzór, gdy ona jest spełnione dla m, to powinno być spełnione dla m+1 przy wykorzystaniu wzoru Szablon:LinkWzór i definicji ${}_{h_{m+i}=\frac{b-a}{2^{m+i}}}$, wtedy: Szablon:CentrujWzór W powyższym dowodzie przyjęto, że pochodna Szablon:Formuła jest równa zero, bo wartości pochodnych fSzablon:Sup(ξ) w punktach ξSzablon:Sub i ξSzablon:Sub są sobie równe, bo błąd kwadratury Szablon:LinkWzór jest zależny od głównego czynnika Szablon:Formuła, zatem na podstawie tego i twierdzenia Lagrange'a pochodna fSzablon:Sup jest równa zero, zatem tylko pochodna fSzablon:Sup(ξ) jest nie równa zero w tym naszym przypadku. Co kończy dowód twierdzenia Szablon:LinkWzór na podstawie dowodu Szablon:LinkWzór.

Całkowanie wielomianu interpolacyjnego Gaussa

Dotychczas rozpatrywaliśmy wielomiany interpolacyjnego, których kwadratury są określone przez wzór Szablon:LinkWzór przy definicji ASzablon:Sub Szablon:LinkWzór. 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 {PSzablon:Sub(x)}=PSzablon:Sub(x), PSzablon:Sub(x),...,PSzablon:Sub(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: Szablon:CentrujWzór Rozpatrzmy twierdzenie mówiące o liczbie pierwiastków wielomianów ortogonalnych: Szablon:Twierdzenie Szablon:Dowód Szablon:Twierdzenie Szablon:Dowód Szablon:Twierdzenie Szablon:Dowód Weźmy sobie punkty Szablon:Formuła, Szablon:Formuła,..,Szablon:Formuła, Szablon:Formuła, Szablon:Formuła,...,Szablon:Formuła, takie że: Szablon:Formuła, ale by zachodziło Szablon:Formuła, wtedy błąd interpolacyjny funkcji interpolującej względem funkcji interpolowanej wyrażamy przy pomocy wzoru Szablon:LinkWzór, zatem piszemy ją: Szablon:CentrujWzór Błąd kwadratury Szablon:LinkWzór, przy pomocy wagi p(x), piszemy podobnym wzorem do Szablon:LinkWzór dla wielomianu Szablon:Formuła, który jest rzędu 2N+2: Szablon:CentrujWzór

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: Szablon:CentrujWzór wtedy współczynniki ASzablon:Sub Szablon:LinkWzór i błąd kwadratury Gaussa-Legendre'a określać będziemy wzorami dla -1<ξ<1: Szablon:ElastycznyWiersz 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>: Szablon:CentrujWzór wtedy korzystając z definicji ωSzablon:Sub, 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) Szablon:LinkWzór na przedziale <a,b> wychodząc od wzoru Szablon:LinkWzór, co jest błędem kwadratury dla tego przedziału, wykorzystując przy tym fakt Szablon:Formuła wynikającą z Szablon:LinkWzór: Szablon:CentrujWzór

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)Szablon:Sup(1+x)Szablon:Sup, 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ą: Szablon:CentrujWzór wtedy możemy zdefiniować współczynniki ASzablon:Sub Szablon:LinkWzór kwadratury Gaussa-Jacobiego zależnego od N, α,β: Szablon:CentrujWzór Funkcja błędu dla kwadratury Gaussa-Jacobiego przy normie p(x) wyżej podanej przedstawiamy w zależności od funkcji interpolowanej f(x): Szablon:CentrujWzór 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-xSzablon:Sup)Szablon:Sup. Kwadratury z tak określoną wagą nazywamy kwadraturami Gaussa-Czebyszewa. Ciąg wielomianów ortogonalnych określamy przez Szablon:LinkWzór, wtedy węzły xSzablon:Sub i współczynniki ASzablon:Sub określamy przez: Szablon:ElastycznyWiersz

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)=eSzablon:Sup, 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: Szablon:CentrujWzór wtedy powyższa waga przy zbieżności całki $\underset{0}{\overset{\mathrm{\infty }}{\int }}p(x)W(x)dx$, gdzie W(x) jest wielomianem dowolnego stopnia, wynika że całka określająca współczynniki ASzablon:Sub jest zbieżna, wtedy możemy określić współczynniki ASzablon:Sub i funkcję błędu zależnego od N i funkcji f: Szablon:ElastycznyWiersz Omówimy tutaj kwadraturę dla przedziału obustronnie nieskończonego przy normie p(x)=eSzablon:Sup, dla której zdefiniujemy ciąg wielomianów ortogonalnych zwanych wzorami Hermite'a: Szablon:CentrujWzór i przy których zdefiniujemy współczynniki ASzablon:Sub i błąd dla kwadratury określonych przy pomocy wielomianów Hermite'a: Szablon:ElastycznyWiersz

Złożone metody całkowania interpolacyjnego Gaussa

Będziemy tutaj stosowali kwadratury Gaussa-Legendre'a, którego błąd kwadratury dla przedziału Szablon:Formuła, który określamy na podstawie wzoru Szablon:LinkWzór, pisaną: Szablon:CentrujWzór We wzorze Szablon:LinkWzór 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: Szablon:CentrujWzór Wyznaczmy teraz pierwiastki wielomianu PSzablon:Sub(x), który określamy poprzez wzór Szablon:LinkWzór, wtedy jego pierwiastki wynoszą według obliczeń: Szablon:CentrujWzór Wyznaczmy czemu odpowiadają węzły w przedziale <a+2hj,a+2h(j+1)> na podstawie Szablon:LinkWzór wiedząc, że węzły w przedziale <-1,1> są określone przez Szablon:LinkWzór: Szablon:CentrujWzór Wtedy kwadraturę złozoną Gaussa, gdy węzły w przedziale wyżej określone są według wzoru Szablon:LinkWzór, piszemy na podstawie Szablon:LinkWzór: Szablon:CentrujWzór Wyznaczmy teraz błąd złożonej kwadratury Gaussa wykorzystując przy tym wzór Szablon:LinkWzór dla N=1 i n=m/2, wtedy: Szablon:CentrujWzór