Решение систем линейных уравнений

В системе 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 находится повторением приведенных выше процедур.