function least_squares(x, y, m)
%least_squares(x, y, m) fits a least-squares polynomial of degree m through
%a set of data where x and y are the coordinates.
%The output to least_squares(x, y, m) is the set of coefficients c0, c1,...,cm
%for the least-square polynomial Pm(x) = c0 + c1*x + c2*x^2 + ...+ cm*x^m
%Hence, there will be m+1 values in the vector c.
%A graph displays the set of points and the fitted polynomial
%This file follows the algorithmic guidelines given in the book of John H. Mathews
%Numerical Methods, for Mathematics, Science, and Engineering, 2nd Ed.
%The polynomial evaluation is an extension of the code given in the divided_diff(x,y,x0)
%m-file.
%Author: Alain G. Kapitho
%Date : Jan. 2006
n = size(x,1);
if n == 1
n = size(x,2);
end
b = zeros(m+1,1);
for i = 1:n
for j = 1:m+1
%right-hand side column vector consisting of sums of
%powers of x multiplied by y's
b(j) = b(j) + y(i)*x(i)^(j-1);
end
end
p = zeros(2*m+1,1);
for i = 1:n
for j = 1:2*m+1
%sums of powers of x
p(j) = p(j) + x(i)^(j-1);
end
end
for i = 1:m+1
for j = 1:m+1
%distributing the sums of powers of x in a matrix A
A(i,j) = p(i+j-1);
end
end
c = A\b
%creating an x-axis and evaluating the least-ssquare polynomial
%this is only needed for data visualization
t = min(x):0.05:max(x);
n = size(t,2);
for i = 1:n
f(i) = c(m+1);
for j = m:-1:1
f(i) = c(j) + f(i)*t(i);
end
end
%data visualization
figure
plot(t,f) %the least_squares polynomial
hold on
plot(x, y, 'r*') %the data points
grid on