1 Wstęp do metod numerycznych Wykład 6 Interpolacja 1 dr inż. Wojciech Bieniecki Instytut Matematyki i Informatyki http://wbieniec.kis.p.lodz.pl/pwsz

2 Sformułowanie zagadnienia interpolacji 2 Wiele zjawisk fizycznych jest opisywanych przez funkcje, których postaci nie znamy, ale potrafimy obliczyć lub zmierzyć wartości tych funkcji oraz ich pochodnych dla określonych wartości argumentu. Na przykład, w określonych chwilach czasowych możemy mierzyć niektóre parametry sygnałów elektrycznych. Interpolacją nazywamy postępowanie prowadzące do znalezienia wartości pewnej funkcji f(x) w dowolnym punkcie przedziału na podstawie znanych wartości tej funkcji w punktach x 0, x 1,..., x n, zwanych węzłami interpolacji (x 0 < x 1

3 Sformułowanie zagadnienia interpolacji 3. Interpolację lub ekstrapolację stosujemy najczęściej w następujących przypadkach: b) gdy obliczanie wartości pewnej funkcji F(x) bezpośrednio z określającego ją wzoru nastręcza zbyt duże trudności rachunkowe. a) gdy nie znamy samej funkcji f(x), a tylko jej wartości w pewnych punktach (tak przeważnie bywa w naukach doświadczalnych); Wtedy zastępujemy ją prostszą funkcją f(x), o której zakładamy, że w punktach x 0, x 1,..., x n ma te same wartości co funkcja F(x). W tym przypadku F(x) nazywamy funkcją interpolowaną, zaś f(x) funkcją interpolującą (rys. poniżej).

4 Sformułowanie zagadnienia interpolacji 4

5 5 Funkcji interpolującej poszukuje się zwykle w pewnej określonej postaci. Najczęściej zakłada się, że jest to wielomian lub funkcja wymierna. Będziemy rozważać interpolację -wielomianami algebraicznymi, -wielomianami trygonometrycznymi -funkcjami sklejanymi. Zastosowanie do obliczeń bardzo szybkich komputerów o dużych pamięciach spowodowało, że interpolacja, będąca niegdyś jedną z podstawowych metod numerycznych, straciła nieco na znaczeniu. Obecnie stosuje się albo proste metody, jak interpolacja liniowa czy kwadratowa, albo też bardziej złożone, wymagające użycia komputera, jak np. interpolacja za pomocą funkcji sklejanych. Wzory interpolacyjne są punktem wyjścia do wyprowadzenia wielu metod stosowanych w różnych działach metod numerycznych (rozwiązywanie równań, różniczkowanie i całkowanie numeryczne, numeryczne rozwiązywanie równań różniczkowych zwyczajnych)

6 Interpolacja liniowa 6 Interpolacja zerowego rzędu polega na przybliżeniu wykresu funkcji pomiędzy węzłami liniami poziomymi. Otrzymana funkcja interpolująca jest nieciągła Interpolacja liniowa polega na przybliżeniu wykresu funkcji pomiędzy węzłami linią prostą

8 Przykład w Pythonie 8 >>> import numpy as np. >>> import matplotlib.pyplot as plt >>> xp = [1,2,3] >>> fp = [3,2,0] >>> x = 2.5 >>> y = np.interp(x, xp, fp) >>> y 1.0 >>> plt.plot(xp,fp,'-') [ ] >>> plt.plot(x,y,'o') [ ] >>> plt.show() >>> import numpy as np. >>> import matplotlib.pyplot as plt >>> xp = [1,2,3] >>> fp = [3,2,0] >>> x = 2.5 >>> y = np.interp(x, xp, fp) >>> y 1.0 >>> plt.plot(xp,fp,'-') [ ] >>> plt.plot(x,y,'o') [ ] >>> plt.show()

9 Interpolacja kwadratowa W tym sposobie interpolacji między węzłami interpolacji poszukujemy wielomianu drugiego stopnia —funkcji kwadratowej: F(x) = a 0 + a 1 x +a 2 x 2

10 Interpolacja kwadratowa Ponieważ funkcja kwadratowa jest jednoznacznie określona przez 3 punkty, więc interpolację przeprowadzamy dla przedziałów [x i-1, x i+1 ], stosując następujący zapis funkcji kwadratowej: F(x)= b 0 + b 1 (x – x i-1 ) + b 2 (x - x i-1 )(x - x i )

11 Interpolacja kwadratowa 11 Wartości współczynników b 0, b 1, b 2 wyznacza się na podstawie funkcji f(x i-1 ), f(x i ), f(x i+1 )

12 Interpolacja kwadratowa 12 Podstawiając wartości b 0 i b 1 mamy Porównując zapisy ze współczynnikami a i oraz b i, można na podstawie wartości współczynników b 0, b 1, b 2 obliczyć współczynniki a 0, a 1, a 2 ze wzorów:

13 Interpolacja kwadratowa w Pythonie 13 >>> import numpy as np >>> from scipy.interpolate import interp1d >>> import matplotlib.pyplot as plt >>> x = np.linspace(0, 10, num=8, endpoint=True) >>> y = np.cos(-x**2/9.0) >>> f = interp1d(x, y, kind='cubic') >>> xnew = np.linspace(0, 10, num=51, endpoint=True) >>> plt.plot(x, y, 'o', xnew, f(xnew), '-',xnew, np.cos(-xnew**2/9.0),'b') >>> plt.legend(['dane','interpolacja','funkcja'],loc='best') >>> plt.show() >>> >>> import numpy as np >>> from scipy.interpolate import interp1d >>> import matplotlib.pyplot as plt >>> x = np.linspace(0, 10, num=8, endpoint=True) >>> y = np.cos(-x**2/9.0) >>> f = interp1d(x, y, kind='cubic') >>> xnew = np.linspace(0, 10, num=51, endpoint=True) >>> plt.plot(x, y, 'o', xnew, f(xnew), '-',xnew, np.cos(-xnew**2/9.0),'b') >>> plt.legend(['dane','interpolacja','funkcja'],loc='best') >>> plt.show() >>>

14 Interpolacja wielomianowa 14 Wielomianem interpolacyjnym W n (x) nazywamy wielomian stopnia co najwyżej n, który w punktach x 0, x 1,..., x n przyjmuje wartości y 0, y 1,..., y n. Twierdzenie Istnieje dokładnie jeden wielomian interpolacyjny, który w punktach x 0, x 1,..., x n przyjmuje wartości y 0, y 1,..., y n. Dowód Węzły interpolacji x 0, x 1,..., x n mogą być rozmieszczone w zupełnie dowolny sposób. Szukany wielomian zapisujemy w postaci: W n (x) = a 0 + a 1 x + a 2 x 2 +... + a n x n Korzystając z definicji wielomianu interpolacyjnego otrzymujemy układ n+1 równań z n+1 niewiadomymi współczynnikami a 0, a 1,..., a n : a 0 + a 1 x 0 + a 2 x 0 2 +... + a n x 0 n = y 0, a 0 + a 1 x 1 + a 2 x 1 2 +... + a n x 1 n = y 1,............................. a 0 + a 1 x n + a 2 x n 2 +... + a n x n n = y n.

15 Interpolacja wielomianowa 15 Macierz współczynników tego układu jest macierzą Vandermonde’a postaci zaś wyznacznik D macierzy A ma postać Przy założeniach, że x i  x j dla i  j, mamy zawsze Zatem układ (3.2) ma dokładnie jedno rozwiązanie, a wartości a i według twierdzenia Cramera są określone wzorem gdzie D ij (j = 0, 1, 2,..., n) są kolejnymi dopełnieniami algebraicznymi elementów i-tej kolumny wyznacznika D.

16 Wielomian interpolacyjny Lagrange’a 16 Francuski matematyk, astronom i fizyk. Genialny samouk, profesor Szkoły Artyleryjskiej w Turynie, członek a potem prezes (1766- 1787) Akademii Nauk w Berlinie. Od 1772 członek francuskiej Akademii Nauk. W dziele "Mechanika analityczna" (1788) usystematyzował dotychczasową wiedzę z zakresu mechaniki teoretycznej, wprowadzając wiele własnych koncepcji, m.in. sformułował Lagrange`a równania mechaniki, wprowadził tzw. funkcję Lagrange`a, rozwinął teorię mechaniki nieba, podał metodę wyznaczania orbit. W optyce zaproponował wzór wyrażający powiększenie układu optycznego za pomocą odpowiedniej aparatury wejściowej i wyjściowej (tzw. wzór Lagrange`a-Helmholza), sformułował równanie ruchu płynu, matematyczne twierdzenie o wartości średniej oraz zajmował się równaniami różniczkowymi (Lagrange`a równanie różniczkowe). W okresie rewolucji francuskiej zaangażowany w proces modernizacji miar i wag. de LAGRANGE Joseph Louis (1736-1813)

17 Wielomian interpolacyjny Lagrange’a 17 Twierdzenie Wielomian W n (x) postaci poniższego wzoru jest wielomianem interpolacyjnym dla dowolnego wyboru n+1 węzłów interpolacji x 0, x 1,..., x n.

18 Wielomian interpolacyjny Lagrange’a 18 Przyjmując oznaczenia:  k (x) = (x - x 0 )(x - x 1 )...(x - x n ), - wartość pochodnej wielomianu  n (x) w punkcie x j będącym zerem tego wielomianu, Wzór ten nazywamy wzorem interpolacyjnym Lagrange’a opartym na węzłach x 0, x 1,..., x n. można wielomian interpolacyjny W n (x) zapisać w postaci:

19 Przykład 19 Znaleźć wielomian interpolacyjny, który w punktach -2, -1, 0, 2 przyjmuje odpowiednio wartości 0, 1, 1, 2. Stosując wzór Lagrange’a dla n = 3 otrzymujemy:

20 Wielomian Lagrange’a w Pythonie 20 >>> import numpy as np >>> from scipy.interpolate import lagrange >>> x = [-2,0,1,2] >>> w = [0, 1, 1, 2] >>> W =lagrange(x,w) >>> W poly1d([ 0.16666667, 0., -0.16666667, 1.]) >>> W(0.5) 0.9375 >>> >>> import numpy as np >>> from scipy.interpolate import lagrange >>> x = [-2,0,1,2] >>> w = [0, 1, 1, 2] >>> W =lagrange(x,w) >>> W poly1d([ 0.16666667, 0., -0.16666667, 1.]) >>> W(0.5) 0.9375 >>>

21 Wzór interpolacyjny Newtona 21 Urodzony w Anglii, studiował na Uniwersytecie w Cambridge, w wieku 21 lat rozpoczął własne badania, które zrewolucjonizowały świat. Pierwszą opublikowaną pracą Newtona były badania nad naturą światła (prawa odbicia i załamania), dzięki którym zaprojektował i zbudował pierwszy teleskop zwierciadlany. W wieku 29 lat odkrycia te oraz wyniki wielu innych doświadczeń z dziedziny optyki przedstawił Brytyjskiemu Towarzystwu Królewskiemu (British Royal Society). Prace w dziedzinie optyki mają jednakże mniejsze znaczenie niż osiągnięcia w dziedzinie czystej matematyki i mechaniki. Jego największym wkładem do matematyki jest odkrycie rachunku całkowego, w dziedzinie mechaniki - poszerzenie nauki o ruchu ciał. W 1687 r. opublikował dzieło: Matematyczne zasady filozofii przyrody (Philosophiae naturalis principia mathematica), w którym przedstawił prawo ciążenia i prawa ruchu, pokazując również, jak można wykorzystać te prawa do precyzyjnego określenia ruchu planet dookoła Słońca. Wniósł znaczny wkład do termodynamiki i akustyki, ogłosił zasady zachowania pędu i momentu pędu, odkrył słynny dwumian Newtona oraz podał pierwsze przekonujące wyjaśnienie pochodzenia gwiazd. NEWTON Isaac (1642 - 1727)

22 Wzór interpolacyjny Newtona 22 Wzór interpolacyjny Lagrange’a jest niewygodny do stosowania w przypadku, gdy zmienia się liczba węzłów. Wtedy do wyznaczenia wielomianu określonego stopnia trzeba powtarzać obliczenia od początku. Zatem poprzez dodanie nowych węzłów interpolacyjnych nie można modyfikować wcześniej wyznaczonego wielomianu Lagrange’a. Wzór interpolacyjny Newtona równoważny wielomianowi Lagrange’a usuwa tę niedogodność.

23 Wzór interpolacyjny Newtona 23 Niech x 0, x 1,..., x n, będą węzłami interpolacji, w których wielomian interpolacyjny przyjmuje odpowiednio wartości y 0, y 1,..., y n. Można wówczas zdefiniować wyrażenia zwane ilorazami różnicowymi: - pierwszego rzędu:- drugiego rzędu (analogicznie): - k-tego rzędu: dla k = 1, 2,... oraz i = 0, 1, 2,....

24 Wzór interpolacyjny Newtona 24 Z ilorazów różnicowych tworzy się zazwyczaj tablicę xixi yiyi rzędu 1rzędu 2rzędu 3rzędu 4rzędu 5 x 0 x 1 x 2 x 3 x 4 x 5... y 0 y 1 y 2 y 3 y 4 y 5... [x0,x1][x1,x2][x2,x3][x3,x4][x4,x5][x0,x1][x1,x2][x2,x3][x3,x4][x4,x5] [x0,x1,x2][x1,x2,x3][x2,x3,x4][x3,x4,x5][x0,x1,x2][x1,x2,x3][x2,x3,x4][x3,x4,x5] [x0,x1,x2,x3][x1,x2,x3,x4][x2,x3,x4,x5][x0,x1,x2,x3][x1,x2,x3,x4][x2,x3,x4,x5] [x0,x1,x2,x3,x4][x1,x2,x3,x4,x5][x0,x1,x2,x3,x4][x1,x2,x3,x4,x5]

25 Wzór interpolacyjny Newtona 25 Korzystając z definicji ilorazów różnicowych, otrzymujemy wzór interpolacyjny Newtona z ilorazami różnicowymi: Wzór interpolacyjny Newtona ma tę własność, że rozszerzenie zadania interpolacji o dodatkowy węzeł sprowadza się do obliczenia i dołączenia do poprzednio wyznaczonego wielomianu dodatkowego składnika.

26 Przykład 26 Znaleźć wielomian interpolacyjny mając dane węzły (−1, −4), (0, −1), (1, 0), (2, 5) xixi yiyi rzędu 1rzędu 2rzędu 3 x 0 = -1 x 1 = 0 x 2 = 1 x 3 = 2 y 0 = -4 y 1 = -1 y 2 = 0 y 3 = 5 Wobec tego szukany wielomian interpolacyjny W(x) ma postać W(x) = −4 + 3(x + 1) − x(x + 1) + x(x + 1)(x − 1) = x 3 − x 2 + x − 1

27 Interpolacja Lagrange’a dla węzłów równoodległych 27 Załóżmy że węzły xi (i=0,1, … n) są równoodległe, czyli xi =x0+ih, to otrzymamy Przy węzłach równoodległych szczególnie wygodna jest postać wielomianu interpolacyjnego, którego współczynniki są wyrażone za pomocą tzw. różnic skończonych funkcji f. Wprowadźmy teraz odpowiednie pojęcia. Różnica zwykła (progresywna) funkcji f to operacja Δ, której wartością jest Δf zdefiniowana Operacja przesunięcia funkcji f Różnica wsteczna funkcji f Różnica centralna funkcji f

28 Interpolacja Lagrange’a (węzły równooległe) 28 Odpowiednie różnice rzędu n określamy jako różnice z różnic rzędu n - 1, np Korzystając z pojęcia różnicy progresywnej i różnicy wstecznej można udowodnić gdzie

29 Przykład 29 Dana jest następująca tablica wartości funkcji: Należy obliczyć wartość wielomianu interpolacyjnego w punkcie x = 3,58. Ponieważ węzły są równoodlegle, możemy użyć powyższego wzoru Sporządzamy najpierw tablice różnic progresywnych

30 Przykład 30 gdzie

31 Interpolacja Czebyszewa 31 Załóżmy, że wartości argumentów funkcji interpolowanej mieszczą się w przedziale [-1, 1]. Postulat ten nie ogranicza możliwości wykorzystania wzoru interpolacyjnego Czebyszewa, ponieważ dowolny przedział argumentów [a, b] poprzez podstawienie normalizuje się do przedziału [-1, 1]. Przykład: Zbiór węzłów: po podstawieniu: sprowadza się do zbioru:

32 Interpolacja Czebyszewa 32 Funkcje bazowe (tzw.bazę Czebyszewa) stanowi zbiór wielomianów określonych wzorem rekurencyjnym: Kilka pierwszych funkcji bazowych określonych wzorem rekurencyjnym wygląda następująco: Współczynniki wzoru interpolacyjnego Czebyszewa wynikają z układu równań:

33 Interpolacja Czebyszewa 33 Przykład: Dla zbioru punktów (węzłów) : (-0.5, 0.25), (0, 0), (1, 1) wyznaczyć wielomian interpolacyjny Czebyszewa. Jest to wielomian stopnia drugiego w postaci: czyli Układ równań sprowadza się do następującego Skąd otrzymujemy

34 Interpolacja Czebyszewa 34 UWAGA ! Przy dowolnym (różne odległości) doborze węzłów błąd zaokrągleń związanych z procedurą odwracania macierzy X jest istotnie mniejszy niż w przypadku interpolacji wielomianami w postaci naturalnej. Przyjęcie siatki węzłów wynikającej z zależności oraz nieco zmodyfikowanej bazy

35 Interpolacja Czebyszewa 35 prowadzi do macierzy X, dla której macierz odwrotną można obliczyć w bardzo prosty sposób Wielomianem interpolacyjnym w tym przypadku będzie więc suma w postaci: czyli: zawierająca n+1 nieznanych parametrów

36 Interpolacja Czebyszewa 36 Przykład: W przedziale [-1, 1] wybieramy cztery węzły interpolacji: x0, x1, x2, x3 ( n = 3 ) w ten sposób, że Otrzymujemy: x 0 = 0.924, x 1 = 0.383, x 2 = -0.383, x 3 = − 0.924 Załóżmy, że wartości funkcji f(x) w tych węzłach wynoszą: y0 = 2.224, y1 = 1.701, y2 = 3.885, y3 = 6.19 Bazę interpolacji stanowi zbiór wielomianów:

37 Interpolacja Czebyszewa 37 Konstruujemy i odwracamy macierz X skąd po obliczeniu otrzymujemy wartości współczynników: Postać wielomianu jest następująca: Po uproszczeniu otrzymujemy:

38 Interpolacja trygonometryczna 38

39 Interpolacja trygonometryczna 39

40 Zadanie interpolacji trygonometrycznej 40

41 Zadanie interpolacji trygonometrycznej 41

42 Zadanie interpolacji trygonometrycznej 42

43 Zadanie interpolacji trygonometrycznej 43

44 Przykład 44

45 Przykład 45

46 Przykład 46

47 Przykład 47

48 Funkcje sklejane 48

49 Funkcje sklejane 49

50 Funkcje sklejane 50

51 Interpolacyjne funkcje sklejane stopnia trzeciego 51

52 Interpolacyjne funkcje sklejane stopnia trzeciego 52

53 Interpolacyjne funkcje sklejane stopnia trzeciego 53

54 Interpolacyjne funkcje sklejane stopnia trzeciego 54

55 Interpolacyjne funkcje sklejane stopnia trzeciego 55

56 Interpolacyjne funkcje sklejane stopnia trzeciego 56

57 Interpolacyjne funkcje sklejane stopnia trzeciego 57

58 Interpolacyjne funkcje sklejane stopnia trzeciego 58

59 Przykład 59

60 Przykład 60 Razem mamy n*(m+1)=3*(3+1)=12 niewiadomych współczynników. Aby je wyznaczyć musimy ułożyć 12 równań. Wartości funkcji w węzłach interpolacji:

61 Przykład 61 Oraz dwa dodatkowe warunki albo na pierwszą albo na drugą pochodną na krańcach przedziału. Ciągłość funkcji i jej pochodnych aż do stopnia m-1=3-1 w wewnętrznych węzłach interpolacji: Mamy zatem tyle równań ile jest niewiadomych!

62 Przykład W postaci rozwiniętej powyższe warunki: 62 Dwa Dwa

63 Przykład 63 Dwa dodatkowe warunki na krańcach przedziału dla pochodnych: Po uwzględnieniu powyższych równań macierz współczynników przyjmie następującą postać:

64 Przykład 64

65 Przykład 65 Po zamianie wierszy macierz już nie posiada zer na głównej przekątnej:

66 Przykład 66 Wektor prawych stron po zmianie kolejności: 2.8345 -1.7804 1.0574 -0.1115 0.9189 0.1351 0.4189 -0.0405 -3.4605 2.7628 -0.1066 -0.0055 Wektor rozwiązań

67 Przykład 67 Ilustracja graficzna

68 Przykład 68 alpha=beha=0., π/4 or - π/4

69 Przykład interpolacji funkcją sklejaną 69