WYKA
AD 13
APROKSYMACJA ÅšREDNIOKWADRATOWA WIELOMIANAMI
RADZAJE APROKSYMACJI
Aproksymacja jest działem analizy numerycznej zajmującym się najbardziej ogólnymi zagadnieniami
f (x)
przybliżania funkcji, polegającymi na wyznaczaniu dla danej funkcji takich funkcji F(x),
f (x).
które w określonym sensie najlepiej przybliżają funkcję
Zadania aproksymacyjne mogą być formułowane bardzo różnie, w zależności od przyjętego sposobu
oszacowania błędów aproksymacji. Wyróżnia się trzy rodzaje aproksymacji:
1) interpolacyjnÄ…,
2) jednostajnÄ…,
3) średniokwadratową.
Rys. 1
W przypadku aproksymacji interpolacyjnej, podobnie jak w zagadnieniu interpolacji, żądamy
f (x) F(x)
spełnienia warunku, aby funkcja dana i funkcja szukana przyjmowały dokładnie te same
wartości na zbiorze z góry ustalonych punktów węzłowych (rys. 1). Warunek ten może być uzupełniony
warunkami wyrażającymi równość pochodnych w węzłach, jeżeli wartości pochodnych zostaną zadane.
Rys. 2
f (x) F(x),
W przypadku aproksymacji jednostajnej funkcję przybliżamy taką funkcją która daje
F(x), f (x)
najmniejsze maksimum różnicy między a w całym przedziale [a, b] - rys. 2
max F(x) - f (x) = min.
(1)
x "[a, b]
F(x)
W przypadku aproksymacji średniokwadratowej jako miarę odchylenia funkcji od danej funkcji
f (x) przyjmujemy wielkość
b
2
S = [F(x) - f (x)] d x,
(2)
+"
a
F(x)
zwanÄ… odchyleniem kwadratowym. Funkcja aproksymujÄ…ca wyznaczana jest z warunku,
aby wartość wyrażenia (2) była możliwie najmniejsza. Geometrycznie warunek ten wyraża żądanie,
f (x)
F(x)
aby pole powierzchni między liniami reprezentującymi funkcję oraz funkcję było minimalne.
Zagadnienia aproksymacji jednostajnej i aproksymacji średniokwadratowej są również formułowane dla
funkcji określonych na dyskretnym zbiorze argumentów. Dla takich funkcji warunek (1) dotyczący
aproksymacji jednostajnej zmienia się w ten sposób, że zamiast ciągłej zmiennej niezależnej x występuje
xi
w nim zmienna dyskretna
max F(xi ) - f (xi ) = min. (i = 0,1, ..., n),
(3)
a w warunku (2) na minimum odchylenia kwadratowego całka jest zastępowana sumą
n
2
S =
"[F(x ) - f (xi )] = min .
i
(4)
i = 0
Aproksymacja średniokwadratowa funkcji określonych na dyskretnym zbiorze argumentów jest
najczęściej wykorzystywana w zastosowaniach praktycznych do wygładzania danych eksperymentalnych
i wyników obliczeń ze względu na mniej skomplikowane algorytmy jej realizacji numerycznej
w porównaniu z algorytmami aproksymacji jednostajnej i możliwość uzyskiwania dobrych przybliżeń
f (x).
funkcji
APROKSYMACJA ÅšREDNIOKWADRATOWA
WIELOMIANAMI
F(x)
W zadaniach aproksymacji średniokwadratowej wielomianami funkcji aproksymującej wygodnie
jest poszukiwać - podobnie jak w zadaniu interpolacji - w postaci wielomianu uogólnionego.
m
F(x) = a Õ (x),
" j j
(5)
j =0
y = f (x) określonej na dyskretnym
Rozważymy najpierw aproksymację średniokwadratową funkcji
aj
zbiorze argumentów. W tym przypadku współczynniki ( j = 0, 1, ..., m) funkcji (5) dobieramy
tak, żeby funkcja
2
n m
îÅ‚ Å‚Å‚
S(a0, a1, ..., am ) = a Õ (xi ) - f (xi )śł
" " j j
ïÅ‚
(6)
i = 0 j = 0
ðÅ‚ ûÅ‚
osiągnęła wartość minimalną.
S(a0, a1,..., am)
Zgodnie z ogólnymi metodami rachunku różniczkowego funkcja
osiÄ…ga minimum wtedy
a0, a1, .., am:
i tylko wtedy, gdy znikają pochodne cząstkowe względem wszystkich zmiennych
n m
1 "S îÅ‚ Å‚Å‚
= a Õ (xi ) - f (xi )śł Õ0 (xi ) = 0,
" "
2 "a0 i =0 ïÅ‚ j j
j=0
ðÅ‚ ûÅ‚
n m
1 "S îÅ‚ Å‚Å‚
= a Õ (xi ) - f (xi )śł Õ1(xi ) = 0,
" "
2 "a1 i = 0 ïÅ‚ j j
j=0
ðÅ‚ ûÅ‚
.......................................................................
n m
1 "S îÅ‚ Å‚Å‚
= a Õ (xi ) - f (xi )śł Õm (xi ) = 0.
" "
2 "am i =0 ïÅ‚ =0 j j
j
ðÅ‚ ûÅ‚
aj
m + 1 m + 1
Stąd otrzymujemy układ równań z niewiadomymi współczynnikami , zwany układem
normalnym:
a0 (Õ0, Õ0 ) + a1(Õ1, Õ0 ) + ... + am (Õm, Õ0 ) = ( f , Õ0 ),
üÅ‚
ôÅ‚
a0 (Õ0, Õ1) + a1(Õ1, Õ1) + ... + am (Õm, Õ1) = ( f , Õ1),
ôÅ‚
żł
(7)
................................................................................
ôÅ‚
ôÅ‚
a0 (Õ0, Õm ) + a1(Õ1, Õm ) + ... + am (Õm , Õm ) = ( f , Õm ),
þÅ‚
w którym wprowadzono skrócone oznaczenie
n
(Õ, È) = Õ(xi ) È (xi ).
"
(8)
i = 0
Õ0(x),
Układ równań (7) ma dokładnie jedno rozwiązanie dla liniowo-niezależnego układu funkcji:
Õ1(x), Õm(x).
..., Macierz współczynników układu (7) jest macierzą symetryczną i dodatnio określoną.
j
x
Dla układu funkcji bazowych tworzących ciąg wielomianów ( j = 0, 1, ..., m) układ równań (7)
przyjmie postać:
a0 s0 + a1s1 + ... + am sm = t0,
üÅ‚
ôÅ‚
a0 s1 + a1s2 + ... + am sm+1 = t1,
ôÅ‚
żł
............................................
ôÅ‚
(9)
ôÅ‚
a0 sm + a1sm+1 + ... + am s2m = tm ,
þÅ‚
gdzie:
k k k
üÅ‚
sk = x0 + x1 + ... + xn ,
ôÅ‚
k k k
tk = x0 y0 + x1 y1 + ... + xn yn ôÅ‚
żł
(10)
ôÅ‚
(k = 0,1, 2, ...).
ôÅ‚
þÅ‚
f (x)
Wielomian aproksymujący daną funkcję w sensie najmniejszych kwadratów
m
j
Qm (x) = a x
" j
(11)
j = 0
f (x)
powinien mieć stopień na tyle wysoki, aby dostatecznie przybliżał funkcję , a jednocześnie mieć
stopień wystarczająco niski, aby wielomian ten wygładzał błędy losowe wynikające np. z pomiarów.
Qm(x)
Jeśli m = n, to wielomian aproksymujący pokrywa się z wielomianem Lagrange a dla układu
x0 , x1, ..., xm
punktów: i wtedy S = 0.
Trudności obliczeniowe związane z aproksymacją średniokwadratową za pomocą wielomianów
wyższych stopni mogą być zmniejszone przy wykorzystaniu wielomianów ortogonalnych.
Õ(x) È(x)
x0, x1,..., xn ,
Dwie dowolne funkcje i nazywamy ortogonalnymi na zbiorze punktów: jeśli
n
Õ(xi ) È (xi ) = 0,
"
(12)
i = 0
przy czym:
n n
2 2
Õ (xi ) > 0, È (xi ) > 0.
" " (13)
i = 0 i = 0
Niech zbiór wielomianów:
G0 (x), G1(x), ..., Gm(x),
(14)
x0 , x1, ..., xn , czyli
będzie danym układem wielomianów ortogonalnych na zbiorze:
n
G (xi ) Gk (xi ) = 0 dla j `" k.
" j (15)
i = 0
Po przedstawieniu wielomianu aproksymującego (11) w postaci kombinacji liniowej wielomianów
układu (14)
Qm(x) = b0 G0(x) + b1G1(x) + ... + bmGm(x)
(16)
odchylenie kwadratowe (6) przyjmuje postać
2
n m
îÅ‚ Å‚Å‚
S = bj G (xi ) - yi śł .
" " j (17)
ïÅ‚
i =0 j = 0
ðÅ‚ ûÅ‚
Podnosząc do kwadratu wyrażenie znajdujące się w nawiasie kwadratowym otrzymujemy
n m m-1 m m
îÅ‚ Å‚Å‚
S = b2 Gj2 (xi ) + 2 bk bj Gk (xi ) Gj (xi ) - 2 bj Gj (xi ) yi + yi2 =
ïÅ‚
"ðÅ‚" " " " śł
j
i=0 j=0 k =0 j=k +1 j=0
ûÅ‚
m n m-1 m n m n n
= b2 Gj2(xi ) + 2 bk bj Gk (xi ) Gj (xi ) - 2 bj Gj (xi ) yi + yi2
" " " " " " " "
j
j=0 i=0 k =0 j=k +1 i=0 j=0 i=0 i=0
i następnie na mocy warunku ortogonalności (15), po wprowadzeniu oznaczeń:
n
2
üÅ‚
s = G (xi ) > 0,
j " j
ôÅ‚
i = 0
ôÅ‚
żł
n (18)
ôÅ‚
c = yi G (xi ),
j " j
ôÅ‚
i = 0
þÅ‚
uzyskujemy
m n
S = (b2 s - 2bj c ) + yi2.
" j j j "
j = 0 i = 0
Uzupełniając wyrażenie znajdujące się w nawiasie pod znakiem pierwszej sumy do pełnego kwadratu
2
ëÅ‚ öÅ‚
c2
j
ìÅ‚bj c j ÷Å‚
b2 s - 2 bj c = s - -
j j j j
ìÅ‚ ÷Å‚
s s
j j
íÅ‚ Å‚Å‚
mamy
2 2
m m n
ëÅ‚ c öÅ‚ c
j j
ìÅ‚
S =
"s j " "
(19)
ìÅ‚bj - s j ÷Å‚ - s j + yi2.
÷Å‚
j =0 j = 0 i=0
íÅ‚ Å‚Å‚
Wynika stąd, że średnie odchylenie kwadratowe S osiąga swą najmniejszą wartość dla współczynników
cj
bj = ( j = 0, 1, ..., m).
(20)
sj
Szukany wielomian aproksymacyjny (16) ma więc postać
m
c
j
Qm (x) = Gj (x).
"
(21)
s
j = 0
j
x0 , x1, ..., xn
m+ 1
Dla układu równoodległych punktów: o stałej odległości h układ wielomianów
ortogonalnych (14) tworzÄ… wielomiany Grama
k
ëÅ‚k öÅ‚ëÅ‚k + söÅ‚
t[s]
ìÅ‚ ÷Å‚ìÅ‚ ÷Å‚
Gk (t) = (-1)s ìÅ‚ ÷Å‚ìÅ‚ ÷Å‚ (k = 0,1, ..., m),
"
(22)
s s
s = 0 n[s]
íÅ‚ Å‚Å‚íÅ‚ Å‚Å‚
gdzie:
pöÅ‚
ëÅ‚
p!
ìÅ‚ ÷Å‚ = ,
q! ( p - q)!
íÅ‚qÅ‚Å‚
q[0] = 1, q[ j] = q (q - 1) ... (q - j + 1) ( j e" 1),
t = (x - x0) h.
Pierwsze wielomiany Grama są następujące:
G0 (t) = 1,
t
G1(t) = 1- 2 ,
n
t t (t -1)
G2 (t) = 1- 6 + 6 ,
n n (n -1)
...........................................
We wzorze aproksymującym (21) opartym na wielomianach Grama współczynnik (18) jest określony
zależnością
[ j+1]
n
(n + j +1)
s = G2 (t) = .
j " j
[ j] (23)
t = 0 (2 j +1) n
....................................
function f(x: Real): Real;
begin
f:=100*x*x*Exp(-10*x);
end;
function Silnia(n: Integer): Real;
begin
if n=0 then Silnia:=1
function G(k,n: Integer; t: Real): Real;
else Silnia:=n*Silnia(n-1);
var
end;
s,p: Integer;
pom: Real;
function Komb(p,q: Integer): Real;
begin
begin
pom:=0; p:=-1;
for s:=0 to k do begin
Komb:=Silnia(p)/Silnia(q)/Silnia(p-q);
if p=-1 then p:=1 else p:=-1;
end;
pom:=pom+p*Komb(k,s)*Komb(k+s,s)*
WCzyn(t,s)/WCzyn(n,s);
function WCzyn(q: Real; j: Integer):
end;
Real;
G:=pom;
var
end;
k: Integer;
....................................
p: Real;
begin
p:=1;
if j>0 then
for k:=0 to j-1 do
p:=p*(q-k);
Wczyn:=p;
end;
.......................................... 2: begin
h:=(b-a)/n; for j:=0 to m do begin
for i:=0 to n do begin s[j]:=0;
x:=a+i*h; for i:=0 to n do
y:=f(x); s[j]:=s[j]+Sqr(G(j,n,i));
xp[i]:=x; end;
yp[i]:=y; for j:=0 to m do begin
xx[i+1]:=x; c[j]:=0;
yy[i+1]:=y; for i:=0 to n do
end; c[j]:=c[j]+yp[i]*G(j,n,i);
Q:=n+2; end;
xx[n+2]:=0; end;
yy[n+2]:=0; end;
case war of odch:=0;
1: begin Writeln(plik,'Wyniki aproksymacji funkcji:');
t[0]:=0; s[0]:=n+1; xd:=xp; Writeln(plik,' i x[i] y[i] blad');
for i:=0 to n do t[0]:=t[0]+yp[i]; for i:=0 to n do begin
for k:=1 to 2*m do begin x:=a+i*h;
t[k]:=0; s[k]:=0; case war of
for i:=0 to n do begin 1: begin
if k<=m then y:=aw[m];
t[k]:=t[k]+xd[i]*yp[i]; for k:=m-1 downto 0 do
s[k]:=s[k]+xd[i]; y:=y*x+aw[k];
xd[i]:=xd[i]*xp[i]; end;
end; 2: begin
end; y:=0;
for i:=0 to m do begin for j:=0 to m do
for j:=0 to m do y:=y+c[j]*G(j,n,i)/s[j];
am[i+1,j+1]:=s[i+j]; end;
am[i+1,m+2]:=t[i]; end;
end; bl:=f(x)-y;
ElimGaussa(m+1,1,am,det); odch:=odch+Sqr(y-yp[i]);
for i:=0 to m do Writeln(plik,i:3,' ',x:13,' ',y:13,' ',bl:13);
aw[i]:=am[i+1,m+2]; Q:=Q+1; xx[Q]:=x; yy[Q]:=y;
end; end;
....................................................
Wyniki obliczeń numerycznych dotyczących aproksymacji średniokwadratowej funkcji
2
1 (1+ 25 x ) za pomocą uogólnionego wielomianu Grama.
Aproksymacja sredniokwadratowa.
Uogolniony wielomian Grama.
Poczatek przedzialu: a = -1.000000E+00
Koniec przedzialu: b = 1.000000E+00
Liczba punktow: n = 100
Stopien wielomianu: m = 15
Wyniki aproksymacji funkcji:
i x[i] y[i] blad
0 -1.000000E+00 2.495641E-02 1.350513E-02
1 -9.800000E-01 6.007476E-02 -2.009075E-02
2 -9.600000E-01 5.404715E-02 -1.244981E-02
3 -9.400000E-01 3.984047E-02 3.468323E-03
4 -9.200000E-01 3.133435E-02 1.379200E-02
5 -9.000000E-01 3.176354E-02 1.529528E-02
6 -8.800000E-01 3.922878E-02 9.887134E-03
7 -8.600000E-01 5.008622E-02 1.222143E-03
8 -8.400000E-01 6.084647E-02 -7.198404E-03
9 -8.200000E-01 6.906730E-02 -1.291907E-02
10 -8.000000E-01 7.360349E-02 -1.477996E-02
11 -7.800000E-01 7.448025E-02 -1.278993E-02
12 -7.600000E-01 7.257915E-02 -7.812310E-03
13 -7.400000E-01 6.926483E-02 -1.191308E-03
14 -7.200000E-01 6.603379E-02 5.599448E-03
15 -7.000000E-01 6.423175E-02 1.123995E-02
16 -6.800000E-01 6.486006E-02 1.475777E-02
...............................................
Odchylenie kwadratowe: 3.300095E-02
KONIEC