|
|
Урок №11. Матричные операции линейной алгебры
Линейная алгебра — область, в которой наиболее часто используются векторы и матрицы. Наряду с операциями общего характера, рассмотренными выше, применятся функции, решающие наиболее характерные задачи линейной алгебры. Они и рассмотрены в данном уроке. Вычисление нормы и чисел обусловленности матрицы Для понимания всего нижеизложенного материала необходимо учесть, что нормы матриц в MATLAB отличаются от норм векторов. Пусть А —матрица. Тогда n=norm(A) эквивалентно п=погп(А,2) и возвращает вторую норму, т. е. самое большое сингулярное число А. Функция n=norm(A, 1) возвращает первую норму, т. е. самую большую из сумм абсолютных значений элементов матрицы по столбцам. Норма неопределенности n=norm(A, inf) возвращает самую большую из сумм абсолютных значений элементов матрицы по рядам. Норма Фробениуса (Frobenius) norm(A, 'fro') = sqrt(sum(diag(A'A))). Пример: » A=[2.3.1:1.9.4:2.6.7] A = 2 3 1 1 9 4 2 6 7 » norm(A.l) ans = 18 Числа обусловленности матрицы определяют чувствительность решения системы линейных уравнений к погрешностям исходных данных. Следующие функции позволяют найти числа обусловленности матриц.
Пример: » d=cond(hilb(4)) d = 1.5514е+004
Большие числа обусловленности означают, что матрица А близка к матрице с кратными собственными значениями. Пример: » d=condeig(rand(4)) d = 1.0766 1.2298 1.5862 1.7540
Пример: » s=rcond(hilb(4)) s = 4.6461е-005 Для нахождения определителя (детерминанта) и ранга матриц в MATLAB имеются следующие функции:
Пример: » А=[2,3,6:1.8.4;3.6,7] А =
» det(A) ans =
-29 [L.U>lu(A): s=det(L): d=s*prod(diag(U)). Ранг матрицы определяется количеством сингулярных чисел, превышающих порог tol=max(size(A))*nprm(A)*eps. При этом используется следующий алгоритм: s=svd(A);tol=max(size(A))*npnri(A)*eps:r=sum(s>tol): Для вычисления ранга используется функция rank:
» rank(hilbdl)) ans = 10 Норма вектора — скаляр, дающий представление о величине элементов вектора. Нужно различать норму матрицы и норму вектора. Функция norm определяет, является ли ее аргументом (входным аргументом в терминологии MATLAB) вектор или матрица, и измеряет несколько различных типов норм векторов:
Определение ортонормированного базиса матрицы Вычисление ортонормированного базиса матрицы обеспечивают нижеприведенные функции:
Пример: » А=[2 4 6:9 8 2:12 23 43] А =
» B=orth(A) В= 0.1453 -0.0414-0.9885 0.1522 -0.98630.0637
0.9776 0.1597 0.1371
Пример: » null(hilb(11)) ans = 0.0000 -0.0000 0.0009 -0.0099 0.0593 -0.2101 0.4606 -0.6318 0.5276 -0.2453 0.0487 Функции приведения матрицы к треугольной форме Треугольной называется квадратная матрица А, если при l>k (верхняя треугольная матрица) или при к>1(нижняя треугольная матрица) элементы матрицы A(l,k) равны нулю. В строго треугольной матрице нули находятся и на главной диагонали. В линейной алгебре часто используется приведение матриц к той или иной треугольной форме. Оно реализуется следующими функциями:
Примеры:
Определение угла между двумя подпространствами Угол между двумя подпространствами вычисляет функция subsрасе:
Пример: » Н = hadamard(20);A = Н(:.2:4);В = Н(:.5:8): » subspace(A.B) ans = 1.5708 След матрицы А — это сумма ее диагональных элементов. Он вычисляется функцией trace:
» а=[2.3.4:5.6,7;8.9,1] а = 234 5.6 7 8 9 1 » trace(a) ans = 9 Разложение Холецкого — известный прием матричных вычислений. Функция chol находит это разложение для действительных и комплексных эрмитовых матриц.
Пример: » c=chol(pascal(4)) с = 1 1 1 1 0 1 2 3 0 0 1 3 0 0 0 1 Обращение матриц — функции inv, pinv Обращение матриц — одна из наиболее распространенных операций матричного анализа. Обратной называют матрицу, получаемую в результате деления единичной матрицы Е на исходную матрицу X. Таким образом, Х^-1=Е/Х. Следующие функции обеспечивают реализацию данной операции:
Пример: » inv(rand(4,4))
ans =
2.2631 -2.3495-0.4696-0.6631
-0.76201.2122 1.7041 -1.2146 -2.04081.4228 1.5538 1.3730 1.3075 -0.0183-2.54830.6344 На практике вычисление явной обратной матрицы не так уж необходимо. Чаще операцию обращения применяют при решении системы линейных уравнений вида Ах=b. Один из путей решения этой системы — вычисление x=inv(A)*b. Но лучшим с точки зрения минимизации времени расчета и повышения точности вычислений является использование оператора матричного деления х=А\b. Эта операция использует метод исключения Гаусса без явного формирования обратной матрицы.
Пример: » pinv(rand(4,3)) ans = -1.3907-0.0485-0.24931.8640 -0.87751.1636 0.6605 -0.0034
2.0906 -0.5921-0.2749-0.5987 Так называемые LU- и QR-разложения реализуются следующими матричными функциями: Функция 1и выражает любую квадратную матрицу X как произведение двух треугольных матриц, одна из которых (возможно, с перестановками) — нижняя треугольная матрица, а другая — верхняя треугольная матрица[В MATLAB 6 аргументом (входным аргументом) функции lu может быть и полная прямоугольная матрица. — Примеч. ред.]. Иногда эту операцию называют LR-разложением. Для выполнения этой операции служит следующая функция:
Функция qr выполняет QR-разложепие матрицы. Эта операция полезна для квадратных и треугольных матриц. Она выполняет QR-разложение, вычисляя произведение унитарной [Квадратная матрица с комплексными элементами, обладающая тем свойством, что обратная матрица ее комплексно сопряженной матрицы равна транспонированной, т. е. (А*)''-А'. — Примеч. ред.] матрицы и верхней треугольной матрицы. Функция используется в следующих формах: [Квадратная матрица с комплексными элементами, обладающая тем свойством, что обратная матрица ее комплексно сопряженной матрицы равна транспонированной, т. е. (А*)''-А'. — Примеч. ред.]
» C=rand(5.4) С= 0.8381 0.5028 0.1934 0.6979 0.0196 0.7095 0.6822 0.3784 0.6813 0.4289 0.3028 0.8600 0.3795 0.3046 0.5417 0.8537 0.8318 0.1897 0.1509 0.5936 » [Q.R]=qr(C) Q= -0.5922-0.11140.5197 0.0743 -0.6011 -0.0139-0.9278 -0.0011 -0.34480.1420 -0.4814-0.11730.0699 0.5940 0.6299 -0.2681-0.1525-0.82680.2632 -0.3898 -0.58770.2997 -0.2036-0.67340.2643 R = -1.4152 -0.7072 -0.5037 -1.4103 0 -0.7541 -0.7274 -0.4819 0 0 -0.3577 -0.4043
0 0
0 0.2573 0 0 0 0
Примеры: » C=rand(3.3) С = 0.0164 0.0576 0.7176 0.1901 0.3676 0.6927 0.5869 0.6315 0.0841 » [Q.R]=qr(C) Q= -0.0265-0.2416-0.9700 -0.3080-0.92120.2378 -0.95100.3051 -0.0500 R = -0.6171-0.7153-0.3123 0 -0.1599-0.7858 0 0 -0.5356 » [Q1.R1]=qrdelete(Q.R.2) Q1 = -0.02650.7459 0.6655 -0.30800.6272 -0.7153 -0.9510-0.22390.2131 R1 = -0.6171-0.3123 0 0.9510 0 0
Примеры: » C=rand(3.3) С = 0.1210 0.8928 0.8656 0.4508 0.2731 0.2324 0.7159 0.2548 0.8049 » [Q,R]-qr(c) Q = -0.14160.9835 0.1126 -0.52750.0213 -0.8493 -0.8377-0.17970.5157 R = -0.8546-0.4839-0.9194 0 0.8381 0.7116 0 0 0.3152 » x=[0.5.-0.3.0.2];[Q2.R2]=qrinsert(Q.R.2,x') Q2 = -0.14160.7995 -0.5838 -0.5275-0.5600-0.6389 -0.83770.2174 0.5010 R2 = -0.8546-0.0801-0.4839-0.9194 0 0.6112 0.6163 0.7369 0 0 -0.5681-0.2505 Вычисление собственных значений и сингулярных чисел Во многих областях математики и прикладных наук большое значение имеют средства для вычисления собственных значений (собственных чисел, характеристических чисел, решений векового уравнения) матриц, принадлежащих им векторов и сингулярных чисел. В новой версии MATLAB собственные вектора нормализуются, иначе, чем в предыдущих. Основной критерий: либо V'V=I, либо V'BV=I, где V — собственный вектор, I — единичная матрица. Поэтому результаты вычислений в новой версии, как правило, отличаются. Ниже дан список средств решения векового уравнения, реализованных в системе MATLAB. Несимметрические матрицы могут быть плохо обусловлены при вычислении их собственных значений. Малые изменения элементов матрицы, такие как ошибки округления, могут вызвать большие изменения в собственных значениях. Масштабирование — это попытка перевести каждую плохую обусловленность собственных векторов матрицы в диагональное масштабирование. Однако масштабирование обычно не может преобразовать несимметрическую матрицу в симметрическую, а только пытается сделать (векторную) норму каждой строки равной норме соответствующего столбца. Масштабирование значительно повышает стабильность собственных значений.
» А=[1 1000 10000:0.0001 1 1000:0.000001 0.0001 1] А = 1.0е+004 * 0.0001 0.1000 1.0000 0.0000 0.0001 0.1000 0.0000 0.0000 0.0001 » [F.G]=balance(A) F = 1.0е+004 * 3.2768 0 0
0 0.0032 0 0 0 0.0000 G = 1.0000 0.9766 0.0095 0.1024 1.0000 0.9766 1.0486 0.1024 1.0000 Величина, связывающая погрешность вычисления собственных значений с погрешностью исходных данных, называется числом обусловленности (собственных значений) матрицы и вычисляется следующим образом: cond(V) = norm(V)*norm(inv(V)) где [V.D]=eig(A).[B=D\A*D, а норма каждого ряда масштабированной матрицы приближается к норме столбца с тем же номером;]
это просто полные матрицы со сравнительно большим [Но небольшим по сравнению с числом нулей разреженной матрицы. Эталонное число нулей разреженной матрицы данного размера можно вычислить, применив к полной матрице этого же размера функцию sparse. — Примеч. ред.] числом нулей), а также во всех случаях, где помимо собственных значений необходимо получать и собственные вектора разреженной матрицы, вместо нее рекомендовано использовать eigs(A);
Нужно использовать [W,D]=e1g(A'); W=W, чтобы вычислить левые собственные вектора, которые соответствуют уравнению W*A=D*W.
» В = [3 -12 -.6 2*eps:-2 48 -1 -eps;-eps/8 eps/2 -1 10;-.5 -.5 .3 1] В = 3.0000 -12.0000 -0.60000.0000 -2.0000 48.0000-1.0000-0.0000 -0.0000 0.0000 -1.0000 10.0000 -0.5000 -0.5000 0.3000 1.0000 » [G.H]=eig(B) G = -0.2548 0.7420 -0.4842 0.1956 0.9670 0.0193 -0.0388 0.0276 -0.0015 -0.6181 -0.8575 0.9780 -0.0075 -0.2588 -0.1694 -0.0676 H = 48.5287 0 0 0 0 3.1873 0 0 0 0 0.9750 0 0 0 0 -1.6909
» F=[23 12;3 5:6 0] F = 23 12 3 5 6 0 » [k,l,m]=svd(F) k=
0.8863 -0.4630 0.4630 0.8863 Приведение матриц к форме Шура и Хессенберга Ниже приводятся функции, обеспечивающие приведение матриц к специальным формам Шура и Хессенберга:
Конкретные столбцы матрицы V больше не являются собственными векторами, но каждая пара векторов связана с блоком размера 2x2 в матрице D. Пример: » А-[2 3 6;-4 0 3:1 5 -2] А = 2 3 6 -4 0 3 1 5 -2
Функция qz обеспечивает приведение пары матриц к обобщенной форме Шура:
Обобщенные собственные значения могут быть найдены из следующего условия: A*V*diag(BB) = B*V*diag(AA) Пример: » А=[1 2 3:6 3 0;4 7 0];В=[1 1 1:0 7 4:9 4 1]; » [aa.bb.f,g.h]=qz(A.B) аа = -2.9395 0.4775 0.8751 0 9.5462 3.5985 0 0 3.2073 bb=
-0.7023 -0.7050 -0.0989 0.6867 -0.6343 -0.3552 -0.1877 0.3174 -0.9295 h = -1.0000 -0.4874 -0.0561 0.9778 -1.0000 0.6238 -0.2673 0.4340 -1.0000 Функция qz(A,B,' real') при заданных матрицах А и В возвращает действительные треугольную матрицу ВВ и квазитреугольную матрицу АА с 2x2 диагональными блоками, соответствующими парам сопряженных комплексных значений. Так как матрица АА квазитреугольная, то необходимо решить проблемы обобщения 2x2 для получения подлинных собственных значений. Пример: » А=[1 2 3:6 3 0:4 7 0];В=[1 1 1:0 7 4;9 4 1]; » [aa.bb.f,g.h]=qz(A.B,'real') аа = -2.9395 0.4775 0.8751
0 9.5462 3.5985 0 0 3.2073
bb =
5.5356 3.5345 -2.2935 0 8.4826 6.7128 0 0 0.7667 f = -0.0367 0.7327 -0.6796 -0.1052 -0.6791 -0.7265 -0.9938 0.0448 0.1020 g= -0.7023 -0.7050 -0.0989 0.6867 -0.6343 -0.3552 -0.1877 0.3174 -0.9295 h = -1.0000 -0.4874 -0.0561 0.9778 -1.0000 0.6238 -0.2673 0.4340 -1.0000
Элементы матрицы Хессенберга, расположенные ниже первой поддиагонали, равны нулю. Если матрица симметричная или эрмитова, то матрица Хессенберга вырождается в трехдиагональную. Эта матрица имеет те же собственные значения, что и оригинал, но для их вычисления необходимо меньшее количество операций. Пример: » f=magic(4) f = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 » hess(f) ans= 16.0000 -8.0577 8.8958 6.1595 -11.0454 24.2131 -8.1984 2.1241 0 -13.5058 -4.3894 -7.8918 0 0 3.2744 -1.8237 В этом уроке мы научились:
|