Metody numeryczne fizyki/Całkowanie numeryczne funkcji interpolacyjnej

Z Wikibooks, biblioteki wolnych podręczników.
Przejdź do nawigacji Przejdź do wyszukiwania
Metody numeryczne fizyki
Metody numeryczne fizyki
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 oraz miejsca pochodzenia książki i jej jakikolwiek części, a także treści, teksty, tabele, wykresy, rysunki, wzory i inne elementy oraz ich części zawarte w książce, i tą książkę, nawet w postaci przerobionej nie można umieszczać w jakikolwiek formie na czasopismach naukowych, archiwach prac, itp.
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 yk jest to właśnie funkcja F(xk)) 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:

(4.1)

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

(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:

(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:

(4.4)

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

(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[edytuj]

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:

(4.6)

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

(4.7)

Punkty xK 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:

(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 niezależna od ilości węzłów "n":
(4.9)

Całkowanie wielomianu interpolacyjnego z ustalonymi punktami[edytuj]

Będziemy tutaj rozpatrywali, że mamy ustalone węzły, tzn. xk dla k=0,1,2,..,n. Poniżej będziemy liczyli współczynniki Ak, 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[edytuj]

Będziemy się tutaj zajmowali kwadraturą Newtona-Cotesa określoną przez wzór początkowy (4.7), w którym zdefiniujemy AK 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ę Ak:

(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)=Ln(x) i I(f)=I(Ln)=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 Ak należy skorzystać ze wzoru (4.10), wtedy poszczególne wspomniane współczynniki dla n=1 liczymy jako:

(4.11)
(4.12)

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

(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:



(4.14)

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


(4.15)

(4.16)

(4.17)

Kwadratura (4.7) przy obliczeniach na A0 (4.14), A1 (4.15) i A2 (4.16) można zapisać wzorem Simpsona przedstawioną:

(4.18)

Obierzmy sobie teraz punkty a, (a+b)/2,b, x4, takie, że x4→(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 rysujemy:



(4.19)

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

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ę:

(4.20)

to przy zmniejszaniu przedziału na p równych cząści (a-b) Cr zmniejszy się pr razy, a błąd zmniejszy się pr-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:

(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:

(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ć:


(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:

(4.24)

Metoda całkowania wielomianu interpolacyjnego Romberga[edytuj]

Niech mamy przedział <a,b> i podzielmy go na 2i 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:

(4.25)
(4.26)
(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:

(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 hi→0 bład kwadratury powninien dążyć do zera:

(4.29)

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

(4.30)

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




(4.31)

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





(4.32)

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

(4.33)

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



(4.34)

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

(4.35)

i udowodnimy, że dla (4.35) współczynniki dm,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 dm,j są dodatnie to dla m+1 współczynniki dm+1,j są też współczynnikami dodatnimi, zatem przejdźmy do sedna dowodu:






(4.36)

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

(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 hi=(b-a)/2i dla m=2i, 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 , wtedy:





(4.38)

W powyższym dowodzie przyjęto, że pochodna jest równa zero, bo wartości pochodnych f(2m+2)(ξ) w punktach ξi i ξ2 są sobie równe, bo błąd kwadratury (4.35) jest zależny od głównego czynnika , 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[edytuj]

Dotychczas rozpatrywaliśmy wielomiany interpolacyjnego, których kwadratury są określone przez wzór (4.7) przy definicji Ak (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 {Pn(x)}=P0(x), P1(x),...,Pn(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:

(4.39)

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

Twierdzenie
Wielomiany ortogonalne Ps(x) mają jednokrotne pierwiastki rzeczywiste należące do przedziału (a,b)
Dowód
Niech pierwiastkami wielomianu Ps(x) są pierwiastkami nieparzystej niejednokrotności, i oznaczmy te pierwiastki symbolami ξ1, ξ2,..,ξn, i niech m<k, wtedy możemy zbudować wielomian W(x)=(x-ξ1)...(x-ξn), który możemy rozłożyć w bazie funkcji ortogonalnych P0(x),.., Pm(x) wedle sposobu:
(4.40)

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

(4.41)
Z drugiej jednak strony wyrazenie (4.41) jest większe od zera, bo Pk posiada pierwiastki o nieparzystej krotności, także i W(x), zatem wyrażenie W(x)Pk(x) posiada potęgi parzyste przy (x-ξi), wtedy powinno być (W,Pk)>0, zatem dochodzimy do sprzeczności, stąd wynika, że wielomian Pk(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 Ak jest rzędu 2(N+1), gdy xk są pierwiastkami wielomianu PN+1 dla ciągu wielomianów ortogonalnych {Pn(x)}.
Dowód
1)Weźmy sobie wielomian , który jest stopnia 2(N+1), który jest zdefiniowany o wielomian ωN+1=(x-0)..(x-xN), wtedy możemy powiedzieć:
(4.42)
(4.43)

(4.43) jest równe zero dla ωn+1 liczona w punktach xk, 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 {Pn(x)} dla aN+1≠0, wtedy mamy:

(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:

(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:

(4.46)

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


(4.47)

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

(4.48)
Twierdzenie
Współczynniki Ak 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:
(4.49)

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

(4.50)
Ponieważ Ri(xi) jest zawsze nieujemne na mocy (4.49), to Ai wedle (4.50) jest dodatnie.

Weźmy sobie punkty , ,..,,,,...,, takie, że , ale by zachodziło , wtedy błąd interpolacyjny funkcji interpolującej względem funkcji interpolowanej wyrażamy przy pomocy wzoru (1.24), zatem piszemy ją:

(4.51)

Błąd kwadratury (4.7), przy pomocy wagi p(x), piszemy podobnym wzorem do (4.8) dla wielomianu , który jest rzędu 2N+2:

(4.52)

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

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:

(4.53)

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

(4.54)
(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>:

(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 wynikającą z (4.56):





(4.57)

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

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ą:

(4.58)

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

(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):

(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[edytuj]

W powyższych wzorach przyjmijmy α=β=-1/2, wtedy powyżej zdefiniowana norma przechodzi w p(x)=(1-x2)-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 xk i współczynniki Ak określamy przez:

(4.61)
(4.62)

Całkowanie interpolacyjne dla przedziału nieskończonego[edytuj]

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:

(4.63)

wtedy powyższa waga przy zbieżności całki , gdzie W(x) jest wielomianem dowolnego stopnia, wynika że całka określająca współczynniki Ak jest zbieżna, wtedy możemy określić współczynniki Ak i funkcję błędu zależnego od N i funkcji f:

(4.64)
(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:

(4.66)

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

(4.67)
(4.68)

Złożone metody całkowania interpolacyjnego Gaussa[edytuj]

Będziemy tutaj stosowali kwadratury Gaussa-Legendre'a, którego błąd kwadratury dla przedziału , który określamy na podstawie wzoru (4.57), pisaną:


(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:

(4.70)

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


(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):


(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):

(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:

(4.74)