Aktualna wersja na dzień 21:49, 11 wrz 2023

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 + i h$ , $i = 0, \dots, n$ :

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

| f (x) - w_{n} (x) | \leq 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}) | \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ść $M$ , np.

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

porównując faktyczny błąd w $[a, a + n h]$ 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 $x^{n}$ na $[0, 1]$

2.2204e-16   0

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

1.1112e-04   2.1857e-04

Dla $\sin (x)$ na $[0, π]$

1.4338e-09   5.4208e-09

Dla $\sin (x)$ na $[0, n π]$

1.00000    296.51659

Dla $x^{n + 1}$ na $[0, n π]$

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 $ϵ_{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 $1 0^{- 8}$ na całym przedziale $[0, \frac{π}{4}]$ . Podaj odpowiedź w przypadku

węzłów równoodległych w $[0, \frac{π}{4}]$ ,
węzłów Czebyszewa w $[0, \frac{π}{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) | \leq \frac{h^{n + 1}}{4 (n + 1)}

,

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

{(\frac{π}{4 n})}^{n + 1} \cdot \frac{1}{4 (n + 1)} \leq 1 0^{- 8},

czyli $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 $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

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} (\begin{matrix} n \\ j \end{matrix}),

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
$⋮$			...
$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 (ξ)$ wielomianu danego w postaci potęgowej jest jednocześnie algorytmem dzielenia tego wielomianu przez jednomian $(x - ξ)$ . Dokładniej, jeśli $w (x) = \sum_{j = 0}^{n} a_{j} x^{j}$ to

w (x) = (\sum_{j = 1}^{n} v_{j} x^{j - 1}) \cdot (x - ξ) + v_{0}

,

gdzie $v_{j}$ są zdefiniowane tak jak w algorytmie Hornera. Zaimplementuj i sprawdź na przykładach.

Rozwiązanie

Ćwiczenie

Niech dane będą $ϵ > 0$ i funkcja $f : [a, b] \to R$ , 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 Π_{n}$ interpolującego $f$ z błędem

‖ f - w_{f, n} ‖_{C ([a, b])} \leq ϵ

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

@@ Linia 20: / Linia 20: @@
 <div class="exercise">
-Wyprowadź wzór na błąd interpolacji w przypadku węzłów równoodległych, tzn. <math>\displaystyle x_i = a + ih</math>, <math>\displaystyle i = 0,\ldots,n</math>:
+Wyprowadź wzór na błąd interpolacji w przypadku węzłów równoodległych, tzn. <math>x_i = a + ih</math>, <math>i = 0,\ldots,n</math>:
-jeśli <math>\displaystyle |f^{(n+1)}(x)| \leq M</math> dla <math>\displaystyle x\in [a, a+nh]</math>, to
+jeśli <math>|f^{(n+1)}(x)| \leq M</math> dla <math>x\in [a, a+nh]</math>, to
-<center><math>\displaystyle |f(x) - w_n(x)| \leq M \frac{h^{n+1}}{4(n+1)}.
+<center><math>|f(x) - w_n(x)| \leq M \frac{h^{n+1}}{4(n+1)}</math></center>
-</math></center>
 <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
-<div style="font-size:smaller; background-color:#f9fff9; padding: 1em"> Wystarczy zauważyć, że jeśli <math>\displaystyle x</math> leży pomiędzy <math>\displaystyle x_j</math> a <math>\displaystyle x_{j+1}</math>, to <math>\displaystyle |(x-x_j)(x-x_{j+1})| \leq h^2/4</math> i skorzystać z twierdzenia o błedzie interpolacji z wykładu. </div>
+<div style="font-size:smaller; background-color:#f9fff9; padding: 1em"> Wystarczy zauważyć, że jeśli <math>x</math> leży pomiędzy <math>x_j</math> a <math>x_{j+1}</math>, to <math>|(x-x_j)(x-x_{j+1})| \leq h^2/4</math> i skorzystać z twierdzenia o błedzie interpolacji z wykładu. </div>
 </div></div>
-Przetestuj eksperymentalnie jakość tego oszacowania dla kilku funkcji, dla których znasz wartość <math>\displaystyle M</math>, np.
+Przetestuj eksperymentalnie jakość tego oszacowania dla kilku funkcji, dla których znasz wartość <math>M</math>, np.
-* wielomianu stopnia <math>\displaystyle n</math>,
+* wielomianu stopnia <math>n</math>,
-* wielomianu stopnia <math>\displaystyle n+1</math>,
+* wielomianu stopnia <math>n+1</math>,
-* funkcji <math>\displaystyle \sin(x)</math> na <math>\displaystyle [0,1]</math>,
+* funkcji <math>\sin(x)</math> na <math>[0,1]</math>,
-* funkcji <math>\displaystyle \sin(x)</math> na <math>\displaystyle [0,n\pi]</math>,
+* funkcji <math>\sin(x)</math> na <math>[0,n\pi]</math>,
-* funkcji <math>\displaystyle e^x</math>.
+* funkcji <math>e^x</math>.
-porównując faktyczny błąd w <math>\displaystyle [a, a+nh]</math> z błędem z powyższego oszacowania.
+porównując faktyczny błąd w <math>[a, a+nh]</math> z błędem z powyższego oszacowania.
 <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
-<div style="font-size:smaller; background-color:#f9fff9; padding: 1em"> "Faktyczną" wartość błędu możesz dobrze przybliżyć, wyznaczając wartość różnicy <math>\displaystyle |f(x) - w_n(x)|</math> 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. </div>
+<div style="font-size:smaller; background-color:#f9fff9; padding: 1em"> "Faktyczną" wartość błędu możesz dobrze przybliżyć, wyznaczając wartość różnicy <math>|f(x) - w_n(x)|</math> 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. </div>
 </div></div>
@@ Linia 80: / Linia 79: @@
 Oto nasze wyniki:
-<div style="font-family: monospace; white-space: pre; border-style: dashed; border-width: thin; border-color: black; margin: 1em; padding:1em; color: #444; background-color:#fdfdfd;"><nowiki>Dla <math>\displaystyle x^n</math> na <math>\displaystyle [0,1]</math>
+Dla <math>x^n</math> na <math>[0,1]</math>
-.2204e-16   0
+<div style="font-family: monospace; white-space: pre; border-style: dashed; border-width: thin; border-color: black; margin: 1em; padding:1em; color: #444; background-color:#fdfdfd;"><nowiki>2.2204e-16   0
-Dla <math>\displaystyle x^{n+1}</math> na <math>\displaystyle [0,1]</math>
+</nowiki></div>
-.1112e-04   2.1857e-04
-Dla <math>\displaystyle \sin(x)</math> na <math>\displaystyle [0,\pi]</math>
+Dla <math>x^{n+1}</math> na <math>[0,1]</math>
-.4338e-09   5.4208e-09
+<div style="font-family: monospace; white-space: pre; border-style: dashed; border-width: thin; border-color: black; margin: 1em; padding:1em; color: #444; background-color:#fdfdfd;"><nowiki>1.1112e-04   2.1857e-04
-Dla <math>\displaystyle \sin(x)</math> na <math>\displaystyle [0,n\,\pi]</math>
+</nowiki></div>
-.00000    296.51659
-Dla <math>\displaystyle x^{n+1}</math> na <math>\displaystyle [0,n\,\pi]</math>
+Dla <math>\sin(x)</math> na <math>[0,\pi]</math>
-.0784e+06   1.1956e+07
+<div style="font-family: monospace; white-space: pre; border-style: dashed; border-width: thin; border-color: black; margin: 1em; padding:1em; color: #444; background-color:#fdfdfd;"><nowiki>1.4338e-09   5.4208e-09
+</nowiki></div>
+Dla <math>\sin(x)</math> na <math>[0,n\,\pi]</math>
+<div style="font-family: monospace; white-space: pre; border-style: dashed; border-width: thin; border-color: black; margin: 1em; padding:1em; color: #444; background-color:#fdfdfd;"><nowiki>1.00000    296.51659
+</nowiki></div>
+Dla <math>x^{n+1}</math> na <math>[0,n\,\pi]</math>
+<div style="font-family: monospace; white-space: pre; border-style: dashed; border-width: thin; border-color: black; margin: 1em; padding:1em; color: #444; background-color:#fdfdfd;"><nowiki>6.0784e+06   1.1956e+07
 </nowiki></div>
-Oczywiście, nie powinniśmy martwić się, że oszacowanie okazało się "fałszywe" dla wielomianu stopnia <math>\displaystyle n</math>, gdyż eksperymentalna wartość błędu wyniosła tylko <math>\displaystyle \epsilon_{ \mbox{mach} }</math> 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.
+Oczywiście, nie powinniśmy martwić się, że oszacowanie okazało się "fałszywe" dla wielomianu stopnia <math>n</math>, gdyż eksperymentalna wartość błędu wyniosła tylko <math>\epsilon_{ \mbox{mach} }</math> 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.
 </div></div></div>
@@ Linia 101: / Linia 108: @@
 Ile węzłów interpolacji wielomianowej Lagrange'a wystarczy, by przybliżyć
-funkcję <math>\displaystyle f(x) = \sin(x)</math> z błędem bezwzględnym <math>\displaystyle 10^{-8}</math> na całym przedziale <math>\displaystyle [0,\frac{\pi}{4}]</math>. Podaj odpowiedź w przypadku
+funkcję <math>f(x) = \sin(x)</math> z błędem bezwzględnym <math>10^{-8}</math> na całym przedziale <math>[0,\frac{\pi}{4}]</math>. Podaj odpowiedź w przypadku
-* węzłów równoodległych w  <math>\displaystyle [0,\frac{\pi}{4}]</math>,
+* węzłów równoodległych w  <math>[0,\frac{\pi}{4}]</math>,
-* węzłów Czebyszewa w  <math>\displaystyle [0,\frac{\pi}{4}]</math>.
+* węzłów Czebyszewa w  <math>[0,\frac{\pi}{4}]</math>.
 Sprawdź eksperymentalnie, czy się nie pomyliłeś.
@@ Linia 111: / Linia 118: @@
 Dla węzłów równoodległych mamy z poprzedniego zadania
-<center><math>\displaystyle |\sin(x) - w_n(x)| \leq \frac{h^{n+1}}{4(n+1)},
+<center><math>|\sin(x) - w_n(x)| \leq \frac{h^{n+1}}{4(n+1)}</math>,</center>
-</math></center>
-gdzie <math>\displaystyle h = \frac{\pi}{4n}</math>, skąd warunek na <math>\displaystyle n</math>:
+gdzie <math>h = \frac{\pi}{4n}</math>, skąd warunek na <math>n</math>:
-<center><math>\displaystyle \left(\frac{\pi}{4n}\right)^{n+1}\cdot \frac{1}{4(n+1)} \leq 10^{-8},
+<center><math>\left(\frac{\pi}{4n}\right)^{n+1}\cdot \frac{1}{4(n+1)} \leq 10^{-8},
 </math></center>
-czyli <math>\displaystyle n \geq 7</math>. Rzeczywiście,
+czyli <math>n \geq 7</math>. Rzeczywiście,
 <div style="font-family: monospace; white-space: pre; border-style: dashed; border-width: thin; border-color: black; margin: 1em; padding:1em; color: #444; background-color:#fdfdfd;"><nowiki>n = 6;
@@ Linia 151: / Linia 157: @@
 <div class="exercise">
-Wyprowadź formułę na pierwszy wzór barycentryczny w przypadku <math>\displaystyle n+1</math> węzłów równoodległych na odcinku <math>\displaystyle [0,1]</math>,
+Wyprowadź formułę na pierwszy wzór barycentryczny w przypadku <math>n+1</math> węzłów równoodległych na odcinku <math>[0,1]</math>,
-<center><math>\displaystyle
+<center><math>
 w_j = (-1)^{n-j}\frac{1}{h^{n} \, j! \, (n-j)!}.
 </math></center>
@@ Linia 162: / Linia 168: @@
 </div></div>
-Podaj algorytm wyznaczania tych wag dla zadanego <math>\displaystyle n</math>.
+Podaj algorytm wyznaczania tych wag dla zadanego <math>n</math>.
 </div></div>
 <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">
-Wyprowadzenie wzoru dla odcinka <math>\displaystyle [0,1]</math> jest trywialne. Aby zaś przejść od odcinka <math>\displaystyle [0,1]</math> do dowolnego <math>\displaystyle [a,b]</math> 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 <math>\displaystyle (b-a)</math> spowoduje zmianę wszystkich wag o wspólny czynnik <math>\displaystyle (b-a)^{-n}</math>, który oczywiście upraszcza się w drugim wzorze barycentrycznym.
+Wyprowadzenie wzoru dla odcinka <math>[0,1]</math> jest trywialne. Aby zaś przejść od odcinka <math>[0,1]</math> do dowolnego <math>[a,b]</math> 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 <math>(b-a)</math> spowoduje zmianę wszystkich wag o wspólny czynnik <math>(b-a)^{-n}</math>, 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
-<center><math>\displaystyle
+<center><math>
 w_j = (-1)^{n-j}\frac{1}{j! \, (n-j)!},
 </math></center>
 ale jeszcze lepszym ---
-<center><math>\displaystyle
+<center><math>
 w_j = (-1)^{n-j}\frac{n!}{j! \, (n-j)!} = (-1)^j \begin{pmatrix}
 n \\
@@ Linia 183: / Linia 189: @@
 gdyż wprowadza współczynniki dwumianu Newtona, które łatwo obliczać.
-<!--
-Narzucająca się implementacja, korzystająca z funkcji <code>silnia</code>, którą sami sobie radośnie napiszemy:
-{{algorytm|(niedobry!)|(niedobry!)|
-<pre>for (j=0; j <= n; j++)
-	w[j] = (pow(-1,n-j)*silnia(n)) / (silnia(j)*silnia(n-j));
-</pre>}}
-jest tak niedobra, że aż przyprawia o drżenie rąk, gdy ją wpisuję. Niestety, tylko w życiu prostolinijność jest cnotą. W numeryce, podobnie jak w całej informatyce, ceni się spryt.
--->
 </div></div></div>
@@ Linia 202: / Linia 196: @@
 Zaimplementuj drugi wzór barycentryczny tak, by w efekcie dostać dwie procedury:
-* <code>c = polyfitb(x,y)</code> --- wyznaczającą współczynniki <math>\displaystyle c</math> wielomianu interpolującego wartości <math>\displaystyle y</math> w węzłach <strong>równoodległych</strong> <math>\displaystyle x</math>;
+* <code>c = polyfitb(x,y)</code> --- wyznaczającą współczynniki <math>c</math> wielomianu interpolującego wartości <math>y</math> w węzłach <strong>równoodległych</strong> <math>x</math>;
-* <code>Y = polyvalb(X,c)</code> --- wyznaczającą, w zadanych węzłach <math>\displaystyle X</math>, wartości <math>\displaystyle Y</math> wielomianu interpolacyjnego o współczynnikach <math>\displaystyle c</math>.
+* <code>Y = polyvalb(X,c)</code> --- wyznaczającą, w zadanych węzłach <math>X</math>, wartości <math>Y</math> wielomianu interpolacyjnego o współczynnikach <math>c</math>.
 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:
@@ Linia 212: / Linia 206: @@
 <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">
-Musimy wyznaczyć współczynniki wagowe <math>\displaystyle w_j</math>. Jeśli z góry wiemy, że np. stopień naszego wielomianu interpolacyjnego nie przekroczy, powiedzmy, <math>\displaystyle N=20</math>, to możemy przyspieszyć naszą procedurę, wykonując <strong>precomputing</strong>: wszak wagi w drugim wzorze barycentrycznym są wyznaczone niezależnie od konkretnych węzłów!
+Musimy wyznaczyć współczynniki wagowe <math>w_j</math>. Jeśli z góry wiemy, że np. stopień naszego wielomianu interpolacyjnego nie przekroczy, powiedzmy, <math>N=20</math>, to możemy przyspieszyć naszą procedurę, wykonując <strong>precomputing</strong>: 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
@@ Linia 219: / Linia 213: @@
 |+ <span style="font-variant:small-caps"> </span>
 |-
-| <math>\displaystyle N=0</math>  ||  1  ||   ||   ||   ||  \
+| <math>N=0</math>  ||  1  ||   ||   ||   ||
 |-
-| <math>\displaystyle N=1</math>  ||  1  ||  -1  ||   ||   ||  \
+| <math>N=1</math>  ||  1  ||  -1  ||   ||   ||
 |-
-| <math>\displaystyle N=2</math>  ||  1  ||  -2  ||  1  ||   ||  \
+| <math>N=2</math>  ||  1  ||  -2  ||  1  ||   ||
 |-
-| <math>\displaystyle N=3</math>  ||  1  ||  -3  ||  3  ||  1  ||  \
+| <math>N=3</math>  ||  1  ||  -3  ||  3  ||  1  ||
 |-
-| \vdots  ||    ||    ||  ...  ||     ||  \
+| <math>\vdots</math>  ||    ||    ||  ...  ||     ||
 |-
-| <math>\displaystyle N=20</math>  ||  1  ||     ||  ...  ||    ||  \
+| <math>N=20</math>  ||  1  ||     ||  ...  ||    ||
 |}
@@ Linia 248: / Linia 242: @@
 <div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">
-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 <math>\displaystyle O(N^2)</math>.
+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 <math>O(N^2)</math>.
-Dla bazy Newtona wystarczy doliczyć jeszcze jeden wiersz w tabelce różnic dzielonych (koszt <math>\displaystyle O(N)</math>), 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 bazy Newtona wystarczy doliczyć jeszcze jeden wiersz w tabelce różnic dzielonych (koszt <math>O(N)</math>), 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 <math>\displaystyle O(N)</math>.
+Dla wzoru barycentrycznego, musimy po prostu aktualizować wszystkie wagi: także kosztem <math>O(N)</math>.
 </div></div></div>
@@ Linia 260: / Linia 254: @@
 <div class="exercise">
-Pokaż, że algorytm Hornera obliczania wartości <math>\displaystyle w(\xi)</math> wielomianu danego w postaci potęgowej jest jednocześnie algorytmem dzielenia tego wielomianu przez jednomian <math>\displaystyle (x-\xi)</math>. Dokładniej, jeśli <math>\displaystyle w(x)=\sum_{j=0}^n a_jx^j</math> to
+Pokaż, że algorytm Hornera obliczania wartości <math>w(\xi)</math> wielomianu danego w postaci potęgowej jest jednocześnie algorytmem dzielenia tego wielomianu przez jednomian <math>(x-\xi)</math>. Dokładniej, jeśli <math>w(x)=\sum_{j=0}^n a_jx^j</math> to
-<center><math>\displaystyle w(x)\,=\,\Big(\sum_{j=1}^n v_jx^{j-1}\Big) \cdot (x-\xi)\,+\,v_0,
+<center><math>w(x)\,=\,\Big(\sum_{j=1}^n v_jx^{j-1}\Big) \cdot (x-\xi)\,+\,v_0</math>,</center>
-</math></center>
-gdzie <math>\displaystyle v_j</math> są zdefiniowane tak jak w algorytmie Hornera. Zaimplementuj i sprawdź na przykładach.
+gdzie <math>v_j</math> są zdefiniowane tak jak w algorytmie Hornera. Zaimplementuj i sprawdź na przykładach.
 </div></div>
@@ Linia 279: / Linia 272: @@
 <div class="exercise">
-Pokaż numeryczną poprawność algorytmu Hornera obliczania wartości wielomianu ze względu na dane współczynniki <math>\displaystyle a_j</math> tego wielomianu.
+Pokaż numeryczną poprawność algorytmu Hornera obliczania wartości wielomianu ze względu na dane współczynniki <math>a_j</math> tego wielomianu.
 </div></div>
@@ Linia 288: / Linia 281: @@
 <div class="exercise">
-Niech dane będą <math>\displaystyle \epsilon>0</math> i funkcja <math>\displaystyle f:[a,b]\toR</math>, która jest nieskończenie wiele razy różniczkowalna i wszystkie jej pochodne są na <math>\displaystyle [a,b]</math> ograniczone przez <math>\displaystyle M</math>. Napisz program znajdujący stopień <math>\displaystyle n</math> oraz współczynniki w bazie Newtona wielomianu <math>\displaystyle w_{f,n}\in\Pi_n</math> interpolującego <math>\displaystyle f</math> z błędem
+Niech dane będą <math>\epsilon>0</math> i funkcja <math>f:[a,b]\to R</math>, która jest nieskończenie wiele razy różniczkowalna i wszystkie jej pochodne są na <math>[a,b]</math> ograniczone przez <math>M</math>. Napisz program znajdujący stopień <math>n</math> oraz współczynniki w bazie Newtona wielomianu <math>w_{f,n}\in\Pi_n</math> interpolującego <math>f</math> z błędem
-<center><math>\displaystyle \|f-w_{f,n}\|_{C([a,b])}\,\le\,\epsilon.
+<center><math>\|f-w_{f,n}\|_{C([a,b])}\,\le\,\epsilon</math></center>
-</math></center>
 Rozpatrz dwa przypadki: gdy węzły interpolacji są równoodległe, oraz gdy węzły są czebyszewowskie.
 </div></div>

MN09LAB: Różnice pomiędzy wersjami

Aktualna wersja na dzień 21:49, 11 wrz 2023

Interpolacja wielomianowa

Menu nawigacyjne

Działania na stronie

Opcje strony

Narzędzia osobiste

Nawigacja

Szukaj

Narzędzia