Решение систем линейных уравнений
В системе MATLAB для решения систем линейных уравнений предусмотрены знаки операций / и \. Чтобы решить систему линейных уравнений вида
A y = B ,
где A – заданная квадратная матрица размером NxN, a B – заданный вектор – столбец длины N, достаточно применить операцию \ и вычислить выражение A\ B . Операция \ называется левым делением матриц и, будучи примененная к матрицам A и B в виде A \ B , примерно эквивалентна вычислению выражения inv(A)*B . Здесь под inv(A) понимается вычисление матрицы, обратной к матрице A .
Операцию / называют правым делением матрицы. Выражение A/ B примерно соответствует вычислению выражения B*inv(A) . Значит, эта операция позволяет решать системы линейных уравнений вида
y A= B .
% Решаем систему Ax=b
Перепишем систему в векторном виде и введём её матрицы
A = [1 2 3; 1 -3 2; 1 1 1];
b = [7;5;3];
% Проверим систему на невырожденность
rank(A)
>> ans = 3
% Ранг системы полный, система невырождена
% Решим систему с помощью обратной матрицы. x=A^(-1)*b
x = inv(A) * b
| >> | x = | 1.0000 |
| 2.0000 |
% Решим систему с помощью средств MATLAB для решения линейных систем
x = A \ b
| >> | x = | |
% Решим систему Ax=b методом Гаусса
% Для этого, сформируем расширенную систему
A = [1 2 3; 1 -3 2; 1 1 1];
b = [7;5;3];
C = [A b];
% Приведём её к ступенчатому виду, выполнив прямой и обратный ход метода Гаусса
D = rref(C)
% Последний столбец матрицы есть решение
x = D(:,4);
% Проверим его
A*x – b
Построить магический квадрат;
A = magic(6), b = [1;2;3;4;5;6];
http://sl-matlab.ru/upload/resources/EDU%20Conf/pp%20546-607%20Revjakin.pdf
Найти число обусловленности в первой, второй и третьей норме.
Cond(A)
Решить методом Гаусса
R = rref(A)
Программа
2. Например A *x =b (4 уравнения с 4 неизвестными)
Прямой ход метода Гаусса.
L = eye(4);
Первый шаг: L(2; 1) = B(2; 1)=B(1; 1); L(3; 1) = B(3; 1)=B(1; 1);
L(4; 1) = B(4; 1)=B(1; 1); B(2; :) = -L(2; 1) * B(1; :) + B(2; :);
B(3; :) = -L(3; 1)*B(1; :) + B(3; :);
B(4; :) = -L(4; 1) _ B(1; :) + B(4; :):
Второй шаг: L(3; 2) = B(3; 2)/B(2; 2); L(4; 2) = B(4; 2)/B(2; 2);
B(3; :) = -L(3; 2) *B(2; :) + B(3; :);
B(4; :) =-L(4; 2) *B(2; :) + B(4; :).
Третий шаг: L(4; 3) = B(4; 3)/B(3; 3);
B(4; :) =-L(4; 3)* B(3; :) + B(4; :).
Выполнением команды U = B( : ; 1 : 4) завершается прямой ход метода Гаусса, в результате которого получается LU- разложение матрицы A. Проверить можно командой A -L * U:
Обратный ход метода Гаусса: bb = B(: ; 5); x = zeros(4; 1); x(4) = bb(4)/U(4; 4);
x(3) = (bb(3) - U(3; 4) *x(4))/ U(3; 3);
x(2) = (bb(2) - U(2; 4) * x(4) - U(2; 3) *x(3))/U(2; 2);
x(1) = (bb(1) -U(1; 4) *x(4) - U(1; 3) *x(3) - U(1; 2) * x(2))/U(1; 1):
3. О решении систем L *c = b1 и U *x = c.
Решить первую систему L*c = b1 можно, например, следующей последовательность команд:
c = zeros(4; 1); x = c; c(1) = b1(1); c(2) = b1(2) -L(2; 1) * c(1);
c(3) = b1(3) - L(3; 1) * c(1) - L(3; 2) * c(2);
c(4) = b1(4) - L(4; 1) * c(1) - L(4; 2) *c(2) - L(4; 3) * c(3).
Решить вторую систему U * x = cможно следующим образом:
x(4) = c(4)/U(4; 4); x(3) = (c(3) - U(3; 4) * x(4))/U(3; 3);
x(2) = (c(2) -U(2; 4) * x(4) - U(2; 3) * x(3))/U(2; 2);
x(1) = (c(1) - U(1; 4)* x(4) -U(1; 3)*x(3) - U(1; 2) * x(2))/U(1; 1).
4. Стратегия выбора ведущего элемента в методе Гаусса с частичным выбором ведущего элемента заключается в том, что на каждом шаге прямого хода в качестве ведущего элемента берется наибольший (по абсолютной величине) элемент в неприведенной части столбца. Отличие указанного метода от метода Гаусса состоит в том, что в расширенной матрице системы переставляются строки.
E
Циклы
for i = 1 : 3
for j = 1 : 3
E(i; j) = sin(i*j);
end
end
Метод простой итерации
Решите систему Ax=b, методом простой итерации с точностью до 0.001, оценив предварительно необходимое для этого число шагов. В процессе итераций постройте для нее матрицу
V1 =[x1 x2 : : : xk]; где xk – получаемое на каждом шаге приближенное решение системы уравнений. По достижении заданной точности постройте график (командой plot(V 1') ). Ответьте на вопрос: что означает построенный график?
Систему уравнений Ax=b (det(A) не равен нулю) можно представить эквивалентной системой x =(inv(S)*T)*x+inv(S)*b; где A =S-T: Выбрав начальное приближенное решение системы (например, x0=b) и обозначив через x предыдущее приближенное решение системы, а через y – последующее приближенное решение, можно построить итерационный процесс y = (inv(S)*T)*x+inv(S)*b: Итерационный сходится, тогда и только тогда, когда каждое собственное значение матрицы inv(S)*T меньше 1: Скорость сходимости определяется величиной спектрального радиуса r (максимального по модулю собственного значения) матрицы inv(S)*T:
Метод прстой итерации. Матрица S определяется S =D; где D - единичная матрица соответствующего размера умноженная на главную диагональ матрицы A ( в командах MATLAB D = diag (diag (A)) ).
Итерационный процесс метода простой итерации для матрицы 5 порядка можно организовать следующим образом:
V =b; k =1; x =b; xx = zeros(5; 1);
Создать m-file, реализующий следующую последовательность операций:
AA =[A b]; for i =1:5 AA(i; :)=AA(i; :)=AA(i; i); end
xx(1)= - AA(1; 2)*x(2) - AA(1; 3) * x(3)-AA(1; 4)*x(4)- AA(1; 5)*x(5)+AA(1; 6);
xx(2)=-AA(2; 1)*x(1)- AA(2; 3)*x(3)-AA(2; 4)*x(4)-AA(2; 5)*x(5)+ AA(2; 6);
xx(3)=-AA(3; 1)*x(1)-AA(3; 2)*x(2)-AA(3; 4)*x(4)-AA(3; 5)*x(5)+AA(3; 6);
xx(4)=-AA(4; 1)*x(1) - AA(4; 2)*x(2)-AA(4; 3)*x(3)-AA(4; 5)*x(5)++AA(4; 6);
xx(5)=-AA(5; 1)*x(1) - AA(5; 2)*x(2) - A(5; 3)*x(3)-AA(5; 4)*x(4)+AA(5; 6);
x = xx; k =k +1; V =[V x]:
Решение системы A x = b находится повторением приведенных выше процедур.