MN09LAB

Interpolacja wielomianowa

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

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

Ć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. $x_i=a+ih$, $i=0,\dots ,n$:

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

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

Wskazówka

Wystarczy zauważyć, że jeśli $x$ leży pomiędzy $x_j$ a $x_{j+1}$, to $|(x-x_j)(x-x_{j+1})|\le 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ść $M$, np.

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

porównując faktyczny błąd w $[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 $|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 <math>\displaystyle x^n</math> na <math>\displaystyle [0,1]</math> 2.2204e-16 0 Dla <math>\displaystyle x^{n+1}</math> na <math>\displaystyle [0,1]</math> 1.1112e-04 2.1857e-04 Dla <math>\displaystyle \sin(x)</math> na <math>\displaystyle [0,\pi]</math> 1.4338e-09 5.4208e-09 Dla <math>\displaystyle \sin(x)</math> na <math>\displaystyle [0,n\,\pi]</math> 1.00000 296.51659 Dla <math>\displaystyle x^{n+1}</math> na <math>\displaystyle [0,n\,\pi]</math> 6.0784e+06 1.1956e+07

Oczywiście, nie powinniśmy martwić się, że oszacowanie okazało się "fałszywe" dla wielomianu stopnia $n$, gdyż eksperymentalna wartość błędu wyniosła tylko ${ϵ}_{\text{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ę $f(x)=\sin (x)$ z błędem bezwzględnym $10^{-8}$ na całym przedziale $[0,\frac{\pi }{4}]$. Podaj odpowiedź w przypadku

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

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

Rozwiązanie

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

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

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

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

czyli $n\ge 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 $n+1$ węzłów równoodległych na odcinku $[0,1]$,

$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 $n$.

Rozwiązanie

Wyprowadzenie wzoru dla odcinka $[0,1]$ jest trywialne. Aby zaś przejść od odcinka $[0,1]$ do dowolnego $[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 $(b-a)$ spowoduje zmianę wszystkich wag o wspólny czynnik $(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

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

ale jeszcze lepszym ---

$w_j=(-1{)}^{n-j}\frac{n!}{j!(n-j)!}=(-1{)}^j\left(\begin{array}{c}n\\ j\end{array}\right),$

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 $c$ wielomianu interpolującego wartości $y$ w węzłach równoodległych $x$;
• Y = polyvalb(X,c) --- wyznaczającą, w zadanych węzłach $X$, wartości $Y$ wielomianu interpolacyjnego o współczynnikach $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 $w_j$. Jeśli z góry wiemy, że np. stopień naszego wielomianu interpolacyjnego nie przekroczy, powiedzmy, $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

$N=0$	1				\
$N=1$	1	-1			\
$N=2$	1	-2	1		\
$N=3$	1	-3	3	1	\
\vdots			...		\
$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 $O(N^2)$.

Dla bazy Newtona wystarczy doliczyć jeszcze jeden wiersz w tabelce różnic dzielonych (koszt $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 $O(N)$.

Ćwiczenie: Algorytm Hornera może więcej

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

$w(x)=\left(\sum _{j=1}^n{v}_j{x}^{j-1}\right)\cdot (x-\xi )+v_0,$

gdzie $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ą $ϵ>0$ i funkcja Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle f:[a,b]\toR} , która jest nieskończenie wiele razy różniczkowalna i wszystkie jej pochodne są na $[a,b]$ ograniczone przez $M$. Napisz program znajdujący stopień $n$ oraz współczynniki w bazie Newtona wielomianu $w_{f,n}\in {\mathrm{\Pi }}_n$ interpolującego $f$ z błędem

$\Vert f-w_{f,n}{\Vert }_{C([a,b])}\le ϵ.$

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