Całkowanie numeryczne

Wstęp

Należy przybliżyć całkę



używając wartości funkcji f w punktach równoległych, przyjmując x₀=a, x_i=x₀+ih, x_n=b; stąd n=(b-a)/h

Wzór prostokątów

Polega na przybliżaniu funkcji f(x) funkcjami, które są przedziałami liniowe

(1)

We wzorze (1) użyto wartości funkcji f w środkach przedziałów (x_i-1,x_i). Jest to istotne dla uzyskania dobrej dokładności, co widać na poniższym rysunku

Wzór trapezów

Podobnie jak wzór prostokątów metoda ta polega na przybliżaniu funkcji f(x) funkcjami, które są przedziałami liniowe

(2)

Wykorzystuje on interpolacje liniową w poszczególnych podziałach (xi-1,xi) (i=1,2,...,n).

Jeśli U jest jednostką na ostatniej zachowanej pozycji dziesiętnej w wartościach funkcji, to dla obu wymienionych wzorów zachodzi nierówność

Błąd obcięcia dla wzorów trapezów

Zarówno dla wzoru prostokątów, jak i dla wzoru trapezów łatwo otrzymać ścisłe oszacowania błędu obcięcia (błędu dyskretyzacji). Załóżmy, że pochodna f"" jest ciągła w [a,b]. Dla jednego kroku wzoru trapezów mamy równość



gdzie


zależy od x.

Ponieważ jednak (x-x_i-1)(x-x_i)<0 dla x჎(x_i-1,x_i), więc z twierdzenia całkowego (w uogólnionej postaci) o wartości średniej wynika, że



, gdzie




Przyjmując, że x=x_i-1+ht, otrzymujemy wzór



Globalny błąd obcięcia jest sumą błędów lokalnych. Dla wzorów trapezów



gdzie


Dla wzorów prostokątów otrzymujemy, korzystając z reszty wzoru Taylora, podobne wyrażenie:



Tak więc wzór prostokątów jest dokładniejszy od wzoru trapezów; ten ostatni jednak jest bardziej ekonomiczny, gdy stosuje się ekstrapolacje iterowaną Richardsona.

Metoda Romberga

Należy ona obecnie do najczęściej stosownych, miedzy innymi dlatego, że daje prostą strategię automatycznego wyboru właściwej długości kroku. Zaczyna się od kroku dość dużego i zmniejsza się go dwukrotnie tyle razy, aby ekstrapolacja dała wartości dostatecznie bliskie. Trzeba przy tym sprawdzić, czy zgodność wartości z tej samej kolumny nie jest przypadkowa, tzn. upewnić się czy krok h jest na tyle mały, że uprawnia do użycia rozwinięcia

2002 Łukasz Przewoźnik  

Schemat Hornera

Częstym zadaniem rachunkowym jest obliczanie w danym punkcie z wartości wielomianu, np.

p(z)=a₀z³+ a₁z²+ a₂z+ a₃

Tę sumę można napisać inaczej, w postaci

p(z)=(( a₀z+ a₁)z+ a₂)z+ a₃

a stąd wynika schemat Hornera, który dla obliczeń ręcznych można sformułować tak:

Przykład

Obliczyć p(8), gdzie p(x)=2x³+x+7



p(8)=1039

Schemat Hornera dla obliczania wartości wielomianu n-tego stopnia,

p(x)=a₀xⁿ+ a₁x^n-1+...+ a_n-1x+a_n

w punkcie z wyraża się wzorem rekurencyjnym

b₀=a₀,         b_i=a_i+zb_i-1        (i=1,2,...,n),     p(z)=b_n

Jeśli pośrednie b_i nie są interesujące, to w większości języków programowania algorytm można opisać bez wskaźnika w b_i, tak jak na poniższym schemacie blokowym

Niekiedy jednak wielkości b_i mają samoistne znaczenie. Wynika to z następującego twierdzenia o dzieleniu wielomianu przez dwumian (tzw. dzieleniu syntetycznym):

Twierdzenie

Zachodzi równość




gdzie b_i są określone przez wzór

b₀=a₀,       b_i=a_i+zb_i-1        (i=1,2,...,n),    p(z)=b_n

Dowód. Oznaczamy prawą stronę tej równości symbolem g(x). Wtedy




stąd




co należało udowodnić.

2002 Łukasz Przewoźnik  

Iteracja

Iteracja jest jednym z najważniejszych narzędzi zaróno praktycznego, jak i teoretycznego badania problemów liniowych i nieliniowych.
Iteracja (z łacińskiego "iteratio", powtarzanie), czyli kolejne przybliżanie oznacza powtarzanie pewnej wzorcowej czynności lub procesu. Iteracja w tym sensie występuje wtedy, gdy wielokrotnie stosuje się jakiś proces numeryczny - być może bardzo złożony i zawierający wiele przykładów użycia iteracji w węższym, niżej opisanym sensie - po to, aby stopniowo ulepszać wcześniejsze wyniki. Aby zilustrować bardziej szczególny użytek z zasady iteracji rozważymy zadanie rozwiązania równania postaci

x=F(x)                                                                                                               (1)

F jest tu funkcją różniczkowalną, której wartości można obliczać dla dowolnej wartości zmiennej rzeczywistej x (z pewnego przedziału). Stosując metodę iteracyjną, zaczynamy od przybliżenia początkowego x₀ i obliczamy ciąg

x₁=F(x₀),      x₂=F(x₁),      x₃=F(x₂),     ...

Każde obliczenie typu

x_n+1=F(x_n)

nazywamy iteracją. Jeżeli ciąg {x_n} jest zbieżny do wartości granicznej , to mamy lim F(x_n)=F(), czyli x równe  spełnia równanie x=F(x). Gdy n rośnie, spodziewamy się, że liczby x_n są coraz lepszymi przybliżeniami szukanego pierwiastka. Iteracja przerywa się po osiągnięciu dostatecznej dokładności.

Interpretacje geometryczną pokazano poniżej

Pierwiastek równania (1) jest odciętą (i rzędną) punktu przecięcia krzywej y=F(x) i prostej y=x . Startując z punktu (x₀,f(x₀)) i używając iteracji, otrzymujemy x₁=F(x₀); punkt x₁ na osi x dostajemy, rysując linie poziomą z punktu (x₀,F(x₀))=(x₀,x₁) do przecięcia z prostą y=x w punkcie (x₁,x₁). Stąd prowadzimy linię pionową do (x₁,F(x₁))=(x₁,x₂) itd. Z powyższego rysunku jest oczywiste, że ciąg {x_n} jest zbieżny monotonicznie do .

Na rysunku poniżej pokazano przypadek, gdy F jest funkcją malejącą

Tym razem zbieżność jest monotoniczna: kolejne przybliżenia x_n leżą na przemian po lewej i po prawej stronie pierwiastka .

Istnieją też jednak pzrypadki rozbieżności, pokazane na przykładach poniżej

Z rysunków tych widać, że wielkością określającą prędkość zbieżności (lub rozbieżności) jest nachylenie krzywej y=F(x) w otoczeniu pierwiastka. Rzeczywiście, z twierdzenia o wartości średniej mamy




gdzie _n leży między x_n-1 i x_n. Wobec tego zbieżność jest tym szybsza, im mniejsze jest |f '(x)| w otoczeniu pierwiastka. Zbieżność jest pewna, jeśli |f '(x)|<1 dla każdego x z otoczenia pierwiastka, które zawiera x₀ i x₁. Jeśli jednak |f '()|>1, to x_n dąży do  tylko w bardzo szczególnych przypadkach, niezależnie od tego, jak blisko  wybierze się do x₀ (różne do ).

Przykład

Szybka metoda obliczania pierwiastków kwadratowych.

Równanie x²=c można napisać w postaci x=F(x), gdzie


             c>0

Pierwiastkiem jest =c^1/2, jest też f '()=0. Przyjmujemy więc




Dla c=2, x₀=1.5 otrzymujemy x₁=1/2(1.5+2/1.5)=1.4167, x₂=1.414216; porównajmy z tym wartość



Dobrą wartość dla x₀ możemy otrzymać za pomocą suwaka logarytmicznego, ale jak możemy sprawdzić - wystarczy znacznie grubsze przybliżenie. Istotnie, dowodzi się, że jeśli x_n ma t cyfr dokładnych, to x_n+1 ma takich cyfr co najmniej 2t-1. Powyższą metodę iteracyjną obliczania pierwiastków kwadratowych stosuje się powszechnie i w kalkulatorach kieszonkowych, i w komputerach.

2002 Łukasz Przewoźnik  

Metoda Newtona

Metoda Newtona jest połączeniem idei iteracji i lokalnej aproksymacji liniowej. W tej metodzie iteracyjnej x_n+1 jest określone jako odcięta punktu przecięcia z osią x stycznej do krzywej y=f(x) w punkcie (x_n,f(x_n)).

Przybliżanie krzywej y=f(x) styczną do niej w punkcie (x₀,f(x₀)) jest równoważne zastąpieniu funkcji dwoma początkowymi składnikami jej szeregu Taylora w punkcie x=x₀. Odpowiednie przybliżanie dla funkcji wielu zmiennych ma również ważne zastosowania.

Do wyznaczenia x_n+1 służy równanie




Metodę Newtona określa następujący wzór iteracyjny:


            gdzie  


Iteracje można przerwać, gdy |h_n| stało się mniejsze od dopuszczalnego błędu pierwiastka (trzeba tu uwzględnić również błędy zaokrągleń i inne błędy popełniane w obliczaniu h_n).

Różnym od kreślenia stycznej sposobem lokalnej aproksymacji krzywej jest wybór dwóch sąsiednich punktów na niej i przybliżanie jej sieczną łączącą te punkty

Przykład

Dana jest funkcja f(x)=sinx-(0.5x)². Wtedy f '(x)=cosx-0.5x. Chcemy obliczyć dodatni pierwiastek z 5 poprawnymi cyframi ułamkowymi.

Przyjmijmy x₀=1.5. Wyniki obliczeń podano w tablicy (jej ostatnią kolumnę można wypełnić dopiero po obliczeniu pierwiastka).




Szukanym pierwiastkiem jest 1.93375. Potrzebne były tylko cztery iteracje, mino że początkowe przybliżenie było dość kiepskie. Widzimy też, że |h_n| maleje coraz szybciej aż do momentu, gdy zaczynają dominować błędy zaokrągleń.

Oczywiście, gdy przybliżenie jest dalekie od pierwiastka, nie trzeba (w obliczeniach ręcznych) używać zbyt wielu cyfr. W powyższym przykładzie trzeba tylko obliczyć f(x_n) z tyloma cyframi ułamkowymi, ile mogłoby być poprawnych w x_n+1, gdyby h_n obliczono dokładnie. Co ważniejsze, ze zbyt dużą dokładnością obliczano pochodną. Wartość f '(x_n) jest tylko środkiem do obliczenia poprawki




Dlatego wystarczy obliczać f '(x_n) z taką dokładnością względną, z jaką oblicza się f(x_n). Ponieważ dokładność względna wartości f(x_n) maleje, gdy x_n zbliża się do pierwiastka, więc w powyższym przykładzie można by użyć wartości f '(x₂) nawet dla n=3 i n=4. Takie uproszczenia mogą być czasem ważne nawet w obliczeniach komputerowych.

Zbieżność metody Newtona

Zakładamy, że funkcje f(x) ma drugą pochodną ciągłą i że szukany pierwiastek  jest pojedynczy. Wobec tego f '()0 i f '(x)0 dla wszystkich x z pewnego otoczenia tego pierwiastka. Niech _n będzie błędem przybliżenia x_n:

_n=x_n-

Przy powyższych założeniach możemy otrzymać zależność między _{n i} _n+1. Rozwijając funkcję f w szereg Taylora w otoczeniu x_n, otrzymujemy

0=f()=f(x_n)+(-x_n)f '(x_n)+0.5(-x_n)²f "() gdzie int(x_n,)

czyli, po podzieleniu przez f '(x_n),




Stąd




a dla x_n

2002 Łukasz Przewoźnik  

Metoda Eulera

Metoda Eulera (metoda siecznych) jest najprostszą metodą przybliżoną dla zagadnienia początkowego


        
                                                                      (1)

z warunkiem początkowym

y(x₀)=y₀

Aby rozwiązać to równanie, musimy znaleźć wartość funkcji y=y(x) w x=x₁=x₀+h.

Załóżmy, że y^*₁=y(x₁) i krzywa (A₀A^*₁) jest rozwiązaniem równania (1). Wykorzystajmy informacje, jakie niesie ze sobą równanie (1). Zauważmy, że znając punkt A₀=(x₀,y₀) (warunek początkowy) możemy policzyć wartość funkcji f w tym punkcie. Jest to współczynnik kierunkowy stycznej do rozwiązania równania (1) w x₀. Możemy więc wyznaczyć tę styczną (prosta (A₀S₀) na rysunku powyżej). Wyznaczmy punkt przecięcia tej stycznej z prostą x=x₁, oznaczmy go przez A₁=(x₁,y₁) i przyjmijmy, że y₁=y(x₁). Przy takim postępowaniu popełnimy błąd w wyznaczaniu wartości y(x₁), na rysunku jest on zaznaczony jako błąd metody. Możemy również przypuszczać, że błąd ten będzie malał wraz z h. Zapiszmy sposób liczenia y₁ w postaci wzoru. Z trójkąta (A₀A₁P) mamy




a stąd


                                                                                         (2)

Chcąc wyznaczyć wartość funkcji y w punkcie x₂=x₁+h=x₀+2h postępujemy analogicznie, tylko zamiast korzystać z warunku początkowego korzystamy z obliczonych uprzednio wartości (x₁,y₁)=(x₁, y(x₁)). Uogólniając to postępowanie wzór powyższy zapiszemy jako

y_i+1=y_i+hf(x_i,y_i), i=0,1,2,...,                                                                                 (3)

gdzie y_i+1 jest wartością funkcji będącej rozwiązaniem równania różniczkowego w x_i+1 czyli y_i+1=y(x_i+1), x_i+1=x₀+(i+1)h dla i=0,2,...

Metode Eulera możemy skonstruować również w inny sposób. A mianowicie zastąpmy funkcję y=y(x) w x=x₀+h jej rozwinięciem w szereg Taylora (zakładamy odpowiednią regularność y).




Jeśli w rozwinięciu Taylora wykorzystamy tylko dwa pierwsze człony, to otrzymujemy




Otrzymaliśmy zatem wzór (2), a powtarzając takie postępowanie dojdziemy do (3).

Błąd metody

Załóżmy, że znamy wartość dokładną rozwiązania w x_i, czyli wartość funkcji y w x_i, i wynosi ona y^*_i=y(x_i). W x_i+1 rozwiązanie równania różniczkowego przyjmuje wartość z (3)


                                                                                         (4)

Z rozwinięcia w szereg Taylora funkcji y w punkcie x_i+1 otrzymujemy wartość dokładną rozwiązania w x_i+1


                                          (5)

Odejmując stronami równania (4) i (5) otrzymujemy błąd lokalny




Wielkość  nazywamy błędem lokalnym, ponieważ założyliśmy, że do obliczenia wartości funkcji y w x_i+1 wykorzystujemy dokładną wartość tej funkcji w x_i. Jeżeli do obliczenia wartości funkcji y w x_i+1 wykorzystamy przybliżoną wartość y_i, to można wykazać, że błąd całkowity metody E=O(h). Przyjmuje się rząd całkowitego błędu jako rząd metody, czyli metoda Eulera jest rzędu pierwszego.

Błąd metody maleje wraz z odpowiednią potęgą h (dla metody Eulera wraz z h). W dotychczasowych rozważaniach na temat błędu nie uwzględniliśmy nieuniknionego błędu zaokrągleń. Okazuje się, że zmniejszenie kroku h powiększa błędy zaokrągleń. Zależność błędu zaokrągleń i błędu metody do kroku h oraz wynikający z nich błąd sumaryczny są przedstawione poniżej

Dla ostatecznie małego h dominującym błędem jest błąd zaokrągleń (nieunikniony). Dalsze zmniejszanie kroku nie będzie poprawiało wyników, może je natomiast pogorszyć.

Przykład

Rozwiązać problem początkowy




metodą Eulera. Rozwiązanie ścisłe ma postać




Weźmy krok h=0.1. Wartość funkcji y będziemy liczyli dla x_i+1=x_i+h, i = 0,1,2,...,9 , i = 0 x₁ = 0.1

y₁ = y₀+hf(x₀,y₀) = y₀+h(-y₀²) = 1+0.1(-1²) = 0.9   x₂=0.2

y₂ = y₁+hf(x₁,y₁) = 0.9+0.1(-0.9²) = 0.819

itd.

Wyniki dalszych obliczeń jak i wartości funkcji będącej ścisłym rozwiązaniem zapisane są w poniższej tabeli

i	x	y	y^*=1/(1+x)
0 1 2 3 4 5 6 7 8	.1 .2 .3 .4 .5 .6 .7 .8 .9 1.0	.9000000 .8190000 .7519239 .6953850 .6470290 .6051643 .5685419 .5362179 .5074649 .4817129	.9090909 .8333333 .7692307 .7142857 .6666667 .6250000 .5882353 .5555555 .5263157 .4999999

2002 Łukasz Przewoźnik  

Metody typu Runge-Kutty

Zajmiemy się równaniem


       
       y(x₀)=y₀                                                     (1)

Podobnie jak przy konstrukcji metody Eulera postawmy sobie następujące zadanie: znaleźć wartość funkcji y=y(x), w x=x₁=x₀+h

Załóżmy, że rozwiązanie zadania (1) jest parabolą i skorzystajmy z następującej własności paraboli:

- współczynnik kierunkowy siecznej (A₀A₁) jest średnią arytmetyczną współczynnika kierunkowego stycznej (A₀S₀) i stycznej (A₁S₁).

Współczynnik kierunkowy stycznej (A₀S₀) umiemy wyznaczyć: jest to f(x₀,y₀). Wyznaczmy przybliżoną wartość funkcji y, w x=x₁ tak, jak to robiliśmy w metodzie Eulera, czyli




Teraz możemy wyznaczyć przybliżoną wartość współczynnika kierunkowego stycznej (A₁S₁), wynosi on




Następnie wyznaczmy współczynnik kierunkowej siecznej (A₀A₁), jest on równy




Jest to przybliżona wartość tego współczynnika ze względu na sposób liczenia y₁^E. Zamiast siecznej (A₀A₁) mamy sieczną (A₀A₁^'). Z trójkąta A₀PA₁^' policzymy






Uporządkujemy dotychczasowe rozwiązania wprowadzając oznaczenia


                                                       (2)

Wzory (2) są wzorami Runge-Kutty II rzędu. Można wykazać, że całkowity błąd metody jest rzędu O(h²).

Spróbujmy dojść do tych samych wzorów (2) inną drogą i wyprowadzić wzory Runge-Kutty wyższych rzędów. Zapiszmy uogólnienie algebraiczne wzorów (2)


                                                                   (3)

Będziemy szukać stałych a,b,c,... , ,,γ,... i w₁,w₂,w₃,... takich, by wartość funkcji y(x) dla x=x₁, czyli y₁ była możliwie bliska wartości dokładnej. Uzyskamy to rozwijając w szereg Taylora funkcję przybliżoną i ścisłą, żądając zgodności (do wyrazu odpowiedniego rzedu) tych rozwinięć. Wykonajmy obliczenie dla metody Runge-Kutty II rzędu, czyli wartość przybliżoną obliczamy następująco


                                                                               (4)

stąd otrzymujemy




Zastąpmy f(x₀+ah,y₀+k₁) jej rozwinięciem w szereg Taylora, czyli


                                                     (5)

gdzie




Z (5) po uporządkowaniu i wstawieniu k₁=hf mamy


                                                 (6)

Załóżmy, że y₀^*=y(x₀) i y₁^*=y(x₁) są dokładnymi wartościami funkcji y i zapiszmy korzystając z rozwinięcia w szereg Taylora


                                             (7)

Żądamy zgodności obu (6, 7) rozwinięć do drugiego wyrazu rozwinięcia w szereg Taylora. Porównajmy więc współczynniki przy h i h² z obu rozwinięć. Otrzymujemy

w₁+w₂=1 ,

w₂a=1/2 ,

w₂=1/2 ,

stąd




Mamy więc całą rodzinę wzorów Runge-Kutty II rzędu zależną od parametru a. I tak dla a=1 mamy =1, w₁=1/2, w₂=1/2, są to wzory (4).

Wartość funkcji y w x=x₂=x₁+h=x₀+2h możemy policzyć stosując (4), o ile (x₁, y(x₁)=y₁) potraktujmy tak, jak warunek początkowy. Czyli ogólnie (4) zapiszemy




Wzory wyższych rzędów można otrzymać analogicznie. W poniższej tabeli są zapisane najczęściej używane współczynniki dla metod typu Runge-Kutty rzędu od I do IV

Rząd metody P	Stałe (wzór 3) w	Wartości współczynników k	Nazwa metody	Błąd lokalny 	Błąd całkowity E
1	w1=1	k1=hf(xi,yi)	Eulera	O(h^2)	O(h)
2	w1=w2=1/2	k1=hf(xi,yi) k2=hf(xi+h,yi+k1)	Heuna	O(h^3)	O(h^2)
3	w1=w3=1/6 w2=2/3	k1=hf(xi,yi) k2=hf(xi+h/2,yi+k1/2) k3=hf(xi+h,yi-k1+2k2)	pokrewna metodzie Simsona	O(h^4)	O(h^3)
4	w1=w4=1/6 w2=w3=1/3	k1=hf(xi,yi) k2=hf(xi+h/2,yi+k1/2) k3=hf(xi+h/2,yi+k2/2) k4=hf(xi+h,yi+k3)	klasyczna metoda Runge-Kutty	O(h^5)	O(h^4)

Metody typu Runge-Kutty są łatwe do zaprogramowani, zmiana kroku całkowania może być dokonywana w dowolnym etapie obliczeń i nie wymaga dużego nakładu pracy, są metodami samostartującymi, tzn. znajomość warunku początkowego wystarcza, by rozpocząć obliczenia. Wymagają one jednak wielokrotnego (dla metody rzędu p p-krotnego) obliczania wartości funkcji f w każdym kroku całkowania, mogą więc być metodami kosztownymi (zwykle najbardziej czasochłonnymi, a więc i kosztownym zadaniem jest obliczanie wartości funkcji f).

Przykład

Znaleźć rozwiązanie zagadnienia początkowego

y'(x)=x+y

y(0)=1

metodą Runge-Kutty IV rzędu. Obliczenia wykonać dla x჎[0,0.2] z krokiem h=0.1 . Porównać otrzymane wyniki z rozwiązaniem ścisłym y(x)=2e^x-x-1.

Wzory Runge-Kutty IV rzędu (patrz (3) i powyższa tabela) mają postać




gdzie




Dla przypomnienia y_i+1=y(x_i+1)=y(x₀+(i+1)h). Dla i=0 liczymy wartości y₁=y(x₁)=y(x₀+h)=y(0.1). Najpierw policzymy




dalej




Analogicznie obliczenia przeprowadzamy dla i=1 (czyli dla x₂=x₀+2h), policzymy y₂=y(x₂). Policzymy wartości rozwiązania ścisłego dla x=x₁ i x=x₂ i wyniki tych obliczeń zestawimy w tabeli

i+1	x	y	k=0.1(x+y)	Y(x)=2e^x-x-1
0	0	1
	0. 0.05 0.05 0.1	1. 1.05 1.055 1.1105	k1=0.1 k2=0.11 k3=0.1105 k4=0.12105
1	0.1	1.11034		1.1103420
	0.1 0.15 0.15 0.2	1.11034 1.170857 1.17638285 1.242978285	k1=0.121034 k2=0.13220857 k3=0.132638285 k4=0.1442978285
2	0.2	1.242803		1.2428055

2002 Łukasz Przewoźnik  

Interpolacja

Niech funkcja f(x) w rozpatrywanym przedziale [a,b] zadana będzie w postaci zbioru dyskretnych wartości

f₀=f(x₀), f₁=f(x₁), f₂=f(x₂), ... , f_n=f(x_n),

wyznaczonych odpowiednio w (n+1) punktach x₀, x₁, x₂, ...,x_n leżących w przedziale [a,b] i zwanych węzłami interpolacji.

Interpolacja funkcji f(x) polega na poszukiwaniu takiej funkcji f^*(x), która w punktach x_i, i=0,1,2,...,n przyjmuje te same wartości co funkcja f(x), tzn.

f^*(x_i)=f(x_i), i=0,1,2,...,n.

Sformułowane tak ogólnie zagadnienie interpolacji może mieć nieskończenie wiele rozwiązań albo też nie mieć ich wcale. Zagadnienie to natomiast staje się jednoznacznie określone, jeśli zamiast dowolnej funkcji f^*(x) szukamy jej w klasie wielomianów określonego stopnia.

Innymi słowy poszukujemy takiego wielomianu P_n(x), zwanego wielomianem interpolacyjnym, stopnia co najwyżej n, który spełnia warunek

P_n(x_i)=f(x_i), i=0,1,2,...,n.

Zauważmy, że na wielomianach interpolacyjnych łatwo możemy wykonywać niektóre działania algebraiczne, takie jak dodawanie, mnożenie, różniczkowanie, całkowanie itp. W kontekście zastosowania interpolacji szczególnie istotne będzie dla nas całkowanie oraz różniczkowanie.

W interpolacji geometrycznej zadanie interpolacyjne sprowadza się do znalezienia krzywej y=P_n(x) przechodzącej przez punkty (x_i,f(x_i)), 1=0,1,2,...,n (co widać na rysunku poniżej)

Postawienie i rozwiązanie zadania interpolacji ma na celu uzyskanie pewnej informacji o funkcji, np. wyznaczenie wartości w punktach różnych od punktów zadanych, Jak się przekonamy, błąd z jakim potrafimy wyznaczyć wartość funkcji w wybranym punkcie x jest mniejszy, gdy leży on wewnątrz rozpatrywanego przedziału [a,b]. W takim przypadku mówimy o interpolacji funkcji. Jeśli natomiast punkt x leży na zewnątrz tego przedziału, błąd taki jest większy. Wówczas mamy do czynienia z ekstrapolacja funkcji.

2002 Łukasz Przewoźnik  

Aproksymacja

Funkcja aproksymowana f(x) określona może być w różny sposób. W naszych rozważaniach będzie zadana w postaci zbioru dyskretnych wartości f₀=f(x₀), f₁=f(x₁), f₂=f(x₂),..., f_m=f(x_m), wyznaczonych odpowiednio w punktach x₀, x₁, x₂,..., x_m, zwanych węzłami aproksymacji. Dla ustalenia uwagi oznaczmy przez f^*(x) funkcję aproksymującą, która przybliża f(x), natomiast przez E(f)=f(x)-f^*(x) błąd takiego przybliżenia. W klasycznym przypadku zatem przez aproksymację rozumieć będziemy poszukiwanie, dla danej funkcji f(x), takiej funkcji f^*(x), aby przyjęta norma błędu ||E(f)|| była najmniejsza.

Do najczęściej stosowanych norm należą

norma średniokwadratowa

oraz norma maksymalna

Jak już wspomnieliśmy funkcję f^*(x) zapisujemy w postaci liniowej kombinacji n funkcji, które będziemy nazywać funkcjami bazowymi.

Zatem:

f^*(x)=c₀₀(x)+c₁₁(x)+ c₂₂(x)+ ... +c_n_n(x),                                                    (1)

gdzie ₀(x),₁(x),₂(x),...,_n(x) są ustalonymi funkcjami bazowymi, c_0,c_1,c_2,...,c_n są poszukiwanymi współczynnikami.

Wielomian (1) nazywać będziemy wielomianem aproksymującym.

Jeżeli np. funkcje bazowe są w postaci jednomianów, tj. _k(x)=x^k dla k=0,1,2,...,n, wtedy

f^*(x)=c₀+c₁x+c₂x²+...+c_nxⁿ.

W takim przypadku mówimy o aproksymacji za pomocą jednomianów. Natomiast jeśli wybrane funkcje bazowe należą do klasy funkcji trygonometrycznych postaci {sinkx, coskx} lub wykładniczych {e^kx}, k=0,1,...,n, to taką aproksymację nazywać będziemy odpowiednio trygonometryczną lub wykładniczą. Zgodnie z powyższym, zagadnienie aproksymacji funkcji f(x) polega na zastąpieniu jej wielomianem aproksymującym f^*(x), rozumiane w sensie przyjętej normy ||E(f)||, było jak najmniejsze. Można to osiągnąć, jak się przekonamy, przez odpowiedni dobór w wyrażeniu (1)współczynników c_0,c_1,c_2,...,c_n..

Załóżmy, że chcemy określić wartości tych współczynników tak, aby spełniony był warunek

f^*(x_k)=f(x_k),                  dla k=0,1,2,...,m.                                                            (2)

Tzn. aby funkcje aproksymujące i aproksymowana miały te same wartości we wszystkich węzłach aproksymacji. Prowadzi to do układu m+1 równań liniowych z n+1 niewiadomymi

                               (3)

Zauważmy, że w tym przypadku nie nakładamy żadnych warunków na zachowanie się funkcji f^*(x) poza punktami x₀, x₁, x₂,..., x_m. W innym podejściu do rozwiązania zadania aproksymacji możemy zrezygnować z warunku (2), a w zamian za to żądać, aby funkcja f^*(x) w całym interesującym nas przedziale przebiegała jak "najbliżej" funkcji f(x). Takie sformułowania zadania przybliżenia funkcji prowadzą do różnych ( w zasadzie trzech) ich typów.

 

1. Przybliżenie interpolacyjne.
   
   Jeżeli n=m, to układ (3) ma dokładnie jedno rozwiązanie, a współczynniki c₀, c₁, c₂,..., c_n mają takie wartości, aby w punktach x₀, x₁, x₂,..., x_m funkcja f(x) była zgodna z funkcją przybliżającą f^*(x), tzn. f(x_k)=f^*(x_k), k-0,1,2,...m. Przybliżenie takie nazywamy interpolacją natomiast punkty x₀, x₁, x₂,..., x_m węzłami interpolacji. Jeśli m > n, to w ogólności nie jest spełniona równość (2), wtedy układ (3) ma więcej równań niż niewiadomych, o takim układzie mówimy, że jest nieokreślony. W tym przypadku trzeba zadowolić się rozwiązaniem przybliżonym.

2. Przybliżenie średniokwadratowe (aproksymacja średniokwadratowa).
   
   Wyrażeniem, którego minimum szukamy jest całka z kwadratu różnicy liczonej między f(x) i jej przybliżeniem f^*(x) w zadanym przedziale [a,b] lub sumą kwadratów takich różnic wyznaczoną w dyskretnych punktach x₀, x₁, x₂,..., x_m.

3. Przybliżenie jednostajne (aproksymacja jednostajna)
   
   Poszukujemy minimum wyrażenia będącego maksymalną różnicą między f(x) a f^*(x) w przedziale [a,b] lub w dyskretnych punktach x₀, x₁, x₂,..., x_m.

Poniższy rysunek jest ilustracją idei zadania aproksymacji.

Zauważmy, że postawienie i rozwiązanie powyższych zagadnień ma na celu uzyskanie pewnej informacji o funkcji, np. wyznaczenie jej wartości w danym punkcie xᄅx_k, k=0,1,2,...m.

Weźmy pod uwagę przybliżenie (1.) . Jeżeli punkt x leży wewnątrz zadanego przedziału [a,b], to mówimy o interpolacji. Jeżeli natomiast punkt leży na zewnątrz tego przedziału, to mamy do czynienia z ekstrapolacją. Zwykle błąd interpolacji jest mniejszy niż błąd ekstrapolacji.

Poniższy rysunek przedstawia zakres formułowania tych zadań.

2002 Łukasz Przewoźnik  

---