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

Помимо прямых методов решения систем линейных алгебраических уравнений MATLAB обладает мощными современными итерационными методами решения больших систем уравнений. К итерационным методам относятся те, решение которых получается в ходе ряда шагов – итераций, постепенно ведущих к получению результата с заданной погрешностью.

В качестве примера рассмотрим двунаправленный метод сопряженных градиентов, который реализуется с помощью функции bigcg.

 

Построим матрицу системы с три диагональной структурой и размером 10 x 10:

 

n = 10;

on = ones(n, 1);

A = spdiags([-2*on 4*on -on], -1: 1, n, n);

 

В качестве элементов вектора правой части возьмем сумму элементов матрицы A в соответствующих строках:

 

b = sum(A, 2);

 

Для решения системы Ax=b применим один из одиннадцати существующих форматов вызова функции bigcg:

x = bicg(A, b);

bicg converged at iteration 10 to a solution with relative residual 8.6e-014

 

В результате система выдает сообщение о сходимости итерационного процесса и приводит относительную погрешность.

 

Обратная матрица и определитель.

Определитель квадратной матрицы находится с помощью функции det:

 

Det(A)

ans =

Обратная матрица находится с помощью функции inv:

 

Inv(A)

ans =

3 -3 1

-3 5 -2

1 -2 1

Факторизация Холецкого.

Если матрица системы является симметричной (или эрмитовой) и положительно определенной, то ее можно представить в виде произведения двух треугольных матриц: A = R’R, где R – верхняя треугольная матрица, R’ – транспонированная. Факторизация Холецкого выполняется с помощью функции chol:

 

R = chol(A)

R =

1 1 1

0 1 2

0 0 1

Легко проверить, что произведение R’R дает исходную матрицу Паскаля. Факторизация Холецкого приводит систему Ax=u к виду R’Rx=u. Поскольку оператор ‘ \ ’ распознает системы с треугольными матрицами, приведенную систему можно решить очень быстро:

 

x = R\(R'\u)

x =

-12

LU факторизация.

LU факторизация, или гауссово исключение, выражает любую квадратную матрицу A как произведение перестановки нижней треугольной матрицы L и верхней треугольной матрицы U, где матрица L имеет единичную главную диагональ. Перестановки важны как с теоретической, так и с практической точки зрения. Так, матрицу невозможно представить в виде произведения двух треугольных матриц без перестановки двух строк. С другой стороны, хотя матрицу можно представить в виде произведения двух треугольных матриц, при малом значении é элементы матриц-сомножителей будут очень большими, что приведет к большой численной погрешности. Для выполнения разложения служит функция lu:

 

· [ L, U ] = lu(X) возвращает верхнюю треугольную матрицу U и нижнюю треугольную матрицу L (точнее, произведение нижней треугольной матрицы и матрицы перестановок), так что X = L*U;

· [ L, U, P] = l u(X) возвращает верхнюю треугольную матрицу U, нижнюю треугольную матрицу L и матрицу перестановок P, так что L*U = P*X;

· B = lu(X) возвращает матрицу B такую что, B(i, j) = U(i, j) для всех индексов jSi и B(i, j) = L(i, j) для всех индексов j< i. Поскольку диагональные элементы матрицы L равны единице, их можно не хранить.

Рассмотрим факторизацию магической матрицы B:

 

[L U] = lu(B)

L =

1.0000 0 0

0.3750 0.5441 1.0000

0.5000 1.0000 0

U =

8.0000 1.0000 6.0000

0 8.5000 -1.0000

0 0 5.2941

Отметим, что строки матрицы L переставлены местами. Легко проверить, что B = L*U.

LU факторизация матрицы B позволяет очень быстро решить систему Bx=u:

 

x = U\(L\u)

x =

0.5528

0.2611

-0.2806

 

Получим представление матрицы перестановок P:

 

[L, U, P] = lu(B)

L =

1.0000 0 0

0.5000 1.0000 0

0.3750 0.5441 1.0000

U =

8.0000 1.0000 6.0000

0 8.5000 -1.0000

0 0 5.2941

P =

1 0 0

0 0 1

0 1 0

 

Строки матрицы L стоят на своих местах. Легко проверить, что L*U = P*B. Определитель и обратную матрицу из LU разложения можно найти по формулам:

det(B) = det(L)*det(U) и inv(B) = inv(U)*inv(L).

QR факторизация.

Функция qr выполняет QR разложение матрицы. Эта операция полезна как для квадратных, так и для прямоугольных матриц. Исходная матрица представляется в виде произведения действительной ортонормальной или комплексной унитарной матрицы Q и верхней треугольной матрицы R. Функция используется в следующих формах:

 

· [Q, R] = qr(A) вычисляет верхнюю треугольную матрицу R того же размера, что и A, и унитарную матрицу Q, так что A = Q*R;

· [Q, R, E] = qr(A) вычисляет матрицу перестановок E, верхнюю треугольную матрицу R c убывающими по модулю диагональными элементами и унитарную матрицу Q, так что A*E = Q*R;

· [Q, R] = qr(A, 0) и [Q, R, E] = qr(A, 0) вычисляют экономное разложение, в котором E – вектор перестановок такой, что Q*R = A(:, E). Вектор E выбирается так, чтобы диагональные элементы матрицы R убывали по модулю;

· X = qr(A) возвращает результат из пакета LAPACK так, что R есть triu(X) – верхняя треугольная часть матрицы X. Матрица R позволяет избежать потери точности при вычислении A'*A с помощью разложения Холецкого: A'*A = R'*R.

 

В качестве примера рассмотрим QR разложение прямоугольной матрицы C:

 

C = fix(10*rand(3, 2))

C =

6 8

9 6

8 1

[Q, R] = qr(C)

Q =

-0.4460 0.7450 0.4961

-0.6690 0.0908 -0.7377

-0.5946 -0.6609 0.4579

R =

-13.4536 -8.1762

0 5.8437

0 0

Во многих случаях последние m-n столбцы матрицы Q можно отбросить, так как они умножаются на последние нулевые строки матрицы R. Используя экономный формат функции разложения, получим:

[Q, R] = qr(C, 0)

Q =

-0.4460 0.7450

-0.6690 0.0908

-0.5946 -0.6609

R =

-13.4536 -8.1762

0 5.8437

Решение системы Cx = u с факторизованной матрицей находится в два этапа:

 

y = Q'*u;

x = R\y

x =

0.3590

-0.0544

Поскольку матрица системы является переопределенной, то найденное решение является наилучшим среднеквадратичным приближением. Отметим, что оператор x = C\u даст такой же результат.

 

Матричная экспонента.

Матричная экспонента expm(X) возводит число e в матричную степень X. Рассмотрим в качестве примера задачу Коши для системы обыкновенных дифференциальных уравнений первого порядка с постоянными коэффициентами:

 

d x/ dt= Ax, x(0)=u,

 

здесь x = x (t) – векторная функция, A – матрица коэффициентов, u – вектор начальных условий. Решение задачи можно выразить через матричную экспоненту x( t )=e ^{t A} x(0).

Положим:

 

A=[0 -6 -1; 6 2 -16; -5 20 -10];

u=[1; 1; 1];

 

Определим решение системы дифференциальных уравнений в узлах равномерной сетки с шагом 0.1 на отрезке [0, 1] и построим трехмерный график фазовой траектории:

 

X = [];

for t = 0:.01: 1

X = [X expm(t*A)*u];

End

plot3(X(1,: ), X(2,: ), X(3,: ), '-o')

Box on

 

 

⇐ Предыдущая3456789101112Следующая ⇒

Поделиться: