## Approximation by polynomials using the least squares method

## Using Matlab function polyfit

>> X = [-1 0 0 1 1 2 4]; % x-coordinates of the points >> Y = [ 5 6 5 7 6 8 11]; % y-coordinates of the points >> step = 0.2; % step-size of the graph of the polynomial >> xx = X(1):step:X(end); % division on x-axis >> deg = 3; % degree of the polynomial >> p = polyfit(X, Y, deg) % computation of coeff. of the polynomial >> yy = polyval(p, xx); % corresponding values of the polynomial >> plot(X, Y, 'ro', xx, yy); >> Yp = polyval(p, X); % approximate values at given points >> D = abs(Y-Yp); % distances from given values >> delta = sqrt( sum( D.^2) ) % quadratic deviation >> max_err = max(D) % maximal error

Experiment with changing the degree of the polynomial from constant (i.e. average value, degree 0), a line (degree 1), up to interpolation polynomial (it can be obtained for distinct x-coordinates only). If you change degree of the polynomial only, you can skip first 4 commands.

## Using a projection

You can use basics of Matlab instead of the function `polyfit`

(the first five and the last six commands are the same as before):

>> X = [-1 0 0 1 1 2 4]; % x-coordinates of the points >> Y = [ 5 6 5 7 6 8 11]; % y-coordinates of the points >> step = 0.2; % step-size of the graph of the polynomial >> xx = X(1):step:X(end); % division on x-axis >> deg = 3; % degree of the polynomial >> %%% computation of coeff. of the polynomial: >> n = length(X); >> A = ones(n,1); % the first column - ones >> for k = 1:deg >> A = [A , (X.^k)' ]; % for every power k, add one column >> end >> %%% A' realizes projection to lin. space >> %%% generated by columns of A: >> C = A' * A; >> b = A' * Y'; >> q = C\b; % computation of coeff. of the polynomial >> p = fliprl(q) % reversing the elements in q >> yy = polyval(p, xx); % corresponding values of the polynomial >> plot(X, Y, 'ro', xx, yy); >> Yp = polyval(p, X); % approximate values at given points >> D = abs(Y-Yp); % distances from given values >> delta = sqrt( sum( D.^2) ) % quadratic deviation >> max_err = max(D) % maximal error