Personal tools

MN09LAB

From Studia Informatyczne


Interpolacja wielomianowa

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Oglądaj wskazówki i rozwiązania
Ukryj wskazówki i rozwiązania

Ćwiczenie: Błąd interpolacji dla węzłów równoodległych

Wyprowadź wzór na błąd interpolacji w przypadku węzłów równoodległych, tzn. \displaystyle x_i = a + ih, \displaystyle i = 0,\ldots,n:

jeśli \displaystyle |f^{(n+1)}(x)| \leq M dla \displaystyle x\in [a, a+nh], to

\displaystyle |f(x) - w_n(x)| \leq M \frac{h^{n+1}}{4(n+1)}.

Wskazówka

Wystarczy zauważyć, że jeśli \displaystyle x leży pomiędzy \displaystyle x_j a \displaystyle x_{j+1}, to \displaystyle |(x-x_j)(x-x_{j+1})| \leq h^2/4 i skorzystać z twierdzenia o błedzie interpolacji z wykładu.

Przetestuj eksperymentalnie jakość tego oszacowania dla kilku funkcji, dla których znasz wartość \displaystyle M, np.

  • wielomianu stopnia \displaystyle n,
  • wielomianu stopnia \displaystyle n+1,
  • funkcji \displaystyle \sin(x) na \displaystyle [0,1],
  • funkcji \displaystyle \sin(x) na \displaystyle [0,n\pi],
  • funkcji \displaystyle e^x.

porównując faktyczny błąd w \displaystyle [a, a+nh] z błędem z powyższego oszacowania.

Wskazówka

"Faktyczną" wartość błędu możesz dobrze przybliżyć, wyznaczając wartość różnicy \displaystyle |f(x) - w_n(x)| w bardzo wielu punktach przedziału. Gdy punkty te są dostatecznie gęsto rozłożone w całym przedziale, a funkcje --- jak u nas --- są gładkie, taka procedura da dobre oszacowanie prawdziwego błędu.

Rozwiązanie

Wartość podanego w zadaniu oszacowania zazwyczaj trochę przeszacowuje faktyczną wartość błędu, co możemy zobaczyć na przykładach, w których zajęliśmy się wielomianem stopnia 7 (czy wyniki różnią się, gdy weźmiemy wielomian stopnia 6?):

n = 7; 
a = 0; b = 1; x = linspace(a,b,n+1); h = x(2)-x(1);
X = linspace(a,b,ceil(1000*(b-a)));

y = x.^n; M = 0;
p = polyfit(x,y,n);
[norm(polyval(p,X) - X.^n, inf), (M*h^(n+1))/(4*(n+1))]

y = x.^(n+1); M = gamma(n+2);
p = polyfit(x,y,n);
[norm(polyval(p,X) - X.^(n+1), inf), (M*h^(n+1))/(4*(n+1))]

y = sin(x); M = 1;
p = polyfit(x,y,n);
[norm(polyval(p,X) - sin(X), inf), (M*h^(n+1))/(4*(n+1))]

a = 0; b = n*pi; x = linspace(a,b,n+1); h = x(2)-x(1);
X = linspace(a,b,ceil(1000*(b-a)));

y = sin(x); M = 1;
p = polyfit(x,y,n);
[norm(polyval(p,X) - sin(X), inf), (M*h^(n+1))/(4*(n+1))]

y = x.^(n+1); M = gamma(n+2);
p = polyfit(x,y,n);
[norm(polyval(p,X) - X.^(n+1), inf), (M*h^(n+1))/(4*(n+1))]

Oto nasze wyniki:

Dla \displaystyle x^n na \displaystyle [0,1]

2.2204e-16 0

Dla \displaystyle x^{n+1} na \displaystyle [0,1]

1.1112e-04 2.1857e-04

Dla \displaystyle \sin(x) na \displaystyle [0,\pi]

1.4338e-09 5.4208e-09

Dla \displaystyle \sin(x) na \displaystyle [0,n\,\pi]

1.00000 296.51659

Dla \displaystyle x^{n+1} na \displaystyle [0,n\,\pi]

6.0784e+06 1.1956e+07

Oczywiście, nie powinniśmy martwić się, że oszacowanie okazało się "fałszywe" dla wielomianu stopnia \displaystyle n, gdyż eksperymentalna wartość błędu wyniosła tylko \displaystyle \epsilon_{ \mbox{mach} } i była spowodowana oczywiście redukcją cyfr przy odejmowaniu (nasz wielomian interpolacyjny przecież pokrywa się z zadanym, więc wartość różnicy wartości wielomianów musi wyjść (prawie) zero.

Ćwiczenie: Aproksymacja funkcji sinus

Ile węzłów interpolacji wielomianowej Lagrange'a wystarczy, by przybliżyć funkcję \displaystyle f(x) = \sin(x) z błędem bezwzględnym \displaystyle 10^{-8} na całym przedziale \displaystyle [0,\frac{\pi}{4}]. Podaj odpowiedź w przypadku

  • węzłów równoodległych w \displaystyle [0,\frac{\pi}{4}],
  • węzłów Czebyszewa w \displaystyle [0,\frac{\pi}{4}].

Sprawdź eksperymentalnie, czy się nie pomyliłeś.

Rozwiązanie

Dla węzłów równoodległych mamy z poprzedniego zadania

\displaystyle |\sin(x) - w_n(x)| \leq \frac{h^{n+1}}{4(n+1)},

gdzie \displaystyle h = \frac{\pi}{4n}, skąd warunek na \displaystyle n:

\displaystyle \left(\frac{\pi}{4n}\right)^{n+1}\cdot \frac{1}{4(n+1)} \leq 10^{-8},

czyli \displaystyle n \geq 7. Rzeczywiście,

n = 6;

a = 0; b = pi/4; x = linspace(a,b,n+1); h = x(2)-x(1); X = linspace(a,b,ceil(1000*(b-a)));

y = sin(x); M = 1; p = polyfit(x,y,n); [norm(polyval(p,X) - sin(X), inf), (M*h^(n+1))/(4*(n+1))] ans =

       1.1720e-08   2.3519e-08

n = 7; a = 0; b = pi/4; x = linspace(a,b,n+1); h = x(2)-x(1); X = linspace(a,b,ceil(1000*(b-a)));

y = sin(x); M = 1; p = polyfit(x,y,n); [norm(polyval(p,X) - sin(X), inf), (M*h^(n+1))/(4*(n+1))] ans =

       1.6663e-10   7.8485e-10

Dla węzłów Czebyszewa musimy skorzystać z oszacowania wykładu. Rachunki pozostawiamy Czytelnikowi.

Ćwiczenie: Wagi wzoru barycentrycznego

Wyprowadź formułę na pierwszy wzór barycentryczny w przypadku \displaystyle n+1 węzłów równoodległych na odcinku \displaystyle [0,1],

\displaystyle  w_j = (-1)^{n-j}\frac{1}{h^{n} \, j! \, (n-j)!}.

Wywnioskuj stąd wzór na wagi w drugim wzorze barycentrycznym.

Wskazówka

Aby uprościć wagi w drugim wzorze barycentrycznym, zauważ, że skalowanie wszystkich wag przez tę samą stałą nie wpływa na jego wartość.

Podaj algorytm wyznaczania tych wag dla zadanego \displaystyle n.

Rozwiązanie

Wyprowadzenie wzoru dla odcinka \displaystyle [0,1] jest trywialne. Aby zaś przejść od odcinka \displaystyle [0,1] do dowolnego \displaystyle [a,b] należy dokonać transformacji liniowej węzłów: przesunięcie węzłów nic nie zmienia, bo we wzorze na wagi występują jedynie odległości między węzłami; natomiast skalowanie o czynnik \displaystyle (b-a) spowoduje zmianę wszystkich wag o wspólny czynnik \displaystyle (b-a)^{-n}, który oczywiście upraszcza się w drugim wzorze barycentrycznym.

Dlatego niezłym wzorem na wagi węzłów równoodległych w drugim wzorze barycentrycznym byłby

\displaystyle  w_j = (-1)^{n-j}\frac{1}{j! \, (n-j)!},

ale jeszcze lepszym ---

\displaystyle  w_j = (-1)^{n-j}\frac{n!}{j! \, (n-j)!} = (-1)^j \begin{pmatrix}  n \\ j \end{pmatrix} ,

gdyż wprowadza współczynniki dwumianu Newtona, które łatwo obliczać.

Ćwiczenie: Implementacja algorytmu barycentrycznego

Zaimplementuj drugi wzór barycentryczny tak, by w efekcie dostać dwie procedury:

  • c = polyfitb(x,y) --- wyznaczającą współczynniki \displaystyle c wielomianu interpolującego wartości \displaystyle y w węzłach równoodległych \displaystyle x;
  • Y = polyvalb(X,c) --- wyznaczającą, w zadanych węzłach \displaystyle X, wartości \displaystyle Y wielomianu interpolacyjnego o współczynnikach \displaystyle c.

Oszacuj koszt twojego algorytmu i porównaj z kosztem algorytmu różnic dzielonych (dla współczynników) i algorytmem Hornera (dla wartości) w dwóch przypadkach:

  • kiedy jest dużo węzłów interpolacji, a potrzebna jest tylko jedna wartość wielomianu interpolacyjnego;
  • na odwrót, kiedy jest mało węzłów interpolacji, a potrzeba wyznaczyć bardzo dużo wartości wielomianu interpolacyjnego.
Rozwiązanie

Musimy wyznaczyć współczynniki wagowe \displaystyle w_j. Jeśli z góry wiemy, że np. stopień naszego wielomianu interpolacyjnego nie przekroczy, powiedzmy, \displaystyle N=20, to możemy przyspieszyć naszą procedurę, wykonując precomputing: wszak wagi w drugim wzorze barycentrycznym są wyznaczone niezależnie od konkretnych węzłów!

Zatem musimy wyliczyć (całkowite!) wartości trójkątnej tabelki

\displaystyle N=0 1
\displaystyle N=1 1 -1
\displaystyle N=2 1 -2 1
\displaystyle N=3 1 -3 3 1
\displaystyle \vdots ...
\displaystyle N=20 1 ...

a potem tylko sięgnąć do odpowiedniego wiersza.

Ćwiczenie: Dodawanie węzła interpolacji

Często w aplikacjach takich jak programy CAD zachodzi potrzeba dodania na bieżąco dodatkowego węzła interpolacji. Podaj, jak to zrobić --- i jaki to będzie miało koszt --- gdy interpolant jest zadany w bazie

  • naturalnej
  • Newtona
  • Lagrange'a
Rozwiązanie

Po raz kolejny okazuje się, jak niedobra jest baza naturalna: koszt wyznaczenia nowych współczynników to koszt rozwiązania układu z macierzą gęstą od nowa. Fakt, dodajemy do niej tylko dodatkowy wiersz i dodatkową kolumnę, więc przy odrobinie szczęścia możemy wykorzystać czynniki rozkładu starej (mniejszej) macierzy, ale i tak, koszt będzie co najmniej \displaystyle O(N^2).

Dla bazy Newtona wystarczy doliczyć jeszcze jeden wiersz w tabelce różnic dzielonych (koszt \displaystyle O(N)), ale pamiętajmy, że zapewne zmienimy w ten sposób uporządkowanie węzłów, co może trochę pogroszyć numeryczne własności tak otrzymanego interpolantu.

Dla wzoru barycentrycznego, musimy po prostu aktualizować wszystkie wagi: także kosztem \displaystyle O(N).

Ćwiczenie: Algorytm Hornera może więcej

Pokaż, że algorytm Hornera obliczania wartości \displaystyle w(\xi) wielomianu danego w postaci potęgowej jest jednocześnie algorytmem dzielenia tego wielomianu przez jednomian \displaystyle (x-\xi). Dokładniej, jeśli \displaystyle w(x)=\sum_{j=0}^n a_jx^j to

\displaystyle w(x)\,=\,\Big(\sum_{j=1}^n v_jx^{j-1}\Big) \cdot (x-\xi)\,+\,v_0,

gdzie \displaystyle v_j są zdefiniowane tak jak w algorytmie Hornera. Zaimplementuj i sprawdź na przykładach.

Rozwiązanie

Zobacz rozdział 3.5 w

  • D. Kincaid, W. Cheney, Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa, 2006.


Ćwiczenie

Niech dane będą \displaystyle \epsilon>0 i funkcja \displaystyle f:[a,b]\toR, która jest nieskończenie wiele razy różniczkowalna i wszystkie jej pochodne są na \displaystyle [a,b] ograniczone przez \displaystyle M. Napisz program znajdujący stopień \displaystyle n oraz współczynniki w bazie Newtona wielomianu \displaystyle w_{f,n}\in\Pi_n interpolującego \displaystyle f z błędem

\displaystyle \|f-w_{f,n}\|_{C([a,b])}\,\le\,\epsilon.

Rozpatrz dwa przypadki: gdy węzły interpolacji są równoodległe, oraz gdy węzły są czebyszewowskie.