Wskaznik uwarunkowania (1)
Ax = b układ o macierzy kwadratowej nieosobliwej
Załóżmy, że A-1 jest zaburzona, co zmienia ją na macierz B
Zaburzenie przenosi się na rozwiązanie x = A-1b
Zamiast wektora x otrzymujemy wektor x = Bb
Jak duże jest bezwzględne i względne zaburzenie rozwiązania ?
Oszacowanie zaburzenia bezwzględnego:
x - x = x - Bb = x - BAx = (I - BA)x Ł I - BA x
Oszacowanie zaburzenia względnego: normy zgodne
x - x
Ł I - BA
x
Wskaznik uwarunkowania (2)
Załóżmy, że zamiast b mamy wektor zaburzony
b
x jest rozwiązaniem równania Ax = b
A jest nieosobliwa
jest rozwiązaniem równania Ax = b
x
Jak różnią się x i bezwzględnie i względnie ?
x
normy zgodne
Oszacowanie zaburzenia bezwzględnego:
x - x = A-1b - A-1b = A-1(b - b) Ł A-1 b - b
Oszacowanie zaburzenia względnego:
b - b b - b
x - x Ł A-1 b - b = A-1 Ax Ł A-1 A x
b b
x - x b - b
Ł k(A)
Stąd gdzie k(A) = A-1 A
x b
k (A) wskaznik uwarunkowania macierzy A (zależy od wybranej normy macierzy)
Wskaznik uwarunkowania. Przykład 1
Przykład 1
1 -1- e
1 1+ e ł
ł
-2
gdzie e > 0
A-1 = e
A =
ę ś
ę1- e 1 ś
1
-1+ e
-2
A-1 Ą = e (2 + e )
A = 2 + e
Ą
norma wierszowa
2
2 + e 4
ć
n
kĄ (A) = A-1 Ą A = >
2
Ą
A = max aij
e
e
Ł ł Ą
1ŁiŁn
j=1
Jeżeli e Ł 0.01 to k (A) > 40 000
Małe zaburzenie względne wektora b może spowodować 40 000 razy większe
zaburzenie względne rozwiązania układu Ax = b
Macierz A o dużym wskazniku k (A) nazywa się x - x b - b
Ł k(A)
zle uwarunkowaną. Jeżeli wskaznik jest niezbyt
x b
duży to macierz jest dobrze uwarunkowana.
Wskaznik uwarunkowania. Przykład 2
Przykład 2
1 3
Dokładne rozwiązanie:
x1 + x2 =
2 2
x1 =1
1 1 5
x2 =1
x1 + x2 =
2 3 6
Zaburzamy nieznacznie wektor b:
1 3
Dokładne rozwiązanie:
x1 + x2 =
2 2
x1 = 0
1 1
x2 = 3
x1 + x2 =1
2 3
Układ jest bardzo czuły na małe zaburzenia wektora b.
norma spektralna
1
1 ł
2
k2 (A) = 19.281
A =
ę1 1 ś
A = max
2
2 3
Funkcja cond
c = cond (A) zwraca wskaznik uwarunkowania k2 macierzy A
c = cond (A, p) zwraca wskaznik uwarunkowania kp macierzy A
p = 1 dla normy kolumnowej
p = 2 dla normy spektralnej
p = inf dla normy wierszowej
Przykład:
1
1 ł
2
A =
ę1 1 ś
2 3
A=[1 0.5; 0.5 1/3];
c = cond (A)
c =
19.2815
Wskaznik uwarunkowania (5)
Rozwiązując numerycznie układ Ax = b, zamiast dokładnego rozwiązania x
otrzymujemy rozwiązanie przybliżone . W celu sprawdzenia jego
x
dokładności obliczamy wektor residualny
r := b - Ax
Różnicę między rozwiązaniem dokładnym x i rozwiązaniem przybliżonym
e := x - x
nazywany wektorem błędu. Wektory te powiązane są zależnością
Ae = r bo Ae = A(x - x) = Ax - Ax = b - Ax = r
jest dokładnym rozwiązaniem układu A x = b z zaburzoną prawą stroną
x
x
b := b - r Jakie są związki pomiędzy błędami względnymi wektorów i b ?
Twierdzenie
Wektory residualny i błędu oraz wskaznik uwarunkowania spełniają nierówność
r e r
1
x - x b - b
e r
Ł Ł k(A)
= =
k(A) b x b
x x b b
Rozwiązywanie układów liniowych metodami iteracyjnymi.
Tworzą ciąg wektorów zbieżny
Metody rozwiązywania
do rozwiązania. Obliczenia
układów Ax = b
przerywamy po osiągnięciu
wymaganej dokładności lub
ustalonej liczbie iteracji.
Metody bezpośrednie
Metody iteracyjne
Gaussa i jej warianty
Richardsona Gaussa-Seidela
Po skończonej liczbie
kroków dałaby dokładne
rozwiązanie gdyby nie
Jacobiego inne
błędy zaokrągleń
Metoda Jacobiego (1).
2x1 - x2 + x3 = -1
x1 = 0.5x2 - 0.5x3 - 0.5
x1 + 2x2 - x3 = 6
x2 = -0.5x1 - 0.5x3 + 3
x1 - x2 + 2x3 = -3
x3 = -0.5x1 + 0.5x2 -1.5
( (
ł ł
2 -1 1 x1 -1
x1k ) 0 0.5 - 0.5 x1k -1) - 0.5
ł ł ł
ł ł
ę ś ę ś
ę ś ę ś
(k ) (k -1)
ę1 2 -1ś ęx ś ę ś
= 0.5 0 0.5 + 3
= 6
2 2
ęx ś ęx ś
2
ę- ś ę ś
ę ś ę ś ę ś
( ( -1)
ęx3k ) ś
ę
ę
- 0.5 0.5 0 ś ęx3k ś ę
-1.5ś
3
1 -1 2 ś ę ś ę 3
x - ś
Wartości początkowe:
0
ł
ę0ś
x(0) =
ę ś
ę
0ś
Metoda Jacobiego (2).
0
Krok 1: ł
ę0ś
x(0) =
( ę ś
ł
x11) 0 0.5 - 0.5 0 - 0.5 - 0.5
ł ł ł ł
ę
ę ś
0ś
ę ś ę0ś ę ś ę ś
(1)
= 0.5 0 0.5 + 3 = 3
2
ęx ś
ę- ś ę ś ę ś ę ś
(
ęx31) ś
ę
- 0.5 0.5 0 ś ę -1.5ś ę
0ś ę -1.5ś
Krok 2:
(
ł
x12) 0 0.5 - 0.5 - 0.5 - 0.5 1.75
ł ł ł ł
ę ś
ę ś ę ś ę ś ę2.50ś
(2)
= 0.5 0 0.5 3 + 3 =
2
ęx ś
ę- ś ę ś ę ś ę ś
(
ęx32) ś
ę
- 0.5 0.5 0 ś ę -1.5ś ę
-1.5ś ę 0.25ś
Dokładne rozwiązanie:
Po 13 iteracjach:
1.0002 1
ł Kryterium zakończenia ł
ę ś ę ś
obliczeń:
x(13) = 2.0001 x(13) = 2
ę ś ę ś
ę ś || x(k) - x(k-1) ||2 < 0.001 ę
- 0.9997 -1ś
Metoda Jacobiego. Algorytm (1)
input n, (aij), (bi), (xi), M,
for k = 1 to M do
for i = 1 to n do
ui Ź bi
for j = 1 to n do
if j ą i then
n
ui Ź ui - aijxj
ui Ź (bi -
ij j
a x ) / aii
end if
j=1, jąi
end do
ui Ź ui / aii
end do
for i = 1 to n do
xi Ź ui
output k, (xi )
end do
end do
Algorytm można ulepszyć
Algorytm jest wykonywany
wykonując wszystkie dzielenia
dla zadanej liczby kroków M.
na początku programu.
input n, (aij), (bi), (xi), M,
Metoda Jacobiego. Algorytm (2)
for i = 1 to n do
d Ź 1/ aii
bi Ź d bi Wszystkie dzielenia są
for j = 1 to n do
wykonywane na początku
aij Ź d aij
programu
end do
end do
for k = 1 to M do
Dzięki temu zasadnicza
for i = 1 to n do
instrukcja podstawienia
ui Ź bi
upraszcza się do postaci
for j = 1 to n do
n
if j ą i then
ui Ź bi -
ij
a xj
ui Ź ui - aijxj
j=1, jąi
end if
end do
end do
for i = 1 to n do
xi Ź ui
output k, (xi )
end do
end do
Metoda Jacobiego (3).
Twierdzenie
Jeżeli macierz A jest dominująca przekątniowo, to dla dowolnego wektora
początkowego metoda Jacobiego tworzy ciąg zbieżny do rozwiązania
układu Ax = b.
Macierze określające metodę:
Macierz dominująca przekątniowo:
Ax = b n
aii > aij dla 1Ł i Ł n
(D - CL - CU) x = b
j=1, jąi
Przykład:
Dx = (CL + CU) x + b
4 -1 0 -1
ł
x = D-1(CL + CU) x + D-1b
ę
ę-1 4 0 -1ś
ś
x(k) = D-1(CL + CU) x(k -1) + D-1b
ę ś
-1 0 4 -1
ę ś
0 -1 -1 4
D - macierz diagonalna
CL - macierz trójkątna dolna bez głównej przekątnej
CU - macierz trójkątna górna bez głównej przekątnej
Metoda Gaussa-Seidela (1).
2x1 - x2 + x3 = -1
x1 = 0.5x2 - 0.5x3 - 0.5
x1 + 2x2 - x3 = 6
x2 = -0.5x1 - 0.5x3 + 3
x1 - x2 + 2x3 = -3
x3 = -0.5x1 + 0.5x2 -1.5
( ( (
x1k ) = 0.5x2k -1) - 0.5x3k -1) - 0.5
( ( (
x2k ) = -0.5x1k ) - 0.5x3k -1) + 3
( ( (
x3k ) = -0.5x1k ) + 0.5x2k ) -1.5
Wartości początkowe:
0
ł
ę0ś
x(0) =
ę ś
ę
0ś
Metoda Gaussa-Seidela (2).
0
ł
Krok 1:
ę0ś
x(0) =
ę ś
(
x11) = 0.5 0 - 0.5 0 - 0.5 = -0.5
ę
0ś
(
x21) = -0.5 (-0.5) - 0.5 0 + 3 = 3.25
(
x31) = -0.5 (-0.5) + 0.5 (3.25) -1.5 = 0.375
Krok 2:
(
x12) = 0.53.25- 0.5 0.375- 0.5 = 0.9375
(
x22) = -0.5 0.9375- 0.5 0.375+ 3 = 2.7188
(
x32) = -0.5 0.9375+ 0.5 2.7188-1.5 = -0.6094
Dokładne rozwiązanie:
Po 10 iteracjach:
1.0001 1
ł Kryterium zakończenia ł
ę ś ę ś
obliczeń:
x(10) = 1.9999 x(13) = 2
ę ś ę ś
ę || x(k) - x(k-1) ||2 < 0.001 ę
-1.0001ś -1ś
Metoda Gaussa-Seidela. Algorytm
input n, (aij), (bi), (xi), M,
for i = 1 to n do
d Ź 1/ aii
Jeżeli wszystkie dzielenia wykonujemy na
bi Ź d bi
for j = 1 to n do
początku programu wówczas zasadnicza
aij Ź d aij
instrukcja podstawienia
end do
n
end do
xi Ź (bi -
ij
a xj ) / aii
for k = 1 to M do
j=1, jąi
for i = 1 to n do
xi Ź bi
upraszcza się do postaci
for j = 1 to n do
n
if j ą i then
xi Ź bi -
ij j
a x
xi Ź xi - aijxj
j=1, jąi
end if
end do
end do
Algorytm jest wykonywany
output k, (xi )
dla zadanej liczby kroków M.
end do
Metoda Gaussa-Seidela (3).
Twierdzenie
Jeżeli macierz A jest dominująca przekątniowo, to dla dowolnego wektora
początkowego metoda Gaussa-Seidela tworzy ciąg zbieżny do rozwiązania
układu Ax = b.
Macierze określające metodę:
Macierz dominująca przekątniowo:
n
Ax = b
aii > aij dla 1Ł i Ł n
(D - CL - CU) x = b
j=1, jąi
Przykład:
(D - CL) x = CU x + b
4 -1 0 -1
ł
x = (D - CL)-1CU x + (D - CL)-1b
ę
ę-1 4 0 -1ś
ś
x(k) = (D - CL)-1CU x(k -1) + (D - CL)-1b
ę ś
-1 0 4 -1
ę ś
0 -1 -1 4
D - macierz diagonalna
CL - macierz trójkątna dolna bez głównej przekątnej
CU - macierz trójkątna górna bez głównej przekątnej
Różnica pomiędzy metodami Gaussa-Seidela i Jacobiego.
Gauss-Seidel Jacobi
Krok k
Krok k
x1 = (b1 - a12x2 - a13x3)/ a11 x1 = (b1 - a12x2 - a13x3)/ a11
x2 = (b2 - a21x1 - a23x3)/ a22
x2 = (b2 - a21x1 - a23x3)/ a22
x3 = (b3 - a31x1 - a23x2)/ a33
x3 = (b3 - a31x1 - a23x2)/ a33
Krok k+1
Krok k+1
x1 = (b1 - a12x2 - a13x3)/ a11 x1 = (b1 - a12x2 - a13x3)/ a11
x2 = (b2 - a21x1 - a23x3)/ a22
x2 = (b2 - a21x1 - a23x3)/ a22
x3 = (b3 - a31x1 - a23x2)/ a33
x3 = (b3 - a31x1 - a23x2)/ a33
Funkcje macierzowe
c = norm (A) zwraca normę spektralną macierzy albo euklidesową wektora A
c = norm (A, p) zwraca normę macierzy albo wektora A
p = 1 dla normy kolumnowej macierzy albo sumy wartości
bezwzględnych współrzędnych wektora
p = 2 dla normy spektralnej macierzy albo euklidesowej wektora
p = inf dla normy wierszowej macierzy albo maksimum wektora
p = -inf dla minimalnej wartości bezwzględnej współrzędnych
wektora
k = rank (A) szacuje ilość liniowo niezależnych wierszy lub kolumn macierzy A
d = det (A) zwraca wyznacznik macierzy kwadratowej A
b = trace (A) zwraca sumę elementów na przekątnej macierzy A
Funkcje algebry liniowej
\ - operator lewostronnego dzielenia macierzy
Jeżeli A i B są macierzami kwadratowymi to A\B jest równe A-1B.
Jeżeli A jest macierzą kwadratową, zaś b wektorem to A\b jest
rozwiązaniem układu równań Ax = b.
/ - operator prawostronnego dzielenia macierzy
Jeżeli A i B są macierzami kwadratowymi to B/A jest równe BA-1.
Jeżeli A jest macierzą kwadratową, zaś b wektorem to b/A jest
rozwiązaniem układu równań xA = b.
[L, U] = lu (A) zwraca macierz trójkątną górną U oraz permutację macierzy
trójkątnej dolnej (tj. iloczyn macierzy trójkątnej dolnej i macierzy
permutacji) tak, że A = LU
[L, U, P] = lu (A) zwraca macierz trójkątną górną U, macierz trójkątną dolną
z jedynkami na głównej przekątnej oraz macierz permutacji P
tak, że LU = PA
Y = inv (A) zwraca macierz odwrotną macierzy kwadratowej A tj. Y = A-1.
Dopasowanie krzywych
Dopasowanie krzywych
Interpolacja Aproksymacja
x0, x1, ..., xn
y0, y1, ..., yn
Interpolacja
x0, x1, ..., xn
Interpolacja
y0, y1, ..., yn
Wielomianowa Ilorazy różnicowe Hermita
Funkcje wymierne
Funkcje sklejane
Twierdzenie
Lagrange a Newtona
Jeżeli liczby x0, x1, ..., xn są parami różne,
to istnieje dokładnie jeden wielomian pn
taki, że pn(xi) = yi dla 0 Ł i Ł n
Interpolacja wielomianami Lagrange a
Dla dwóch punktów (x0, y0) i (x1, y1) wielomian Lagrange a jest linią prostą i
przechodzi przez te punkty:
x - x1 x - x0
p1(x) = y0 + y1
x0 - x1 x1 - x0
Dla trzech punktów (x0, y0) i (x1, y1) i (x2, y2) wielomian Lagrange a jest
parabolą i przechodzi przez te punkty:
(x - x1)(x - x2) (x - x0)(x - x2) (x - x0)(x - x1)
p2(x) = y0 + y1 + y2
(x0 - x1)(x0 - x2) (x1 - x0)(x1 - x2) (x2 - x0)(x2 - x1)
Ogólna postać wielomianów Lagrange a
Dla n+1 danych
pn (x) = L0 y0 + L1y1 + ...+ Lk yk + ...+ Ln yn
punktów dostajemy
gdzie
wielomian stopnia n.
(x - x0)...(x - xk -1)(x - xk +1)...(x - xn )
Lk (x) =
(xk - x0)...(xk - xk -1)(xk - xk +1)...(xk - xn )
Lk (xk ) =1
Wielomiany Lk zależą od węzłów x0, x1, ..., xn,
Lk (xj ) = 0 dla j ą k
ale nie zależą od wartości y0, y1, ..., yn.
Przykład 1. Interpolacja Parabolą.
Dane punkty:
Wstawiając do ogólnego wzoru dostajemy:
(x0, y0) = (-2, 4),
(x - 0)(x - 2) (x - (-2))(x - 2) (x - (-2))(x - 0)
p2(x) = 4 + 2 + 8
(x1, y1) = (0, 2),
(-2 - 0)(-2 - 2) (0 - (-2))(0 - 2) (2 - (-2))(2 - 0)
(x2, y2) = (2, 8)
8
Po uproszczeniu:
7
p2(x) = x2 + x + 2
6
5
4
Uwaga:
3
Dla innych punktów
dostaniemy inny
2
wielomian Lagrange a.
1
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
Interpolacja wielomianami Lagrange a. Algorytm 1.
x wektor węzłów function c = Lagrange_Coef (x, y)
y wektor wartości n = length(x);
c wektor współczynników for k = 1 : n
d(k) = 1;
for i = 1 : n
if i ~= k
d(k) = d(k)*(x(k)-x(i));
end
Obliczanie współczynników ck
c(k) = y(k)/d(k);
wielomianów Lagrange a
end
Uwaga:
przedstawionych następująco:
end
Numeracje ck
pn-1(x) = c0N0 + c1N1 + ...+ cn-1Nn-1
we wzorze
yk
oraz funkcji
gdzie ck =
są różne.
(xk - x0)...(xk - xk -1)(xk - xk +1)...(xk - xn-1)
Lk yk = ck Nk
Nk (x) = (x - x0)...(x - xk-1)(x - xk+11)...(x - xn-1)
Interpolacja wielomianami Lagrange a. Algorytm 2.
Obliczanie wartości function p = Lagrange_Eval (t, x, c)
wielomianów Lagrange a m = length(x);
w punkcie t. for i = 1 : length(t)
p(i) = 0;
for j = 1 : m
t punkt albo wektor
N(j) = 1;
zmiennej niezależnej
for k = 1 : m
x wektor węzłów
if j ~= k
c wektor współczynników
N(j) = N(j)*(t(i)-x(k));
p wektor wartości
end
w punktach t
end
p(i) = p(i)+N(j)*c(j);
end
end
Przykład 2. Stężenie produktu reakcji chemicznej
0.35
Czas Stężenie
0.0 0.0
0.3
0.5 0.19
0.25
1.0 0.26
1.5 0.29
0.2
2.0 0.31
0.15
0.1
0.05
0
0 0.5 1 1.5 2
współczynniki ck: Czas reakcji
0 -0.5067 1.0400 -0.7733 0.2067
Udzial molowy produktu
Przykład 2. Dodatkowe dane
Czas Stężenie
0.4
0.0 0.0
0.1 0.06
0.35
0.4 0.17
0.3
0.5 0.19
0.25
0.6 0.21
0.2
0.9 0.25
0.15
1.0 0.26
0.1
1.1 0.27
1.4 0.29
0.05
1.5 0.29
0
1.6 0.30
-0.05
0 0.5 1 1.5 2
1.9 0.31
Czas reakcji
2.0 0.31
Udzial molowy produktu
Zalety i wady wielomianów Lagrange a
Zalety:
" Niezależność wielomianów Lk od rzędnych yi , co przydaje się w rozważaniach
analitycznych
" Wygodne w sytuacji gdy mamy kilka serii danych pomiarowych różniących
się tylko wartościami zmiennej zależnej, natomiast wartości zmiennej
niezależnej (węzłów) są dla wszystkich serii takie same.
Wady:
" Niewygodne gdy musimy dołączyć dodatkowe punkty danych lub gdy stopień
wielomianu interpolacyjnego jest nieznany (tzn. gdy być może lepiej użyć
niepełnego zbioru danych pomiarowych).
Interpolacja wielomianami Newtona
Dla dwóch punktów (x0, y0) i (x1, y1) wielomian Newtona jest linią prostą
przechodzącą przez te punkty:
p1(x) = a0 + a1(x - x0)
Dla trzech punktów (x1, y1), (x2, y2) i (x3, y3) wielomian Newtona jest
parabolą przechodzącą przez te punkty:
p2(x) = a0 + a1(x - x0)+ a2(x - x0)(x - x1)
Dla n punktów:
pn(x) = a0 + a1(x - x0)+ a2(x - x0)(x - x1)+...+ an(x - x0)(x - x1)...(x - xn-1)
Obliczanie współczynników wielomianu Newtona
pn(x) = a0 + a1(x - x0)+ a2(x - x0)(x - x1)+...+ an(x - x0)(x - x1)...(x - xn-1)
(x0, y0)
a0 = pn(x0) = y0
y1 - y0
y2 - y1 y1 - y0
a0 i (x1, y1)
a1 =
-
x1 - x0
x2 - x1 x1 - x0
a2 =
x2 - x0
a0, a1, i (x2, y2)
itd.
Przykład 1. Interpolacja Newtona Parabolą.
Dane punkty:
(x1, y1) = (-2, 4), Wstawiając a0, a1, a2, x1 i x2 do wzoru
(x2, y2) = (0, 2),
p2(x) = a0 + a1 x - x0 + a2 x - x0 x - x1
( ) ( ) ( )
(x3, y3) = (2, 8)
po uproszczeniu dostajemy:
a0 = y0 = 4
p2(x) = x2 + x + 2
y1 - y0 2 - 4
a1 = = = -1
x1 - x0 0 - (-2)
y2 - y1 y1 - y0
8 - 2 2 - 4
-
-
x2 - x1 x1 - x0
2 - 0 0 - (-2)
a2 = = =1
x2 - x0 2 - (-2)
Ilorazy różnicowe (1)
Jeżeli liczby x0, x1, ..., xn są parami różne, to wielomian pn taki, że p (xi) = yi
dla 0 Ł i Ł n można wyrazić na pomocą wzoru Newtona:
p(x) = a0 + a1(x - x0)+ a2(x - x0)(x - x1)+...+ an(x - x0)(x - x1)...(x - xn-1)
Współczynniki ak zależą tylko od węzłów x0, x1, ..., xk i od wartości funkcji w
tych węzłach. Zależność tę przyjęto oznaczać następująco:
Iloraz różnicowy rzędu 0 dla
a0 = f [x0 ]
funkcji f i węzła x0
a1 = f [x0, x1]
a2 = f [x0, x1, x2]
M
Iloraz różnicowy rzędu k dla
ak = f [x0, x1,K, xk -1, xk ]
funkcji f i węzłów x0, x1, ... , xk
M
an = f [x0, x1,K, xn-1, xn ]
Ilorazy różnicowe (2)
a0 = y0
f [x0] = f (x0)
y1 - y0
f (x1) - f (x0)
a1 =
f [x0, x1] =
x1 - x0
x1 - x0
y2 - y1 y1 - y0
-
f [x1, x2]- f [x0, x1]
x2 - x1 x1 - x0
a2 =
f [x0, x1, x2] =
x2 - x0
x2 - x0
Twierdzenie
Ilorazy różnicowe spełniają zależność
f [x1, x2,...,xk ] - f [x0, x1,...,xk-1]
f [x0, x1,...,xk ] =
xk - x0
Ilorazy różnicowe (3)
Oznaczenia węzłów nie są istotne. Ogólne zależności są następujące:
f [xi+1, xi+2,...,xi+ j ] - f [xi , xi+1,...,xi+ j-1]
f [xi , xi+1,...,xi+ j ] =
xi+ j - xi
Znając węzły xi i wartości funkcji f (xi), czyli ilorazy f [xi] zerowego
rzędu można za pomocą powyższego wzoru utworzyć tablicę ilorazów
różnicowych wyższych rzędów.
Tablica dla 4 węzłów:
x0 f [x0] f [x0, x1] f [x0, x1, x2] f [x0, x1, x2, x3]
x1 f [x1] f [x1, x2] f [x1, x2, x3]
x2 f [x2] f [x2, x3]
f (x1) - f (x0)
f [x0, x1] =
x1 - x0
x3 f [x3]
f [x1, x2]- f [x0, x1]
dane f [x0, x1, x2] =
x2 - x0
Ilorazy różnicowe. Przykład (1)
Utworzyć tablicę ilorazów różnicowych i podać postać wzoru interpolacyjnego
Newtona.
xi f (xi ) f (x1) - f (x0) - 3 -1
f [x0, x1] = = = 2
3 7
x1 - x0 1- 3
3 1 2 -
8 40
5 3
1 - 3
4 20
f (x2) - f (x1) 2 - (-3) 5
f [x1, x2] = = =
5 2 2
x2 - x1 5 -1 4
6 4
f (x3) - f (x2) 4 - 2
f [x2, x3] = = = 2
x3 - x2 6 - 5
dane
f [xi ] = f (xi )
Ilorazy różnicowe. Przykład (2)
Utworzyć tablicę ilorazów różnicowych i podać postać wzoru interpolacyjnego
Newtona.
5
- 2
xi f (xi )
f [x1, x2]- f [x0, x1] 3
4
f [x0, x1, x2] = = = -
3 7
3 1 2 -
8 40
x2 - x0 5 - 3 8
5 3
1 - 3
4 20
5
2 -
5 2 2
f [x2, x3]- f [x1, x2] 3
4
f [x1, x2, x3] = = =
6 4
x3 - x1 6 -1 20
dane
33
ć
- -
7
f [x1, x2, x3]- f [x0, x1, x2]
20 8
Ł ł
f [x0, x1, x2, x3] = = =
x3 - x0 6 - 3 40
37
p(x) =1+ 2 x - 3 - x -3 x -1 + x - 3 x -1 x -5
( ) ( ) ( ) ( ) ( ) ( )
8 40
Ilorazy różnicowe. Przykład (3)
p(x) = 1 + 2(x-3) - (3/8)(x-3)(x-1) + (7/40)(x-3)(x-1)(x-5)
10
8
6
xi f (xi )
4
3 1
2
1 - 3
0
5 2
-2
6 4
-4
-6
-8
-10
0 1 2 3 4 5 6 7
x
p(x)
Ilorazy różnicowe. Algorytm
cij = f [xi, xi+1,...,xi+ j ]
Oznaczmy
x0 c00 c01 c02 K c0,n-1 c0n
Wystarczy użyć tablicy d
z jednym wskaznikiem.
x1 c10 c11 c12 K c1,n-1
Uwaga:
M M M
wskaznik tablicy d zaczyna
xn-1 cn-1,0 cn-1,1 się od 0.
xn cn0
for i = 0 to n do
di Ź f (xi)
end do
ci0 = f [xi ] = f (xi )
for j = 1 to n do
for i = n to j step - 1 do
Początkowa wartość zmiennej di = ci0 = f (xi).
di Ź (di - di-1) / (xi - xi-j)
Następnymi wartościami są cn-1,1, ... , c11, c01, c00.
end do
Tablicę tworzymy kolumnami, a każdą kolumnę
end do
od dołu do góry.
Zalety wielomianów Newtona
Zalety:
" Najbardziej użyteczna w obliczeniach numerycznych jest wersja wzoru
interpolacyjnego Newtona zawierająca ilorazy różnicowe.
" Wersja ta pozwala na dołączanie dodatkowych punktów (xi, yi) bez potrzeby
przeliczania wcześniej wyznaczonych współczynników ci.
" Wartość wielomianu łatwo policzyć stosując schemat Hornera
" Oszacowanie błędu może być w prosty sposób włączone do algorytmu
wn(x) = anxn + an-1xn-1 + ... + a1x + a0
Schemat Hornera:
w0 = an
wk = wk -1x + an-k dla k = 1, 2, ..., n
Błędy interpolacji wielomianowej
Twierdzenie
m+1
Jeśli f C [a, b], a wielomian p Pn interpoluje wartości funkcji f w n+1
różnych punktach x0, x1, ... , xn przedziału [a, b], to dla każdego x [a, b]
istnieje takie xx (a, b), że
(n+1)
f (xx )
f (x) - p(x) = (x - x0)(x - x1)K(x - xn )
(n +1)!
Przykład
Jaki jest błąd przybliżenia funkcji sin x wielomianem interpolacyjnym 9 stopnia
w przedziale [0, 1] do którego należą węzły?
sin(10) xx Ł1 oraz (x - x0)(x - x1)K(x - x9) Ł1
więc
m
C przestrzeń funkcji zmiennej
1
sin x - p(x) Ł < 2.810-7
rzeczywistej mających m-tą pochodną
10!
ciągłą
Pn zbiór wszystkich wielomianów
stopnia co najwyżej n-tego
Funkcje sklejane. Wady interpolacji wielomianowej.
0.35
0.3
Wielomian 4-go stopnia
0.25
0.2
0.15
0.4
0.1
0.05
0.3
0
0 0.5 1 1.5 2
0.2
Czas reakcji
0.1
Wielomian 12-go stopnia
0
-0.1
0 0.5 1 1.5 2
Czas reakcji
Udzial molowy produktu
Udzial molowy produktu
Funkcje sklejane. Definicja.
Mamy n+1 węzłów x0, x1, ... , xn takich, że x0< x1< ... < xn. Niech k będzie liczbą
całkowitą i nieujemną. Funkcją sklejaną stopnia k nazywamy taką funkcję S,
która:
1. W każdym z przedziałów [xi, xi+1) dla 0 Ł i Ł n-1 jest wielomianem
klasy Pk.
2. Ma ciągłą (k-1) pochodną w przedziale [x0, xn].
Liniowe fukcje sklejane
3
xi yi
3.0 2.5
2
4.5 1.0
7.0 2.5
1
9.0 0.5
0
2 3 4 5 6 7 8 9 10
Kwadratowe funkcje sklejane
Si (x) = aix2 + bix + ci
Dla n+1 węzłów (i = 0, 1, ... , n) mamy n
przedziałów 3n niewiadomych (ai, bi, ci)
" wartości funkcji przyległych
wielomianów mają być równe w węzłach
Si-1(xi ) = yi
dla 1Ł i Ł n -1
Si (xi ) = yi
2n-2 warunki
" pierwsza i ostatnia funkcja muszą przechodzić
2 warunki
przez punkty końcowe
n-1 warunków
S1(x0) = y0, Sn(xn) = yn
Dodatkowy warunek:
" pierwsze pochodne w węzłach wewnętrznych
mają być równe
óó
S1(x0) = 0
ó ó
Si-1(xi ) = Si (xi ) dla 1Ł i Ł n -1
Kwadratowe funkcje sklejane. Przykład (1)
Si (x) = aix2 + bix + ci
20.25a1 + 4.5b1 + c1 =1.0
xi yi
Si-1(xi ) = yi
20.25a2 + 4.5b2 + c2 =1.0
3.0 2.5
49a2 + 7b2 + c2 = 2.5 Si (xi ) = yi
4.5 1.0
49a3 + 7b3 + c3 = 2.5
7.0 2.5
S1(x0) = y0
9a1 + 3b1 + c1 = 2.5
9.0 0.5
Sn (xn ) = yn
81a3 + 9b3 + c3 = 0.5
9a1 + b1 = 9a2 + b2
ó ó
Si-1(xi ) = Si (xi )
ó
Si (x) = 2aix + bi
14a2 + b2 =14a3 + b3
óó óó
Si (x) = 2ai S1(x0) = 0
a1 = 0
Mamy układ 9 równań z 9 niewiadomymi. Ponieważ a1 = 0 więc zostaje układ
8 równań z 8 niewiadomymi.
Kwadratowe funkcje sklejane. Przykład (2)
4.5 1 0 0 0 0 0 0 b1 1
ł ł ł
ę ś
0 0 20.25 4.5 1 0 0 0ś ęc1 ś ę 1
ę ś ę ś ę ś
ę ś ę ś ę ś
0 0 49 7 1 0 0 0 a2 2.5
ę ś ę ś ę ś
0 0 0 0 0 49 7 1ś ęb2 ś ę2.5ś
ę
=
ę
3 1 0 0 0 0 0 0ś ęc2 ś ę2.5ś
ę ś ę ś ę ś
0 0 0 0 0 81 9 1ś ęa3 ś ę0.5ś
ę
ę ś
1 0 - 9 -1 0 0 0 0ś ęb3 ś ę 0
ę ś ę ś ę ś
0 0 14 1 0 -14 -1 0 c3 0
ę ś ę ś ę ś
a1 = 0 b1 = -1 c1 = 5.5
a2 = 0.64 b2 = -6.76 c2 =18.46
a3 = -1.6 b3 = 24.6 c3 = -91.3
S1(x) = -x + 5.5 3.0 Ł x Ł 4.5
S2(x) = 0.64x2 - 6.76x +18.46 4.5 Ł x Ł 7.0
S3(x) = -1.6x2 + 24.6x - 91.3 7.0 Ł x Ł 9.0
Kwadratowe funkcje sklejane. Przykład (3)
Kwadratowe funkcje sklejane
4
3
2
1
0
2 3 4 5 6 7 8 9 10
Liniowe fukcje sklejane
3
2
1
0
2 3 4 5 6 7 8 9 10
Sześcienne funkcje sklejane (1)
Si (x) = aix3 + bix2 + cix + di
Dla n+1 węzłów (i = 0, 1, ... , n) n przedziałów 4n niewiadomych (ai, bi, ci, di)
" wartości funkcji przyległych wielomianów mają być równe w węzłach
Si-1(xi ) = yi
dla 1Ł i Ł n -1
2n-2 warunki
Si (xi ) = yi
" pierwsza i ostatnia funkcja muszą przechodzić przez punkty końcowe
S1(x0) = y0, Sn(xn) = yn
2 warunki
" pierwsze pochodne w węzłach wewnętrznych mają być równe
ó ó
Si-1(xi ) = Si (xi ) dla 1Ł i Ł n -1
n-1 warunków
" drugie pochodne w węzłach wewnętrznych mają być równe
óó óó
Si-1(xi ) = Si (xi ) dla 1Ł i Ł n -1
n-1 warunków
" drugie pochodne w węzłach końcowych są równe zero
2 warunki
óó óó
S1(x0) = 0, Sn(xn) = 0
Taką funkcję nazywamy naturalną funkcją sklejaną.
Razem: 4n warunków
Sześcienne funkcje sklejane (2)
Można pokazać, że
ć ć
zi zi+1 yi+1 zi+1hi yi zihi
Si (x) = (xi+1 - x)3 + (x - xi )3 + -
(x - xi ) + hi - 6
(xi+1 - x)
6hi 6hi hi 6
Ł ł Ł ł
óó
gdzie hi = xi+1 - xi, zi = S (xi )
zi możemy wyznaczyć z układu równań
6 6
hi-1zi-1 + 2(hi-1 + hi ) + hizi+1 = (yi+1 + yi ) - (yi + yi-1) dla 1Ł i Ł n -1
hi hi-1
Mamy n-1 równań i n+1 niewiadomych (zi). Można przyjąć z0 = zn = 0. Oznaczamy
6
hi = xi+1 - xi, ui = 2(hi-1 + hi ), bi = (yi+1 - yi ), vi = bi -bi-1
hi
u1 h1 z1 v1
ł ł ł
ęh u2 h2 ś ę
z2 ś ę v2 ś
1
ę ś ę ś ę ś
ę ś ę ś ę ś
=
ę Macierz jest symetryczna,
hn-3 un-2 hn-2 ś ęzn-2 ś ęvn-2 ś
ę ś ę ś ę ś
trójprzekątniowa
ę
i przekątniowo dominująca
hn-2 un-1 ś ę zn-1 ś ęvn-1 ś
Sześcienne funkcje sklejane (3)
Twierdzenie
Jeżeli f C2[a, b], a = x0< x1< ... < xn = b, a S jest naturalną funkcją sklejaną,
interpolującą f w węzłach xi dla 0 Ł i Ł n, to
b b
2 2
[Sóó(x)] dx Ł [f óó(x)] dx
a a
Znaczy to, że w pewnym sensie naturalna funkcja sklejana sześcienna jest
najgładszą funkcją interpolującą.
Wyszukiwarka
Podobne podstrony:
Wyklad studport 9Wyklad studport 5Wyklad studport 11Wyklad studport 7Wyklad studport 12Wyklad studport 14Wyklad studport 4Wyklad 12 Podstawowe typy zwiazkow chemicznych blok s i p PCHN SKP studportWyklad 5 Uklad okresowy pierwiastkow studportWyklad 10 Elektrolity, woda, kwasy i zasady PCHN SKP studportWYKŁAD 4 Elektronowa struktura atomu studportSieci komputerowe wyklady dr FurtakWykład 05 Opadanie i fluidyzacjaWYKŁAD 1 Wprowadzenie do biotechnologii farmaceutycznejwięcej podobnych podstron