Metody Numeryczne w Fizyce 1 Aproksymacja i Interpolacja 1
Aproksymacja / Interpolacja
Definicje
Należy dokonać rozróżnienia pomiędzy zagadnieniami interpolacji, aproksymacji i ekstrapolacji. W zależności
od podręcznika można spotkać dość różne określenia. Najwygodniejsze są chyba takie:
1. Interpolacja jest zagadnieniem polegającym na znalezieniu funkcji (łatwej do póz
niejszego manipulowania),
która na tyle dobrze przybliżałaby inną, nieznaną, funkcję, której wartości w pewnych punktach xi znamy, aby
umożliwić wykorzystanie jej do obliczenia wartości f(x) dla dowolnego x nie będącego xi. Najczęściej punkty xi
są nam dane i nie mamy wpływu na ich rozmieszczenie (bo np. są wynikiem pomiarów). Jeśli interesujący nas
punkt x leży gdzieś pomiędzy punktami xi to zagadnienie takie nazywamy zagadnieniem interpolacyjnym, jeśli
zaś leży na zewnątrz ich zbioru to mamy do czynienia z ekstrapolacją.
2. Aproksymacja jest zagadnieniem analogicznym do poprzedniego, z tą jednak różnicą, że najczęściej znamy
funkcję wyjściową. Zaś funkcja aproksymująca jest nam potrzebna np. dlatego, że obliczanie wartości funkcji
wyjściowej jest zbyt czasochłonne. W tej sytuacji możemy sami zadecydować o rozmieszczeniu punktów, w
których wartości funkcji wykorzystamy do utworzenia funkcji aproksymującej.
Istnieją też inne określenia. Np.: zagadnienie interpolacji polega na znalezieniu funkcji (najczęściej wielomianu
odpowiednio wysokiego stopnia), która będzie przebiegać przez wszystkie znane punkty. Ma to pewną
niedogodność, ponieważ w przypadku wielu punktów (np. wyników pomiarów) wielomian musiałby być bardzo
wysokiego stopnia. Poza tym wartości wielomianu w punktach znajdujących się pomiędzy punktami znanymi
mogłyby się znacznie różnić od wartości nieznanej funkcji f (np. funkcja liniowa, dla której zmierzono 10 wartości
(obarczonych błędem pomiarowym) byłaby opisana wielomianem stopnia 9-go). Z kolei zagadnienie aproksymacji
polega na znalezieniu takiej funkcji, której wartości będą z zadaną dokładnością zgodne z pewną liczbą znanych
wartości innej funkcji (najczęściej nieznanej)  nie muszą się więc dokładnie z nimi pokrywać. Funkcją
aproksymującą jest najczęściej wielomian, ale mogą to też być funkcje trygonometryczne (najczęściej sin i cos) albo
funkcja wykładnicza. Może nią też być iloraz dwóch wielomianów. Najczęściej możliwe jest wybranie stopnia
wielomianu aproksymującego (i jest to stopień znacznie niższy od liczby danych punktów). Warto jednak pamiętać,
że w rzeczywistości wartości funkcji trygonometrycznych czy wykładniczych (które miałyby być użyte do
aproksymacji) będą obliczane za pomocą sum szeregów będących w istocie wielomianami aproksymującymi te
funkcje. Aproksymacja wielomianowa ma dodatkowo jedną ważną zaletę:
Udowodniono, że dla każdej funkcji ciągłej w zadanym przedziale i dla dowolnej dokładności (za-
danej przez dopuszczalne odchylenie µ) istnieje wielomian stopnia n (zależnego od µ), którego
wartoÅ›ci mogÄ… różnić siÄ™ od wartoÅ›ci funkcji o mniej niż µ.
W niektórych podręcznikach nie wprowadza się rozróżnienia pomiędzy interpolacją a aproksymacją określając
oba te zagadnienia jednakowo.
W przypadku interpolacji/aproksymacji zawsze mamy do czynienia z sytuacją, w której nie jest możliwe
naprawdę dokładne przewidzenie popełnianego błędu  zawsze możemy wyobrazić sobie funkcję, której
nieciągłość jest widoczna w przedziale mniejszym niż odległość punktów danych  wówczas funkcja taka będzie
wyglądać na ciągłą i tak zostanie przybliżona, dając w pewnych punktach nieskończenie wielkie błędy.
Wielomiany interpolacyjne
Jeśli danych jest n+1 punktów (x0...xn) to mówimy, że mamy do czynienia z interpolacją rzędu n. Do interpolacji
najczęściej wykorzystywane są wielomiany  będzie to więc wielomian rzędu n.
Wielomian interpolacyjny Lagrange a
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
2 Metody Numeryczne w Fizyce Andrzej Brozi
Załóżmy, że znamy n wartości nieznanej funkcji fi w n punktach xi, a chcielibyśmy poznać jej wartość w
pewnym punkcie x. Istnieje wielomian stopnia n, który będzie mieć we wszystkich punktach xi wartości fi:
n
(x x0) (x xi 1)(x xi 1) (x x )
n
L (x ) fi
n
(xi x0) (xi xi 1)(xi xi 1) (xi x )
i 0
n
Poniżej mamy dwa przykłady bezpośredniego zastosowania wzoru interpolacyjnego Lagrange'a do funkcji
2
f x e x xoraz f x ex :
Na powyższych rysunkach przedstawiono błędy interpolacji Lagrange a dwóch funkcji: po lewej stronie
2
f x e x x, a po prawej f x ex . W obu przypadkach interpolowano wartości funkcji dla 37 punktów
równomiernie rozłożonych na przedziale <0,1.6>. Górne rysunki odpowiadają interpolacji opartej na 5
równomiernie rozłożonych punktach, a dolne  na 10 punktach.
Można udowodnić, że rozwiązanie zadania interpolacyjnego jest tylko jedno, i że da się je przedstawić w
postaci:
p(x ) c0 c1(x x0) c2(x x0)(x x1) c (x x0)(x x1) (x x )
n n 1
co można zapisać również tak:
n i 1
p x c0 ci x xj
i 1 j 0
zaś występujące w nim współczynniki ci można obliczyć rekurencyjnie z następującego układu równań:
Aproksymacja / Interpolacja
Metody Numeryczne w Fizyce 1 Aproksymacja i Interpolacja 3
p(x0) c0
p(x1) c0 c1(x1 x0)
p(x2) c0 c1(x2 x0) c2(x2 x0)(x2 x1)
p(x3) c0 c1(x3 x0) c2(x3 x0)(x3 x1) c3(x3 x0)(x3 x1)(x3 x2)
....................
Wzór ogólny:
n i 1
p(x ) c ci x xj
n o n
i 1 j 0
zaś wzór ogólny na współczynnik cn:
n 1 i 1
p(x ) c ci x xj
n o n
i 1 j 0
c
n
n 1
x xj
n
j 0
Jeśli nie chodzi nam o zapisanie postaci wielomianu, a jedynie o poznanie jego wartości w konkretnym punkcie
x możemy skorzystać z algorytmu Neville a:
Załóżmy, że L0 jest wartością wielomianu stopnia zero (wielkości stałej) przechodzącego przez punkt (x1,f1)  a
więc L1 = f1. Analogicznie zdefiniujemy wielkości L2...Ln. Teraz zdefiniujmy L12 jako wielomian stopnia pierwszego
przechodzÄ…cy przez punkty (x1,f1) i (x2,f2) i analogicznie L23, L34, ... L(n 1)n. Analogiczne definicje wprowadzimy dla
wielomianów wyższych stopni aż do L123...n będącego wartością wielomianu interpolującego przechodzącego przez
wszystkie n punktów, czyli szukaną wielkością. Kolejne L tworzą tablicę, dla n = 4 będzie to:
x1: f1 = L1
L12
x2: f2 = L2 L123
L23 L1234
x3: f3 = L3 L234
L34
x4: f4 = L4
Algorytm Neville a jest rekurencyjnym sposobem znajdowania wartości z powyższej tablicy, kolumnami od
lewej do prawej:
(x xi m)Li (i 1) (i m 1) (xi x )L(i 1)(i 2) (i m)
Li(i 1) (i m)
xi xi m
W przypadku ostatniego elementu powyższej przykładowej tablicy mamy i = 1, m = 3 :
(x x4)L123 (x1 x )L234
L1234
x1 x4
Przykład: Chcemy znalez
ć wartość pewnej funkcji w punkcie x = 0,61. O funkcji tej wiemy, że w punktach x = 0,5;
0,55; 0,6; 0,65; 0,7; 0,75 przyjmuje znane nam wartości (w tabeli poniżej znajdują się wartości funkcji
f(x) = exp( x)  x  znając funkcję będziemy mogli obliczyć dokładnie popełniany błąd). Posłużymy się
schematem Neville a. Zapiszemy następującą tablicę:
Interpolacja funkcji: exp(-x) x, w punkcie x0=0.61 wartość dokł. =  0.06664913
xi fi (=Li) błędy względne
0.5 0.10653066 2.84787e-02 -9.38222e-05 -9.22437e-07 -1.62817e-08 -3.71975e-10
-0.06854721
0.55 0.02694981 -0.06664288
-0.06681600 -0.06664907
0.6 -0.05118836 -0.06665132 -0.06664913
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
4 Metody Numeryczne w Fizyce Andrzej Brozi
-0.06654154 -0.06664918 -0.06664913
0.65 -0.12795422 -0.06664597 -0.06664913
-0.06758584 -0.06664902
0.7 -0.20341470 -0.06669181
-0.06982094
0.75 -0.27763345 4.75897e-02 6.40281e-04 -1.62569e-06 1.98768e-08
błędy względne
Jak widać wartość błędu względnego tak obliczonej wartości w stosunku do wartości dokładnej zależy od ilości
wykorzystanych do obliczeń punktów, a nawet od tego gdzie, w stosunku do wartości znanych, znajduje się wartość
szukana.
W metodzie tej odległości pomiędzy punktami, w których wartości funkcji są znane nie muszą koniecznie być
równe  poniżej znajduje się przykład dla punktów rozmieszczonych nierównomiernie:
Interpolacja funkcji: exp(-x)-x dla punktów rozmieszczonych nierównomiernie
Pozostałe parametry jak w poprzednim przykładzie
xi fi (=Li) błędy względne
0.49 0.1226263942 0.0472303195 -0.000783151 -0.000011459 -0.00000018 -4.9566e-09
-0.069796991
0.52 0.074520548 -0.066596935
-0.067930291 -0.066648367
0.56 0.0112090638 -0.066674084 -0.066649119
-0.065836612 -0.06664962 -0.06664913
0.67 -0.158291422 -0.066627875 -0.066649147
-0.067893896 -0.066648254
0.69 -0.188423931 -0.066717541
-0.070050546
0.78 -0.321593989 0.0510346539 0.001026427 -0.000013163 0.0000002426
błędy względne
Jak widać błędy poszczególnych wartości rozłożone są nieco inaczej niż poprzednio.
Pewną modyfikacją tej metody jest obliczanie niewielkich różnic pomiędzy wartościami L z kolejnych kolumn.
Zdefiniujmy:
C Li (i m) Li (i m 1)
m,i
D Li (i m) L(i 1) (i m)
m,i
Wykorzystując podany wyżej wzór rekurencyjny łatwo można je przekształcić do postaci:
(xi m 1 x )(C D ) (xi x )(C D )
m,i 1 m,i m,i 1 m,i
D ; C
m 1,i m 1,i
xi xi m 1 xi xi m 1
Aproksymacja / Interpolacja
Metody Numeryczne w Fizyce 1 Aproksymacja i Interpolacja 5
Ostateczne rozwiązanie L1...n jest sumą któregokolwiek fi oraz tych C i/lub D, które leżą na drodze od wybranego
fi do ostatecznego L.
Jednak ta metoda interpolacji nie jest zbyt wygodna w praktyce  nie znaleziono dla niej skutecznej metody
szacowania błędów.
Interpolacja funkcjami wymiernymi
Niektóre funkcje, niezbyt dobrze aproksymowane wielomianami, dają się dobrze aproksymować funkcjami
wymiernymi (za pomocą wielomianów nie da się aproksymować funkcji mających w pewnych miejscach
mianownik równy zeru). Istnieje algorytm analogiczny do algorytmu Neville a, autorstwa Bulirscha i Stoera,
umożliwiający znalezienie wartości aproksymowanej funkcją wymierną. Wzór rekurencyjny (analogiczny jak
poprzednio) ma postać:
R(i 1) (i m) Ri (i m 1)
Ri (i 1) (i m) R(i 1) (i m)
x xi R(i 1) (i m) Ri (i m 1)
1 1
x xi m R(i 1) (i m) R(i 1) (i m 1)
Wartościami początkowymi są Ri = fi. Ponieważ jednak wzory te wykorzystują dodatkowo jeszcze jedną
kolumnę (m 1) zatem dla zainicjowania obliczeń trzeba przyjąć R = 0 wszędzie tam, gdzie zachodziłaby
konieczność  sięgnięcia do kolumny m =  1 (która przecież nie istnieje).
Interpolacja funkcji: exp(-x)-x w punkcie x0 = 0,61 - punkty rozmieszczone równomiernie
xi fi (=Li) błąd względny
0.50 0.10653066 -1.21322e+00 -4.21185e-05 4.49466e-08 6.61065e-09 -2.07659e-09
0.01421083
0.55 0.02694981 -0.06664632
-0.03240018 -0.06664913
0.60 -0.05118836 -0.06665016 -0.06664913
-0.05816791 -0.06664913 -0.06664913
0.65 -0.12795422 -0.06664759 -0.06664913
-0.09867111 -0.06664913
0.70 -0.20341470 -0.06667072
-0.13733217
0.75 -0.27763345 1.06052e+00 3.23968e-04 3.15098e-08 -8.00456e-09
błąd względny
Jak widać przybliżenia występujące w pierwszej kolumnie (tzn. oparte tylko na dwóch punktach) obarczone są
bardzo dużymi błędami, wynika to z przyjęcia w ich przypadku wartości z kolumny m 1 wynoszących 0.
Dla punktów rozmieszczonych nierównomiernie mamy:
Interpolacja funkcji: exp(-x)-x - punkty rozmieszczone nierównomiernie
xi fi (=Li) błąd względny
0.49 0.12262639 -1.51362e+00 -3.43613e-04 6.06882e-07 7.31603e-08 -2.50924e-08
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
6 Metody Numeryczne w Fizyce Andrzej Brozi
0.03423260
0.52 0.07452055 -0.06662623
0.00543607 -0.06664917
0.56 0.01120906 -0.06666065 -0.06664914
0.02183867 -0.06664911 -0.06664913
0.67 -0.15829142 -0.06663886 -0.06664912
-0.10697129 -0.06664915
0.69 -0.18842393 -0.06668409
-0.13772840
0.78 -0.32159399 1.06647e+00 5.24474e-04 2.56192e-07 -9.77641e-08
błąd względny
Analogicznie jak poprzednio możemy zdefiniować różnice:
C Ri (i m) Ri (i m 1)
m,i
D Ri (i m) R(i 1) (i m)
m,i
Różnice te spełniają zależność (którą można wykorzystać do udowodnienia wzorów rekurencyjnych poniżej):
C D C D
m 1,i m 1,i m,i 1 m,i
A wzory rekurencyjne na różnice:
x xi
D C Dm,i
C C Dm,i x xi m 1 m,i m,i 1
m,i 1 m,i 1
D ; C
m 1,i m 1,i
x xi x xi
D C D C
x xi m 1 m,i m,i 1 x xi m 1 m,i m,i 1
Metoda obliczeń jest analogiczna jak poprzednio.
Zjawisko Rungego
Jest to zjawisko występujące w przypadku wielomianowej interpolacji funkcji przy dużej ilości równoodległych
węzłów.
Aproksymacja / Interpolacja
Metody Numeryczne w Fizyce 1 Aproksymacja i Interpolacja 7
W powyższym przykładzie wykorzystano funkcję postaci:
2
x
1
2Ã2
f (x ) e
2Ä„Ã
(czyli rozkład Gaussa) dla à = 0,05. Wszystkie wykresy (tzn. wartości dokładnych i interpolowanych) oparte są na
21 punktach  odległych od siebie o 0,1. Punktów, na podstawie których przeprowadzono interpolację było 11 
były odległe od siebie o 0,2). Na wykresie przedstawiony jest przebieg badanej funkcji, jej aproksymacji
wielomianem Lagrange a oraz wynik analogicznej interpolacji przeprowadzonej dla punktów nierównoodległych
(punktów Czebyszewa).
Jak widać przebieg przybliżenia wykorzystującego większą liczbę punktów (a zatem będącego wielomianem
wyższego rzędu) może bardzo odbiegać od przebiegu funkcji przybliżanej. Zjawisko to występuje szczególnie silnie
w przypadku węzłów rozłożonych równomiernie. Dużo lepszą zgodność (choć i tak daleką od ideału) można
uzyskać rozmieszczając węzły nierównomiernie. Na wykresie pokazano wynik zastosowania tzw. węzłów Czeby-
szewa, których położenia opisane są wzorem:
2i 1 Ä„
xi cos i 0,1, ,m
m 1 2
(w naszym przykładzie m = 10  liczba przedziałów, na których oparto wielomian interpolacyjny).
Interpolacja odwrotna
Oba powyżej opisane schematy interpolacyjne można wykorzystać w  przeciwnym kierunku . Załóżmy, że
mamy daną pewną funkcję f(x) i chcemy dowiedzieć się dla jakiej wartości x przyjmuje ona konkretną wartość f (w
szczególności wartością tą może być zero  będziemy wówczas szukać miejsca zerowego tej funkcji). Musimy
założyć oczywiście, że w ogóle istnieje taka wartość x, dla której funkcja przyjmuje zadaną wartość f. W takim
przypadku możemy obliczyć wartości funkcji dla jakiegoś zbioru wartości x, takiego, aby pożądana wartość f
znalazła się między nimi. Z formalnego punktu widzenia nie ma różnicy między zmienną niezależną (x) a zależną
(f(x))  można je więc zamienić miejscami i szukać wartości x odpowiadającej  argumentowi f.
Przykład: Załóżmy, że chcemy rozwiązać równanie e x  x = 0 (jest to ta sama funkcja, którą wykorzystaliśmy w
przykładach interpolacji schematami Neville'a i Bulirscha-Stoera). Tablica wartości funkcji jest następująca:
x1 = 0,50 f1 = 0,10653
x2 = 0,55 f2 = 0,02695
x3 = 0,60 f3 =  0,05119
x4 = 0,65 f4 =  0,12795
x5 = 0,70 f5 =  0,20341
x6 = 0,75 f6 =  0,27763
Szukając miejsca zerowego tej funkcji zamienimy miejscami wartości funkcji i jej argumentu (tym samym
problem poszukiwania miejsca zerowego zmienimy w problem znalezienia za pomocą interpolacji wartości funkcji
do niej odwrotnej w zerze).
Oczywiście także i w tym przypadku można zwiększyć liczbę znanych punktów (podwyższając tym samym
stopień wielomianu interpolacyjnego) i w ten sposób zwiększyć dokładność.
wartość dokł. = 0.5671432904
Szukanie miejsca zerowego funkcji: exp(-x)-x
(interpolacja odwrotna)
f(wart.dokł.)=5.15538e-17
fi xi błędy względne
0.10653066 0.50 -3.71941e-04 3.82197e-07 5.45703e-08 -1.47749e-09 -1.07303e-10
0.56693235
0.02694981 0.55 0.56714351
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
8 Metody Numeryczne w Fizyce Andrzej Brozi
0.56724497 0.56714332
-0.05118836 0.60 0.56714310 0.56714329
0.56665942 0.56714323 0.56714329
-0.12795422 0.65 0.56714422 0.56714329
0.56521772 0.56714395
-0.2034147 0.70 0.56714540
0.56296273
-0.27763345 0.75 -7.37125e-03 3.72778e-06 1.15687e-06 3.46360e-09
błędy względne
Interpolacja za pomocÄ… funkcji sklejanych trzeciego stopnia
Metoda ta jest znana jako splines. Skupmy uwagę na pojedynczym przedziale pomiędzy xi i xi+1. Interpolacja
liniowa w tym przedziale da nam wzór:
f (x) Afi Bfi 1
gdzie:
xi 1 x x xi
A ; B 1 A
xi 1 xi xi 1 xi
Taka funkcja ma na krańcach przedziału nieokreśloną drugą pochodną. Celem interpolacji metodą funkcji
sklejanych jest znalezienie takiej postaci wzoru interpolacyjnego, która miałaby ciągłą drugą pochodną i gładką
pierwszą pochodną, zarówno wewnątrz jak i na krańcach przedziałów.
Załóżmy (co nie jest prawdą), że oprócz wartości funkcji fi mamy dane również wartości f . Moglibyśmy
i
wówczas do prawej strony dodać wielomian trzeciego stopnia taki by jego druga pochodna zmieniała się w sposób
ciągły od wartości f do wartości f . Jeśli na dodatek będzie on tak skonstruowany by w punktach xi i xi+1 miał
i i+1
wartość zero to nie popsuje on zgodności ze znanymi wartościami funkcji fi i fi+1. Otrzymamy wówczas:
f (x) Afi Bfi 1 Cfi Dfi 1
gdzie A i B są już zdefiniowane, zaś C i D:
2
1 3
C A A xi 1 xi
6
2
1 3
D B B xi 1 xi
6
Wykorzystując definicje A, B, C i D możemy obliczyć pierwszą i drugą pochodną f(x):
2
fi 1 fi 3A 2 1
d f 3B 1
xi 1 xi fi xi 1 xi fi 1
dx xi 1 xi 6 6
2
d f
Afi Bfi 1
2
dx
Problemem jest oczywiście fakt, że wartości f nie są znane. Nie wykorzystaliśmy jednak jeszcze postulatu
i
ciągłości pierwszej pochodnej  zażądamy, aby pierwsza pochodna w punkcie xi, obliczona na przedziale (xi 1,xi),
była równa pierwszej pochodnej w punkcie xi obliczonej na przedziale (xi,xi+1). Powtarzając tę operację dla
wszystkich przedziałów otrzymamy układ równań:
xi xi 1 xi 1 xi 1 xi 1 xi fi 1 fi fi fi 1
fi 1 fi fi 1
6 3 6 xi 1 xi xi xi 1
czyli układ n 2 równań z n niewiadomymi (którymi są nieznane wartości drugich pochodnych). Aby uzyskać
jednoznaczne rozwiązanie musimy dodać jeszcze dwa warunki (zazwyczaj podaje się brzegowe wartości f i f ).
1 n
Można przyjąć obie te wartości za zero albo narzucić wartości pierwszych pochodnych na brzegach i z wzorów na
Aproksymacja / Interpolacja
Metody Numeryczne w Fizyce 1 Aproksymacja i Interpolacja 9
pierwszą pochodną obliczyć odpowiednie wartości drugiej pochodnej (można to zrobić na obu albo tylko na jednym
brzegu). Otrzymany w ten sposób układ równań jest opisany macierzą trójdiagonalną (jak widać każda z szukanych
wartości i związana jest tylko ze swymi najbliższymi sąsiadami: i 1 i i+1) co umożliwia szybkie rozwiązanie
(zagadnieniem rozwiązywania układów równań zajmiemy się w dalszej części wykładu).
Zagadnienie to można potraktować też nieco inaczej. Jeśli dane składają się z n punktów to zechcemy je
przedstawić jako połączone odcinki parabol (w liczbie n 1):
3 2
fi x ai x xi bi x xi ci x xi di
Dla każdego z tych segmentów musimy określić cztery parametry  ai, bi, ci i di. Z faktu, że uzyskana przez nas
linia musi przechodzić przez zadane punkty xi wynika, że:
yi fi xi di
yi 1 fi xi 1 ai"xi3 bi"xi2 ci"xi di
gdzie "xi jest szerokością i-tego przedziału ("xi = (xi+1 xi).
Potrzebne będą jeszcze dwa dodatkowe warunki  będą one dotyczyć gładkości otrzymanej krzywej oraz
połączenia (sklejenia) sąsiednich odcinków:
fi 1 xi 1 fi xi 1
fi 1 xi 1 fi xi 1
Ponieważ poszczególne segmenty są parabolami, ich drugie pochodne, które muszą być sobie równe w węzłach,
można wyrazić przez funkcję liniową:
x xi
fi x fi xi fi xi fi xi
"xi 1 1
korzystając z wcześniej zapisanej postaci funkcji fi(x) można zapisać jako:
x xi
fi x 2bi 2 bi 1 bi
"xi
Całkując f (x) (i wykorzystując fakt, że fi (xi) ci) otrzymamy fi (x):
x
x xi 2
fi x fi xi fi x dx ci 2bi x xi bi bi
"xi 1
xi
Wykorzystując zgodność nachyleń (pierwszych pochodnych) w węzłach otrzymamy:
fi 1 xi 1 fi xx 1 ci 1 ci 2bi"xi "xi bi 1 bi
Całkując ponownie otrzymamy funkcje fi(x):
3
x
x xi
2
fi x fi xi fi x dx yi ci x xi bi x xi bi bi
3"xi 1
xi
WykorzystujÄ…c warunek fi(xi+1) = yi+1 otrzymamy:
"yi 1
ci "xi bi 1 2bi
"xi 3
Podstawiając to wyrażenie do wzoru na pochodne otrzymujemy równanie (obowiązujące dla i = 1,...,n 2):
"yi 1 "yi
1
bi "xi 2bi 1 "xi "xi 1 bi 2"xi 1
3 "xi 1 "xi
Z równania tego można obliczyć wszystkie wartości bi (za wyjątkiem b0 i bn, w przypadku których założymy, że
równe są zero). Mając bi możemy z jednego z powyższych równań obliczyć ci, natomiast wyrażenie na ai będziemy
mogli zapisać jako:
bi 1 bi
ai
3"xi
Znając wszystkie współczynniki będziemy mogli wykreślić całą linię.
Aproksymacja
Żądanie aby funkcja interpolacyjna przebiegała przez wszystkie znane punkty może czasami być nierozsądne.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
10 Metody Numeryczne w Fizyce Andrzej Brozi
Jeśli np. wartości punktów są wynikami pomiarów to z pewnością są obarczone pewnymi błędami. W takiej sytuacji
interpolacja (wymagająca dokładnej zgodności wszystkich znanych punktów) mijałaby się z celem. W przypadku
aproksymacji akceptujemy pewne niezgodności funkcji aproksymującej także w punktach o zadanych wartościach i
zadowalamy się ograniczeniem błędów w tych punktach.
Problemem pozostaje jednak sposób zadecydowania, która z możliwych funkcji aproksymujących jest najlepsza.
W tym celu obliczamy odchylenia dla wszystkich znanych punktów a następnie wybieramy tę funkcję interpolującą,
dla której suma kwadratów odchyleń jest najmniejsza (metoda najmniejszych kwadratów) albo tę, dla której
odchylenie maksymalne jest najmniejsze (aproksymacja jednostajna) (ang. minimax).
Aproksymacja metodą najmniejszych kwadratów
Najczęściej chyba stosowaną metodą aproksymacji jest metoda najmniejszych kwadratów (znana chyba
wszystkim z pierwszej pracowni fizycznej. Warto najpierw wyjaśnić dlaczego wybieramy właśnie sumę kwadratów
odchyleń jako  wskaz
nik jakości funkcji aproksymującej.
Założymy, że każda spośród zmierzonych wartości jest obarczona błędem pomiarowym i wykazuje losowy
rozkład Gaussa wokół wartości dokładnej oraz że odchylenia standardowe są we wszystkich punktach jednakowe i
wynoszÄ… Ã. Wówczas prawdopodobieÅ„stwo tego, że wszystkie punkty danych sÄ… odlegÅ‚e od wartoÅ›ci dokÅ‚adnych o
mniej niż "y jest proporcjonalne do następującego iloczynu:
N
yi y xi 2
1
P exp "y
2 Ã
i 1
Maksymalizowanie wartości takiego wyrażenia sprowadza się do zminimalizowania jego logarytmu ze znakiem
minus:
N
yi y xi 2
N log "y
i 1 2Ã2
Natomiast ponieważ N, à i "y sÄ… wielkoÅ›ciami staÅ‚ymi sprowadza siÄ™ to do zminimalizowania:
N
yi y xi;a1...aM 2
i 1
(ai są parametrami wybranej przez nas funkcji przybliżającej dane eksperymentalne, mogą więc być to np.
współczynniki a i b funkcji liniowej). Musimy więc rozwiązać następujący układ równań:
f
0
a1
f
0
a2
f
0
aM
Jak wszyscy powinni pamiętać z zajęć w pracowni fizycznej rozwiązanie takiego układu nie stanowi problemu
gdy poszukiwaną funkcją jest funkcja liniowa czy kwadratowa, jednak gdybyśmy wybrali funkcję o
skomplikowanej postaci, mogłoby się okazać, że rozwiązania symbolicznego znalez
ć nie możemy. Konieczne
byłoby zastosowanie numerycznych metod rozwiązywania układów równań  będziemy o nich mówić póz
niej.
Aproksymacja jednostajna
Celem aproksymacji jednostajnej jest podanie wyrażenia minimalizującego maksymalną wartość błędu, nie
będziemy zajmować się szczególnymi postaciami tego typu aproksymacji, jako że aproksymacja wielomianami Cze-
byszewa daje faktycznie najmniejszą możliwą wartość maksymalnego błędu.
Wielomiany Czebyszewa
Aproksymacja / Interpolacja
Metody Numeryczne w Fizyce 1 Aproksymacja i Interpolacja 11
Aproksymacja za pomocą wielomianów ma pewną niedogodność, mianowicie jej błędy są bardzo małe w środku
przedziału (gdzie mamy do czynienia z potęgowaniem liczb zbliżonych do zera) i znacznie wzrastają w miarę
zbliżania się do krańców przedziału. Zawsze możemy przekształcić każdy przedział badany do zakresu ( 1,1) 
niech a < z < b, wówczas jeśli:
2z b a
x
b a
to  1 < x < 1. Możemy więc, bez utraty ogólności rozważań, ograniczyć się do przedziału ( 1,1).
Wielomiany Czebyszewa zdefiniowane sÄ… jako:
T x cos n arccosx
n
Od razu widać, że T0(x) = 1 i T1(x) = x. WykorzystujÄ…c wzór: cos2¸ = 2cos2¸  1 otrzymamy T2(x) = 2x2  1.
Dalsze wielomiany można obliczać stosujÄ…c wzór: cosn¸ = 2cos¸ cos(n 1)¸  cos(n 2)¸, co prowadzi do wzoru
rekurencyjnego:
T x 2xT x T x
n n 1 n 2
Kolejne wielomiany Czebyszewa mają postać:
T0(x) = 1
T1(x) = x
T2(x) = 2x2  1
T3(x) = 4x3  3x
T4(x) = 8x4  8x2 + 1
T5(x) = 16x5  20x3 + 5x
T6(x) = 32x6  48x4 + 18x2  1
T7(x) = 64x7  112x5 + 56x3  7x
T8(x) = 128x8  256x6 + 160x4  32x2 + 1
T9(x) = 256x9  576x7 + 432x5  120x3 + 9x
Z samej definicji wielomianów Czebyszewa wynika, że wszystkie one mają wartości ekstremalne wynoszące ą1
rozmieszczone w różnych punktach przedziału ( 1,1), co zapewnia bardziej równomierny rozkład błędów
aproksymacji.
Wielomiany Czebyszewa wykorzystuje się do  skracania rozwinięć w szeregi potęgowe. Rozpatrzmy to na
następującym przykładzie:
Chcemy znalez
ć aproksymację funkcji cosinus w przedziale ( 1,1) z dokładnością nie gorszą niż 5 10 5.
Wyjdziemy od zwykłego rozwinięcia w szereg potęgowy:
2 4 6 8 10
x x x x x
cosx 1
2! 4! 6! 8! 10!
W przypadku tego rodzaju szeregów przemiennych błąd spowodowany obcięciem szeregu po którymś składniku
ma wartość bezwzględną mniejszą od wartości bezwzględnej składnika obciętego. Jeśli więc odrzucimy składnik
x8/8! 2.48 10 5 to otrzymamy zadowalającą nas dokładność. Zakładając, że wartości odwrotności silni mamy
stablicowane oraz, że mamy dość miejsca w pamięci na przechowanie wyników pośrednich, musimy wykonać 8
mnożeń i 4 dodawania. Spróbujemy teraz zapisać otrzymane wyrażenie wykorzystując wielomiany Czebyszewa.
Składając odpowiednio wielomiany Czebyszewa możemy za ich pomocą zapisać kolejne potęgi x:
1 = T0
x = T1
x2 = (T0 + T2) / 2
x3 = (3T1 + T3) / 4
x4 = (3T0 + 4T2 + T4) / 8
x5 = (10T1 + 5T3 + T5) / 16
x6 = (10T0 + 15T2 + 6T4 + T6) / 32
x7 = (35T1 + 21T3 + 7T5 + T7) / 64
x8 = (35T0 + 56T2 + 28T4 + 8T6 + T8) / 128
x9 = (126T1 + 84T3 + 36T5 + 9T7 + T9) / 256
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
12 Metody Numeryczne w Fizyce Andrzej Brozi
Wykorzystując powyższe wyrażenia otrzymamy:
1 1 1 1
x T0 T0 T2 3T0 4T2 T4
2! 2 4! 8
1 1
10T0 15T2 6T4 T6
6! 32
1 1
35T0 56T2 28T4 8T6 T8
8! 128
Grupując współczynniki przy Ti otrzymamy:
c o s x 0,76519775T0 0,22980686T2 0,0049533419T4
4,185265 10 5T6 1,937624 10 7T8
Ponieważ maksymalna wartość wielomianów Czebyszewa w przedziale ( 1,1) wynosi 1 odrzucenie dwóch
ostatnich składników nie spowoduje błędu większego niż ich suma, wynosząca 4,206 10 5, co po dodaniu do
maksymalnego błędu obciętego rozwinięcia w szereg Taylora (2,76 10 7) da w sumie maksymalny błąd mniejszy niż
4,234 10 5  czyli jeszcze do przyjęcia. Możemy się więc ograniczyć do pierwszych trzech składników. Teraz
wielomiany T0, T2 i T4 pozostające w rozwinięciu zamienimy znów na potęgi x i otrzymamy:
2 4
c o s x 0,99995795 0,49924045x 0,03962674x
Jak widać konieczne było wykonanie 4 mnożeń i 2 dodawań. Gdyby oryginalne rozwinięcie obciąć po wyrazie z
x4 błąd wynosiłby 1/6! 1,39 10 3, czyli byłby ok. 28 razy większy niż założona wartość.
Wielomiany Newtona
Najpierw zajmiemy się punktami równoodległymi o wzajemnej odległości h. Różnicami skończonymi funkcji
f(x) będziemy nazywać wielkości:
"yi = yi+1  yi  różnice pierwszego rzędu
"2yi = "yi+1  "yi  różnice drugiego rzędu
...........
"kyi+1 = "k 1yi+1  "k 1yi  różnice k-tego rzędu
Różnice te można zapisać w postaci tablicy:
x y "y "2y "3y "4y
x0 y0 "y0 "2y0 "3y0 "4y0
x1 y1 "y1 "2y1 "3y1 "4y1
x2 y2 "y2 "2y2 "3y2
x3 y3 "y3 "2y3
x4 y4 "y4
x5 y5
Pierwszy wzór interpolacyjny Newtona
Wykorzystując tak zdefiniowane różnice możemy zapisać pierwszy wzór interpolacyjny Newtona:
q q 1 q q 1 q n 1
y x y0 q"y0 "2y0 "ny0
2! n!
x x0
gdzie q . Warto zauważyć, że wzór ten wykorzystuje górny wiersz tabeli powyżej.
h
q q 1 q n
n 1 n 1
Reszta dana jest wzorem: R x h f ¾
n
n 1 !
Aproksymacja / Interpolacja
Metody Numeryczne w Fizyce 1 Aproksymacja i Interpolacja 13
gdzie ¾ jest jakimÅ› punktem leżącym wewnÄ…trz najmniejszego przedziaÅ‚u zawierajÄ…cego wszystkie punkty xi oraz
punkt x. Jeśli znamy również punkt xn+1 możemy posłużyć się wygodniejszym wzorem:
"n 1y0
R x q q 1 q n
n
n 1 !
Wartość n należy dobrać tak aby kolejne różnice "nyi były praktycznie stałe.
Dla n = 1 i n = 2 mamy przypadki szczególne: interpolacji liniowej i interpolacji kwadratowej:
y x y0 q"y0
q q 1
y x y0 q"y0 "2y0
2!
Dla przykładu rozważmy interpolację funkcji ln(x), załóżmy, że potrzebna jest nam wartość ln(1001) podczas
gdy znamy z zadowalającą dokładnością wartości dla x = 1000, 1010, 1020, 1030, 1040, 1050, 1060, 1070 i 1080
(być może nie wykorzystamy ich wszystkich). Różnice kolejnych rzędów będą następujące:
xi yi "y "2y "3y "4y "5y "6y "7y "8y
1000 6.907755 0.0099503 -0.000098 0.00000191 -5.545e-08 2.12e-09 -1.0e-10 5.68e-12 -3.8e-13
1010 6.917706 0.0098523 -0.0000961 0.00000186 -5.333e-08 2.02e-09 -9.5e-11 5.30e-12
1020 6.927558 0.0097562 -0.0000943 0.0000018 -5.130e-08 1.93e-09 -9.0e-11
1030 6.937314 0.0096619 -0.0000925 0.00000175 -4.938e-08 1.84e-09
1040 6.946976 0.0095695 -0.0000907 0.0000017 -4.754e-08
1050 6.956545 0.0094787 -0.000089 0.00000166
1060 6.966024 0.0093897 -0.0000873
1070 6.975414 0.0093024
1080 6.984716
Wykorzystując te dane otrzymamy następujące kolejne przybliżenia wartości ln(1001):
reszta faktyczny
kolejne przyrost
y(x) y(x)-y0 obliczona błąd
składniki: względny
Rn względny
y0 = 6.907755 6.9077553 0
1 = 0.000985 6.9087405 0.00098523 0.00014261 0.00000147 0.000002
2 = 0.000004 6.9087448 0.00098947 0.00000061 0.00000005 0.000001
3 = 5.45e-08 6.9087448 0.00098953 7.8912e-09 4.5828e-09 0.000001
4 = 1.15e-09 6.9087448 0.00098953 1.6583e-10 6.8416e-10 0.000001
5 = 3.42e-11 6.9087448 0.00098953 4.9514e-12 1.5890e-10 0.000001
6 = 1.32e-12 6.9087448 0.00098953 1.9168e-13 5.2896e-11 0.000001
7 = 6.30e-14 6.9087448 0.00098953 9.1277e-15 2.4386e-11 0.000001
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
14 Metody Numeryczne w Fizyce Andrzej Brozi
Decyzję o tym, w którym momencie (tzn. po uwz-
ględnieniu ilu składników) zaprzestać sumowania
podejmiemy na podstawie obliczanej na bieżąco reszty.
Jak widać rozsądną dokładność uzyskujemy już po 3-
-cim składniku.
Warto zwrócić uwagę, że drugi z wcześniej
podanych wzorów na resztę sumowania ma identyczną
postać jak wzory na kolejne składniki, można więc bez
dodatkowego nakładu obliczeń zadecydować o
momencie ich przerwania.
Drugi wzór interpolacyjny Newtona
W drugim wzorze Newtona wykorzystuje się dolne wartości różnic zapisanych w tabeli:
q q 1 q q 1 q n 1
y x yn q"yn 1 "2yn 2 "ny0
2! n!
Reszta jest postaci:
q q 1 q n
n 1 n 1
R x h f ¾
n
n 1 !
Oznaczenia q i ¾ jak poprzednio.
Wzory Gaussa, Stirlinga i Bessela
Różnice można zapisać także w inny sposób:
x y "y "2y "3y "4y "5y "6y
x y 4
 4
"y 4
x y 3 "2y 4
 3
"y 3 "3y 4
x y 2 "2y 3 "4y 4
 2
"y 2 "3y 3 "5y 4
x y 1 "2y 2 "4y 3 "6y 4
 1
"y 1 "3y 2 "5y 3
x0 y0 "2y 1 "4y 2 "6y 3
"y0 "3y 1 "5y 2
x1 y1 "2y0 "4y 1 "6y 2
"y1 "3y0 "5y 1
x2 y2 "2y1 "4y0
"y2 "3y1
x3 y3 "2y2
"y3
x4 y4
Aproksymacja / Interpolacja
Metody Numeryczne w Fizyce 1 Aproksymacja i Interpolacja 15
Pierwszy wzór Gaussa
q q 1 q 1 q q 1 q 1 q q 1 q 2
y x y0 q"y0 "2y 1 "3y 1 "4y 2
2! 3! 4!
q 2 q 1 q q 1 q2 q n 1 q n 1 q n 1 q n
"5y 2 "2n 1y n 1 "2ny n
5! 2n 1 ! 2n !
Pierwszy wzór interpolacyjny Gaussa (wykorzystujący powyższe różnice) ma postać:
gdzie q jest określone tak samo jak w przypadku wzorów Newtona. Różnice wykorzystywane w tym wzorze
zaznaczono w tabeli pogrubieniem.
Drugi wzór Gaussa
q 1 q q 1 q q 1 q 2 q 1 q q 1
y x y0 q"y 1 "2y 1 "3y 2 "4y 2
2! 3! 4!
q n 1 q n 1 q n q n 1 q n 1
2n 1
" y n "2ny n
2n 1 ! 2n !
Wzór na resztę jest wspólny dla obu wzorów Gaussa:
2n 1 2n 1
h f ¾
2 2 2 2
R2n q q 12 q 22 q n
2n 1 !
Wzór Stirlinga
Wzór Stirlinga jest w gruncie rzeczy średnią arytmetyczną obu wzorów Gaussa:
2 2 2
"y 1 "y0 q 2
q q 12 "3y "3y q q 12
2 1
y x y0 q "2y 1 "4y 2
2 2 3! 2 4!
2 2 2 2
q q 12 q 22 "5y "5y q q 12 q 22
3 2
"6y 3
5! 2 6!
2 2 2
"2n 1y n "2n 1y n 1 q 2 q 2 12 q 2 n 1 2
q q 12 q n 1
"2ny n
2n 1 ! 2 2n !
Wzór na resztę jest taki sam jak dla wzorów Gaussa.
Wzór Bessela
y0 y 1 "2y 1 "2y0
q q 1
y x q 0,5 "y0
2 2 2
"4y 2 "4y 1
q 0,5 q q 1 q q 1 q 1 q 2
"3y 1
3! 4! 2
q 0,5 q q 1 q 1 q 2
"5y 2
5!
"6y 3 "6y 2
q q 1 q 1 q 2 q 2 q 3
6! 2
"2ny n "2ny n 1
q q 1 q 1 q 2 q 2 q n q n 1
2n ! 2
q 0,5 q q 1 q 1 q 2 q 2 q n q n 1
"2n 1y n
2n 1 !
Wzór na resztę:
2n 2
h
2n 2 2 2 2 2
R x f ¾ q q 12 q 22 q n q n 1
n
2n 2 !
Interpolacja w dwóch wymiarach
Omówienie interpolacji w dwóch wymiarach ograniczymy do przypadku, gdy znane wartości funkcji dwóch
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
16 Metody Numeryczne w Fizyce Andrzej Brozi
zmiennych rozmieszczone są w węzłach siatki prostokątnej, a przynajmniej gdy punkt (x,y), w którym poszukujemy
wartości funkcji, otoczony jest przez cztery punkty (xi,yj), w których wartości te znamy: xi < x < xi+1 oraz yj < y < yj+1.
Najprostszym sposobem jest interpolacja liniowa (i do tego tylko sposobu ograniczymy nasze rozważania).
Można wprowadzić wielkości t i u:
x xi y yj
t ; u
xi 1 xi yj 1 yj
wówczas wartość funkcji w interesującym nas punkcie można przedstawić w postaci:
f x,y 1 t 1 u f xi,yj t 1 u f xi 1,yj tuf xi 1,yj 1 1 t uf xi,yj 1
Często metoda taka daje zadowalające wyniki jednak jej podstawową wadą jest fakt, że podczas  przekraczania
granic pomiędzy prostokątami wyznaczanymi przez punkty (xi,yj) gradient interpolowanej funkcji zmienia się w
sposób nieciągły.
Aproksymacja / Interpolacja