Исследование систем линейных уравнений неполного ранга .

% SLAE
% The decision of System of the linear algebraic equations
% Решение системы линейных уравнений с минимизацией
% вектора решения по евклидовой норме.
%
% Входные параметры:
%                                      A - матрица коэффициентов системы
%                                      B - вектор столбец решения системы
% Выходные параметры:
%                                      X - вектор решений (A * X = B),
минимизированный по норме
%                                      N - Евклидова норма
%                                      Eps - невязка B - A*X



function [X, N] = SLAE(A, B)

if (nargin < 2)
   error('Необходимо ввести матрицу системы и вектор свободных
коэффициентов');
end;

%Если матрица коэффициентов системы нулевая,
%то вывод сообщения об ошибки и выход
if (A == 0)
   error('Неправильное задание параметров');
end

% m - число строк, n - число столбцов
[m, n] = size(A);

%Проверка на совместность системы
if rank(A) ~= rank([A, B])
   disp('Система не совместна');
   for i = 1 : n
      X(i) = NaN;
   end
   X = X';
   N = 0;
   return
end

% Если высота матрицы а и столбца b не совпадают
% то выдача диагностирующего сообщения
if m ~= length(B)
   error('Высота матрицы A и столбца B не совпадают');
end

% Приведение расширенной матрицы A|B к диагональному виду
A = rref([A, B]);
B = A(:, n + 1);
A = A(:, 1 : n);
%m - число базисных строк
m = rank(A);
%Расчет коэффициентов С(1)..С(n-m), при которых вектор решения Х
%будет минимальным по евклидовой норме. Приравнивая частные производные
%нулю, составляем матрицу коэффициентов D и матрицу свободных коэффициентов
E.
%Соответствующие формулы смотрите в описании к программе.
% i - номер строки, j - номер элемента в строке (номер столбца)
for i=1:(n-m)
   for j=1:(n-m)
      D(i,j) = 0;
      for k=1:m
         D(i,j)=D(i,j)+A(k,i+m)*A(k,j+m);
      end
      if i==j
         D(i,j)=D(i,j)+1;
      end
   end
   E(i)=0;
   for k=1:m
      E(i)=E(i)+B(k)*A(k,i+m);
   end
end


%Транспонирование вектора-строки E в вектор-столбец и
%вычисление коэффициентов С(1)..С(n-m)
E = E';
\ E;

%Вычисление вектора решений в соответствии с найденными коэффициентами
for k = m+1 : n
   X(k) = C(k-m);
end
for k = 1 : m
   X(k) = B(k);
   for j = 1 : (n - m)
      X(k) = X(k) - A(k, j+m)*X(j+m);
   end
end

%Транспонирование вектора-строки X в вектор-столбец
X = X';

%Вывод РЕЗУЛЬТАТОВ
disp('Вектор решения минимизированный по евклидовой норме');
disp(X);
N = norm(X, 'fro');
disp('Евклидова норма вектора решений');
disp(N);
%disp('Невязка Eps =');
%Eps = B - A*X


return

1 2 3 4 5 6 7