% Least Squares_Method
% МЕТОД НАИМЕНЬШИХ КВАДРАТОВ ДЛЯ БАЗИСНЫХ ФУНКЦИЙ x^n
% Синтаксис вызова:
%                       [A]=LSM01(F,nV)
% 
%  F - Исходная матрица табличная функция вида:
%
%    x1 x2 x3 x4 .... xM
%    y1 y2 y3 y4 .... yM
%  
%  nV - Длина рапроксимирующего степенного ряда
%
%  Например (подготовка тестовых данных): 
% 
%    x=[0:pi/100:2*pi];
%    y=eval('x+sin(x)-log(x+1)');
%
%    F=[x;y];
%
%  Вызов метода:
%  [A]=LSM01(F,10)
%
%  Результат:
%  Максимальная абсолютная ошибка : 0.00028535
%  A =
%   -0.0003
%    1.0068
%    0.4619
%   -0.4028
%    0.1001
%   -0.0232
%    0.0070
%   -0.0013
%    0.0001
%   -0.0000    
%    
function [A]=lsm01(F,nV)
if isempty(F)
   disp('Матрица [F] не задана'); 
   A=NaN;
   return
end;
% Минимальный входной контроль
if (size(F,1) ~= 2)
   disp('Размеры матрицы [F] заданы неправильно'); 
   A=NaN;
   return
end;
if (size(F,2) < nV)
   disp('Недостаточно точек в [F] для определения коэффициентов для [V]'); 
   A=NaN;
   return
end;
% Подготовка рабочих матриц
C=zeros(nV,nV);
B=zeros(nV,1);
% Построим систему линейных уравнений для аппроксимации
% Для всех строк матрицы 'C'
for row=1:nV
    % Для всех 'x' матрицы 'F'
    B(row,1)=0;
    for k=1:size(F,2)
        x=F(1,k);
        B(row,1)=B(row,1)+ (F(1,k)^(row-1))*F(2,k);  
    end;
    % Для всех колонок матрицы 'C'
    for col=1:nV
        % Для всех 'x' матрицы 'F'
        C(row,col)=0;
        for k=1:size(F,2)
            x=F(1,k);
            C(row,col)=C(row,col)+(F(1,k)^(row-1))*(F(1,k)^(col-1));  
        end;
    end;
end;
% Вычислим коэффициенты для аппроксимирующего выражения
ab=cat(2,C,B);
A = fgauss01(ab);
%  Ниже представлены варианты, которые иногда выдают предупреждения 
%  A=inv(C)*B; % A=C\B;

% Вычислим значения функции по аппроксимирующему выражению
% Для всех 'x' матрицы 'F'
for k=1:size(F,2)
    x=F(1,k);
    y(1,k)=0;
    % Для всех строк матрицы 'V'
    for row=1:nV
        y(1,k)=y(1,k)+A(row,1)*(F(1,k)^(row-1)); 
    end;   
end;

clf;
% Построим график исходной табличной функции
plot(F(1,:),F(2,:));
grid on;
repl=input('Для подолжения нажмите Enter...');
% Построим график по результатам аппроксимации
plot(F(1,:),y);
grid on;
repl=input('Для подолжения нажмите Enter...');


% Построим график абсолютной ошибки
delta=F(2,:)-y;
maxdelta=max(abs(delta));
disp(strcat('Максимальная абсолютная ошибка : ', num2str(maxdelta)));
plot(F(1,:),delta);
grid on;