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

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

%Тестирование функции решения систем линейных алгебраических уравнений SLAE

%Пример 1

% Матрица коэффициентов при неизвестных

A = [ 1 -3 6 -5 0; 4 2 1 10 2; 2 0 -9 1 6 ]

% Матрица свободных членов

B = [ 3; 5; 7 ]

% --== 1 ==--

disp('- - = = 1 = = - -');

disp('Стандартное решение посредствам системы MatLab X = A\B');

X = A\B;

disp('X = ');

disp(X);

disp('Невязка Eps = ')

disp(B - A*X);

disp('Евклидова норма N = ')

disp(norm(X, 'fro'));

% --== 2==--

disp('- - = = 2 = = - -');

disp('Решение MatLab c первоначальной диагонализацией по методу Гауса');

% Приведение расширенной матрицы A|B к диагональному виду

[m, n] = size(A);

A = rref([A, B]);

B = A(:, n + 1);

A = A(:, 1 : n);

X = A\B

disp('Невязка Eps = ');

disp(B - A*X);

disp('Евклидова норма N = ');

disp(norm(X, 'fro'));

% --== 3 ==--

disp('- - = = 3 = = - -');

disp('Решение системы функцией SLAE');

% Повторный ввод параметров

A = [ 1 -3 6 -5 0; 4 2 1 10 2; 2 0 -9 1 6 ];

B = [ 3; 5; 7 ];

[X, N3] = SLAE(A, B);

disp('Невязка Eps = ');

disp(B - A*X);