X =
3	-3	1
-3	5	-2
1	-2	1

Опять таки, поскольку A является симметричной матрицей целых чисел и имеет единичный определитель, то же самое справедливо и для обратной матрицы. С другой стороны, для

B = magic(3)

B =

8	1	6
3	5	7
4	9	2

d = det(B) X = inv(B)

d =

-360

X =

0.1472 -0.1444 0.0639 -0.0611 0.0222 0.1056 -0.0194 0.1889 -0.1028

Внимательное изучение элементов матрицы X, или использование формата rational , показывает, что они являются целыми числами, разделенными на 360.

Если матрица A является квадратной и несингулярной, то, пренебрегая ошибками округления, выражение X = inv(A)*B теоретически означает то же, что и X = A\B , а Y = B*inv(A) теоретически есть то же, что и Y = B/A. Однако вычисления включающие операторы \ и / более предпочтительны, поскольку требуют меньше рабочего времени, меньшей памяти и имеют лучшие свойства с точки зрения определения ошибок.

Псевдообратные матрицы

Прямоугольные матрицы не имеют детерминантов и обратных матриц. Для таких матриц по крайней мере одно из уравнений AX = I или XA = I не имеет решения. Частично данный пробел восполняется так называемой псевдообратной матрицей Мура-Пенроуза, или просто псевдообратной матрицей, которая вычисляется при помощи функцииpinv. На практике необходимость в этой операции встречается довольно редко. Желающие могут всегда обратиться к соответствующим справочным пособиям.

Степени матриц и матричные экспоненты

Положительные целые степени

Если А есть некоторая квадратная матрица, а р – положительное целое число, то A^p эквивалентно умножению A на себя р раз.

X = A^2

35

X =

3	6	10
6	14	25
10	25	46

Отрицательные и дробные степени

Если А является квадратной и невырожденной, то A^(-p) эквивалентно умножению inv(A) на себя p раз.

Y = B^(-3)

Y =

0.0053 -0.0068 0.0018 -0.0034 0.0001 0.0036 -0.0016 0.0070 -0.0051

Дробные степени, например A^(2/3), также допускаются; результаты при этом зависят от распределения собственных значений матрицы А.

Поэлементное возведение в степень

Оператор .^ (с точкой !) осуществляет поэлементное возведение в степень. Например,

X = A.^2

A =
1	1	1
1	4	9
1	9	36

Вычисление корня квадратного из матрицы и матричной экспоненты

Для невырожденных квадратных матриц А функция sqrtm вычисляет главное значение квадратного корня , т.е. если X = sqrtm(A) , то X*X = A . Буква m в sqrtm означает, что выполняется матричная операция. Это отличает данную функцию отsqrt(A), которая, подобно A.^(1/2) (обратите внимание на точку!), выполняет операцию извленчения корня поэлемен-

тно.

Система обыкновенных линейных дифференциальных уравнений первого порядка может быть записана в виде

dx/dt = Ax

где x = x(t) есть векторная функция от t, а A есть постоянная матрица не зависящая от t. Решение данной системы может быть выражено в виде матричной экспоненты.

x(t) = ℮Atx(0)

Функция expm(A) вычисляет матричную экспоненту. Рассмотрим пример системы дифференциальных уравнений со следующей 3х3 матрицей коэффициентов

A =

36

0	-6	-1
6	2	-16
-5	20	-10

и начальными условиями x(0)

x0 = [ 1 1 1]’.

Использование матричной экспоненты для вычисления решения дифференциального уравнения в 101 точке с шагом 0.01 на интервале 0 ≤ t ≤ 1 записывается в виде

X = [ ];

for t = 0 : 0.01 : 1

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

end

Трехмерный график решения в фазовом пространстве может быть получен при помощи специальной функции

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

Решение имеет вид спиральной функции сходящейся к началу координат (см. рис. ниже). Такое решение обусловлено комплексными собственными значениями матрицы коэффициентов А.

Собственные значения и собственные векторы

Собственным значением и собственным вектором квадратной матрицы А называются скаляр λ и вектор v, удовлетворяющие условию

Av = λv

37

Диагональная декомпозиция

Имея диагональную матрицу Λ, составленную из собственных значений λ матрицы А и матрицу V , составленную из соответствующих собственных векторов v, можно записать

AV = VΛ

Если матрица V несингулярная, на основании данного выражения получаем спектральное разложение матрицы А

А = VΛV-1

Неплохой пример использования спектрального разложения дает рассмотренная выше матрица коэффициентов линейного дифференциального уравнения. Ввод выражения

lambda = eig(A)

дает следующий вектор-столбец собственных значений(два из них являются комплексносопряженными)

lambda =

-3.0710 -2.4645 + 17.6008i -2.4645 - 17.6008i

Действительные части всех собственных значения являются отрицательными, что обеспечивает устойчивость процессов в системе. Ненулевые мнимые части комплексно-сопряженных собственных значений обуславливают колебательный характер переходных процессов.

При двух выходных аргументах, функция eig вычисляет также собственные векторы и выдает собственные значения в виде диагональной матрицы

.

	[V,D] = eig(A)
V =
-0.8326	0.2003 - 0.1394i	0.2003 + 0.1394i
-0.3553 -0.2110 - 0.6447i		-0.2110 + 0.6447i
-0.4248	-0.6930		-0.6930
D =
-3.0710	0		0
0	-2.4645+17.6008i		0
0	0	-2.4645-17.6008i

Первый собственный вектор (первый столбец матрицы V) является действительным, а два других являются комплексно-сопряженными. Все три вектора являются нормализованными по длине, т.е. их Евклидова норма norm(v,2), равна единице.

Матрица V*D*inv(V), которая в более сжатой форме может быть записана как V*D/V, равна, в пределах погрешностей округления, матрице А. Аналогично, inv(V)*A*V, или V\A*V, равна, в пределах погрешностей округления, матрице D.

38