Материалы сайта
Это интересно
Исследование систем линейных уравнений неполного ранга
% 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'; C = D \ 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