LABORATORIUM

METODY NUMERYCZNE

Aproksymacja średniokwadratowa

Prowadzący: dr inż. Tomasz Kozłowski

Autor: Przemysław Wojciechowski

Treść zadania:

Dla funkcji f(x) = 2cos(3x) w punktach{1,2,5} wyprowadź ręcznie współczynniki funkcji: f(x) = ax² + bx + c stosując metodę aproksymacji średnio kwadratowej i porównaj
z wyliczeniami numerycznymi stosując odpowiednie skrypty Matlaba. Następnie wyznacz postać wielomianu interpolacyjnego metodą Newtona używając skrypty Matlaba. Wykreśl
za pomocą komendy plot błędy lokalne E₁(x) = abs(y(x)−f(x)) oraz
E₂(x) = abs(y(x)−P₂(x)) i skomentuj rezultat. Wykreśl podane punkty {x_i,y(x_i)} przebieg funkcji y(x), f(x) oraz przebieg wyznaczonego wielomianu interpolacyjnego używając komenty plot. Wydruki wykresów (wraz z kodem) dołącz do sprawozdania

Wstęp teoretyczny:

Interpolacja polega na znajdowaniu krzywej przechodzącej przez wszystkie węzły. Zdarzają się jednak sytuacje, w których dane te mogą być obarczone błędem, np. gdy są one odczytywane z pomiarów. W takim wypadku bardziej interesujące jest zachowanie ogólnego charakteru krzywej, która te punkty przybliża niż przeprowadzanie jej przez wszystkie węzły. Do tego zadania wykorzystywana jest aproksymacja metod¡ najmniejszych kwadratów.
W aproksymacji metodą najmniejszych kwadratów mamy z góry zadaną bazę funkcji
za pomocą których chcemy przybliżyć funkcję poszukiwaną postaci:

F(x) =  c₁f₁(x) +  c₂f₂(x) +  . . .  +  c_nf_n(x)

dla zadanego zbioru punktów (x_j,  y_j), j = 1,  . . . , m. Jak wynika z tego wzoru, poszukiwana funkcja F(x) jest liniową kombinacją funkcji bazowych f_i(x), i = 1,  . . . ,  n. Funkcje bazowe mogą być dowolne, dopóki są liniowo niezależne. Niewiadomymi w tak postawionym zadaniu są współczynniki c_i. Ponieważ, w ogólnym przypadku, tak zadana funkcja poszukiwana nie będzie pokrywała się ze wszystkimi węzłami. W zależności od doboru współczynników c_i może lepiej przybliżać pewne węzły a gorzej inne. Należy zatem dobrać odpowiednie kryterium porównawcze. Możemy obliczyć pionową odległość funkcji
od każdego z węzłów:

r_j = y_j − F(x_j)

Ale ponieważ tak zapisane r_j może być zarówno dodatnie jak i ujemne, to proste zsumowanie tych wartości nie nadaje się do porównywania różnych krzywych. Lepszym wyborem jest użycie r_j² wtedy:

$$p = \sum_{j = 1}^{m}r_{j}^{2}$$

jest wartością, której minimum będzie określało krzywą najlepiej dopasowaną do węzłów.

Dopasowanie metod¡ najmniejszych kwadratów dla dowolnej bazy funkcji

W ogólnym przypadku doboru bazy funkcji f_i(x), i = 1,  . . . ,  n będziemy mieli:

y = F(x)=c₁f₁(x) +  c₂f₂(x) +  . . .  +  c_nf_n(x)

Niewiadomymi będą tu współczynniki c_i. Celem jest znalezienie takich współczynników c_i, …,  c_n, które minimalizują kwadraty pionowych odległości funkcji F(x) od węzłów (x_j, y_j), j = 1,  . . . , m. Dla tak postawionego zadania możemy zapisać układ równań:

c₁f₁(x₁) +  c₂f₂(x₁) +  . . .  +  c_nf_n(x₁) = y₁

c₁f₁(x₂) +  c₂f₂(x₂) +  . . .  +  c_nf_n(x₂) = y₂

. . .

c₁f₁(x_n) +  c₂f₂(x_n) +  . . .  +  c_nf_n(x_n) = y_n

lub w postaci macierzowej:

Ac = y

gdzie:

$$\left\lbrack \begin{matrix} f_{1}(x_{1}) & f_{2}(x_{1}) & \ldots \\ \ldots & \ldots & \ldots \\ f_{1}(x_{m}) & f_{2}(x_{m}) & \ldots \\ \end{matrix}\text{\ \ \ \ }\begin{matrix} f_{n}(x_{1}) \\ \ldots \\ f_{n}(x_{m}) \\ \end{matrix} \right\rbrack \bullet \begin{bmatrix} c_{1} \\ \ldots \\ c_{n} \\ \end{bmatrix} = \begin{bmatrix} y_{1} \\ \ldots \\ y_{m} \\ \end{bmatrix}$$

Dochodzimy do podobnej postaci układu jak w poprzednim podrozdziale różniącym się tylko definicją macierzy A. Naszym zadaniem jest zminimalizowanie wartości ρ = ||y−Ac||₂². Postępując w przedstawiony wcześniej sposób otrzymujemy równania normalne postaci:

(A^TA)c = A^Ty

Wyznaczając wektor c z tego równania otrzymujemy nieznane współczynniki.

Definicja iloczynu skalarnego 

Iloczyn skalarny dwóch wektorów  i  jest to liczba równa iloczynowi modułów (długości) tych wektorów i cosinusa kąta między nimi w przypadku, gdy są to wektory niezerowe i równa zeru, gdy jeden lub drugi wektor jest wektorem zerowym. Iloczyn skalarny oznaczamy następująco:

$$\left( \overset{\overline{}}{a},\overset{\overline{}}{b} \right) = \overset{\overline{}}{a} \bullet \overset{\overline{}}{b} \bullet cos(\overset{\overline{}}{a},\overset{\overline{}}{b})$$

Przykład obliczeniowy:

Tab.1 Tabela z wartościami punktów

i	0	1	2
x	1	2	5
y(x)	2	1,99	1,93

Wzór wielomianu:

Q(x) = cφ₀(x) + bφ₁(x) + aφ₂(x)

gdzie:
φ₀(x) = 1, φ₁(x) = x, φ₂(x) = x²

Rozpisując macierz otrzymujemy:

$$\begin{bmatrix} \left( \varphi_{0},\varphi_{0} \right) & \left( \varphi_{0}{,\varphi}_{1} \right) & \left( {\varphi_{0},\varphi}_{2} \right) \\ \left( \varphi_{1},\varphi_{0} \right) & \left( \varphi_{1},\varphi_{1} \right) & \left( \varphi_{1}{,\varphi}_{2} \right) \\ \left( \varphi_{2},\varphi_{0} \right) & \left( \varphi_{2},\varphi_{1} \right) & \left( \varphi_{2},\varphi_{2} \right) \\ \end{bmatrix} \bullet \begin{bmatrix} a \\ b \\ c \\ \end{bmatrix} = \begin{bmatrix} \left( f,\varphi_{0} \right) \\ \left( f,\varphi_{1} \right) \\ \left( f,\varphi_{2} \right) \\ \end{bmatrix}$$

Obliczam elementy macierzy:

$$\left( f,\varphi_{0} \right) = \sum_{i = 0}^{2}{f\left( x_{i} \right) \bullet}\varphi_{0}\left( x_{i} \right) = 2 \bullet 1 + 1,99 \bullet 1 + 1,93 \bullet 1 = 5,92$$

$$\left( f,\varphi_{1} \right) = \sum_{i = 0}^{2}{f\left( x_{i} \right) \bullet}\varphi_{1}\left( x_{i} \right)2 \bullet 1 + 1,99 \bullet 2 + 1,93 \bullet 5 = 15,63$$

$$\left( f,\varphi_{2} \right) = \sum_{i = 0}^{2}{f\left( x_{i} \right) \bullet}\varphi_{2}\left( x_{i} \right) = 2 \bullet 1^{2} + 1,99 \bullet 2^{2} + 1,93 \bullet 5^{2} = 58,21$$

$$\left( \varphi_{0},\varphi_{0} \right) = \sum_{i = 0}^{2}{\varphi_{0}\left( x_{i} \right) \bullet}\varphi_{0}\left( x_{i} \right) = 3$$

$$\left( \varphi_{0},\varphi_{1} \right) = \sum_{i = 0}^{2}{\varphi_{1}\left( x_{i} \right) \bullet}\varphi_{0}\left( x_{i} \right) = 8$$

$$\left( \varphi_{0},\varphi_{2} \right) = \sum_{i = 0}^{2}{\varphi_{2}\left( x_{i} \right) \bullet}\varphi_{0}\left( x_{i} \right) = 30$$

$$\left( \varphi_{1},\varphi_{1} \right) = \sum_{i = 0}^{2}{\varphi_{1}\left( x_{i} \right) \bullet}\varphi_{1}\left( x_{i} \right) = 30$$

$$\left( \varphi_{1},\varphi_{2} \right) = \sum_{i = 0}^{2}{\varphi_{2}\left( x_{i} \right) \bullet}\varphi_{1}\left( x_{i} \right) = 134$$

$$\left( \varphi_{2},\varphi_{2} \right) = \sum_{i = 0}^{2}{\varphi_{2}\left( x_{i} \right) \bullet}\varphi_{2}\left( x_{i} \right) = 642$$

Po podstawieniu:

$$\begin{bmatrix} 3 & 8 & 30 \\ 8 & 30 & 134 \\ 30 & 134 & 642 \\ \end{bmatrix} \bullet \begin{bmatrix} a \\ b \\ c \\ \end{bmatrix} = \begin{bmatrix} 5,92 \\ 15,63 \\ 58,21 \\ \end{bmatrix}$$

Wielomian interpolacyjny będzie miał następującą postać:

Q(x)=−1,27x^2+7,9x−8,7

Kod źródłowy:

clc; clear all; close all; format compact

%% Dane

xi=[1 2 5];

q=inline('2*cos(3*xi)');

yi=q(xi);

%% Aproksymacja średnio kwadratowa

n=length(xi);

ai=polyfit(xi,yi,n-1);

a=min(xi); b=max(xi);

x=linspace(a,b,1000000);

Ap=polyval(ai,x);

disp(sprintf('Wspolczynniki aproksymacji średniokwadratowej'))

ai

E1=abs(q(x)-Ap);

%% Newton

[cn,DD] = newtonp(xi,yi);

P2 = polyval(cn,x);

disp(sprintf('Wspolczynniki wielomianu Newtona'))

cn

E2=abs(q(x)-P2);

%% Rysownie wykresu

grid on

plot(xi,yi,'mo',x,q(x),'-k',x,Ap,'.g',x,E1,'xb',x,P2,'m',x,E2,'ok')

legend('punkty zadane','funkcja bazowa','aproksymacja średnio kwadratowa','błąd aprok.','wielomian Newtona','błąd w. Newtona')

Skrypt obliczenia wielomianu Newtona:

function [cn,DD] = newtonp(x,y)

N = length(x) - 1;

DD = zeros(N+1,N+1);

DD(1:N+1,1)=y';

for k=2:N+1

for m=1:N+2-k

DD(m,k) = (DD(m+1,k-1) - DD(m,k-1))/(x(m+k-1)-x(m));

end

end

r = DD(1,:);

cn = r(N+1);

for k=N:-1:1

cn = [cn r(k)] - [0 cn*x(k)];

end

end

Wyniki obliczeń:

Tab.2 Wyniki obliczeń

Wielomian interpolacyjny wyliczony ręcznie metodą aproksymacji średniokwadratowej	Q(x) = −1, 27x^2 + 7, 9x − 8, 7
Wielomian interpolacyjny wyliczony w programie Matlab	M3(x) = −1, 2617x^2 + 7, 6855x − 8, 4038

Wnioski:

Wyliczone wielomiany interpolacyjne są bardzo do siebie zbliżone. Różnice między nimi mogą wynikać z błędów przybliżeń zastosowanych w metodzie obliczania metodą aproksymacji średniokwadratowej obliczanej ręcznie, w celu ułatwienia obliczeń.