Учебное пособие для студентов высших учебных заведений
 Скачать 5.41 Mb.
|
|
1.4. Функции прикладной численной математики 56 » disp(chol(A)) 1. 0000 2. 0000 3. 0000 0 3. 3166 0. 6030 0 0 19. 7645 Функция lu(A) осуществляет LU-разложение матрицы А в виде произведе- ния нижней треугольной матрицы L (возможно, с перестановками) и верхней тре- угольной матрицы U, так что A = L * U. Обращение к этой функции вида [ L, U, P ] = lu(A) позволяет получить три составляющие этого разложения - нижнюю треугольную матрицу L, верхнюю треугольную U и матрицу перестановок P такие, что P * A =L * U. Приведем пример: A = 1 2 3 2 15 8 3 8 400 » disp(lu(A)) 3.0000 8.0000 400.0000 -0. 6667 9. 6667 -258. 6667 -0. 3333 0. 0690 -148. 1724 » [ L, U, P] = lu(A); » L L = 1.0000 0 0 0. 6667 1. 0000 0 0. 3333 -0. 0690 1. 0000 » U U = 3.0000 8.0000 400.0000 0 9. 6667 -258. 6667 0 0 -148. 1724 » P P = 0 0 1 0 1 0 1 0 0 Из него вытекает, что в первом, упрощенном варианте обращения функция выдает комбинацию из матриц L и U. Обращение матрицы осуществляется с помощью функции inv(А): » disp(inv(A)) 1. 3814 -0. 1806 -0. 0067 -0. 1806 0. 0910 -0. 0005 -0. 0067 -0. 0005 0. 0026 Процедура pinv(А) находит матрицу, псевдообратную матрице А, которая имеет размеры матрицы А’ и удовлетворяет условиям A * P * A = A; P * A * P = P. Например: A = 1 2 3 4 5 5 -1 4 6 0 1.4. Функции прикладной численной математики 57 » P = pinv(A) P = -0.0423 0.0852 0.0704 -0.0480 0.0282 0.0372 0. 0282 0. 0628 0. 1408 -0. 0704 » A*P*A, % проверка 1 ans = 1. 0000 2. 0000 3. 0000 4. 0000 5. 0000 5. 0000 -1. 0000 4. 0000 6. 0000 0. 0000 » P*A*P % проверка 2 ans = -0.0423 0.0852 0.0704 -0.0480 0. 0282 0. 0372 0. 0282 0. 0628 0.1408 -0. 0704 Для квадратных матриц эта операция равнозначна обычному обращению. Процедура [ Q, R, P ] = qr(A) осуществляет разложение матрицы А на три - унитарную матрицу Q, верхнюю треугольную R с диагональными элементами, уменьшающимися по модулю, и матрицу перестановок P - такие что A * P = Q * R. Например: A = 1 2 3 4 5 5 -1 4 6 0 » [Q,R,P] = qr(A) Q = -0.5547 -0.8321 -0.8321 0.5547 R = -7.2111 -2.7735 -4.9923 -4.7150 -0.2774 0 -4.1603 -0.2774 1.9415 -2.2188 P = 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 Определение характеристического полинома матрицы A можно осуще- ствить с помощью функции poly(A). Обращение к ней вида p=poly(A) дает воз- можность найти вектор-строку p коэффициентов характеристического полинома p(s) = det(s*E - A) = p 1 *s n + ... + p n *s + p n+1 , где Е - обозначение единичной матрицы размером (n*n). Например : » A = [1 2 3; 5 6 0; -1 2 3] A = 1 2 3 5 6 0 -1 2 3 » p = poly(A) p = 1. 0000 -10. 0000 20. 0000 -36. 0000 Вычисление собственных значений и собственных векторов матрицы осуществляет процедура eig(А). Обычное обращение к ней позволяет получить вектор собственных значений матрицы А, т. е. корней характеристического полинома матрицы. Если же обращение имеет вид: 1.4. Функции прикладной численной математики 58 [ R, D ] = eig(A), то в результате получают диагональную матрицу D собственных значений и матрицу R правых собственных векторов, которые удовлетворяют условию A*R = R * D. Эти векторы являются нормированными так, что норма любого из них рав- на единице. Приведем пример: A = 1 2 3 -1 8 16 -5 100 3 » disp(eig(A)) 1.2234 45. 2658 -34. 4893 » [ R,D] = eig(A) R = 0.9979 -0.0798 -0.0590 0.0492 -0.3915 -0.3530 0.0416 -0.9167 0.9338 D = 1. 2234 0 0 0 45. 2658 0 0 0-34. 4893 Сингулярное разложение матрицы осуществляет процедура svd(А). Уп- рощенное обращение к ней позволяет получить сингулярные числа матрицы А. Более сложное обращение вида: [ U, S, V ] = svd(A) позволяет получить три матрицы - U, которая состоит из ортонормированных собственных векторов, отвечающих наибольшим собственным значениям матри- цы А*А T ; V - из ортонормированных собственных векторов матрицы А T *А и S - диагональную матрицу, которая содержит неотрицательные значения квадрат- ных корней из собственных значений матрицы А T *А (их называют сингулярными числами). Эти матрицы удовлетворяют соотношению: A = U * S * VT. Рассмотрим пример: » disp(svd(A)) 100.5617 15. 9665 1. 1896 » [U,S,V] = svd(A) U = -0.0207 0.1806 -0.9833 -0.0869 0.9795 0.1817 -0.9960 -0.0892 0.0045 S = 100.5617 0 0 0 15.9665 0 0 0 1.1896 V = 0. 0502 -0. 0221 -0. 9985 -0. 9978 -0. 0453 -0. 0491 -0. 0442 0. 9987 -0. 0243 1.4. Функции прикладной численной математики 59 Приведение матрицы к форме Гессенберга осуществляется процедурой hess(А). Например: A = 1 2 3 -1 8 16 -5 100 3 » disp(hess(A)) 1. 0000 -3. 3340 -1. 3728 5. 0990 25. 5000 96. 5000 0 12. 5000 -14. 5000 Более развернутое обращение [P,H] = hess(A) дает возможность получить, кроме матрицы Н в верхней форме Гессенберга, также унитарную матрицу преоб- разований Р, которая удовлетворяет условиям: A = P * H * P; P' * P = eye(size(A)). Пример: » [P,H] = hess(A) P = 1.0000 0 0 0 -0.1961 -0.9806 0 -0.9806 0.1961 H = 1.0000 -3.3340 -1.3728 5. 0990 25. 5000 96. 5000 0 12. 5000 -14. 5000 Процедура schur (А) предназначена для приведения матрицы к форме Шура. Упрощенное обращение к ней приводит к получению матрицы в форме Шура. Комплексная форма Шура - это верхняя треугольная матрица с собст- венными значениями на диагонали. Действительная форма Шура сохраняет на диагонали только действительные собственные значения, а комплексные изо- бражаются в виде блоков (2*2), частично занимая нижнюю поддиагональ. Обращение [U,T] = schur(A) позволяет, кроме матрицы Т Шура, получить также унитарную матрицу U, удовлетворяющую условиям: A = U * H * U'; U' * U = eye(size(A)). Если начальная матрица А является действительной, то результатом будет действительная форма Шура, если же комплексной, то результат выдается в виде комплексной формы Шура. Приведем пример: » disp(schur(A)) 1.2234 -6.0905 -4.4758 0 45. 2658 84. 0944 0 0. 0000 -34. 4893 » [U,T] = hess(A) U = 1.0000 0 0 0 -0.1961 -0.9806 0 -0.9806 0.1961 T = 1. 0000 -3. 3340 -1. 3728 5. 0990 25. 5000 96. 5000 0 12. 5000 -14. 5000 1.4. Функции прикладной численной математики 60 Функция [U,T] = rsf2csf(U,T) преобразует действительную квазитреуголь- ную форму Шура в комплексную треугольную: » [U,T] = rsf2csf(U,T) U = -0.9934 -0.1147 0 -0.0449 0.3892 -0.9201 -0.1055 0.9140 0.3917 T = 1. 4091 -8. 6427 10. 2938 0 45. 1689 -83. 3695 0 0 -34. 5780 Процедура [AA, BB, Q, Z, V] = qz(A,B) приводит пару матриц А и В к обобщенной форме Шура. При этом АА и ВВ являются комплексными верхними треугольными матрицами, Q, Z - матрицами приведения, а V - вектором обобщен- ных собственных векторов такими, что Q * A * Z = AA; Q * B * Z = BB. Обобщенные собственные значения могут быть найдены, исходя из такого условия: A * V * diag(BB) = B * V * diag(AA). Необходимость в одновременном приведении пара матриц к форме Шура возникает во многих задачах линейной алгебры - решении матричных уравнений Сильвестра и Риккати, смешанных систем дифференциальных и линейных алгеб- раических уравнений. Пример. Пусть задана система обычных дифференциальных уравнений в неявной форме Коши с одним входом и одним выходом такого вида: u y Q x R x b c ⋅ + ⋅ = ⋅ = ⋅ + ⋅ & ; u y x d u причем матрицы Q, R и векторы b, c и d равны соответственно Q = 1. 0000 0 0.1920 1. 0000 R = 1. 1190 -1. 0000 36.4800 1. 5380 b = 31. 0960 0. 1284 c = 0.6299 0 d = -0. 0723 Необходимо вычислить значения полюсов и нулей соответствующей пере- даточной функции. Эта задача сводится к отысканию собственных значений λ, которые удовле- творяют матричным уравнениям: 1.4. Функции прикладной численной математики 61 R r Q r R b c d r Q r ⋅ = − ⋅ ⋅ − ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⋅ = ⋅ ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ⋅ λ λ ; 0 0 0 Решение первого уравнения позволяет вычислить полюсы передаточной функции , а второго - нули. Ниже приведена совокупность операторов, которая приводит к расчету по- люсов: » [AA, BB] = qz(R,-Q) % Приведение матриц к форме Шура AA = 5.5039 + 2.7975i 24.8121 -25.3646i 0.0000 - 0.0000i 5.5158 - 2.8036i BB = -0. 6457 + 0. 7622i -0. 1337 + 0. 1378i 0 -0. 6471 - 0. 7638i » diag(AA) . /diag(BB) % Расчет полюсов ans = -1. 4245 - 6. 0143i -1. 4245 + 6. 0143i Расчет нулей осуществляется таким образом: » A = [-R b % Формирование c d] % первой матрицы A = -1.1190 1.0000 0.1284 -36. 4800 -1. 5380 31. 0960 0.6299 0 -0. 0723 » B = [ -Q zeros(size(b)) % Формирование zeros(size(c)) 0 ] % второй матрицы B = -1.0000 0 0 -0. 1920 -1. 0000 0 0 0 0 » [AA,BB] = qz(A,B) % Приведение матриц к форме Шура AA = 31.0963 -0.7169 -36.5109 0.0000 1.0647 0.9229 0 0.0000 0.5119 BB = 0 0.9860 -0.2574 0 0. 0657 0. 9964 0 0 -0. 0354 » diag(AA) . /diag(BB) % Вычисление нулей ans = Inf 16. 2009 -14. 4706 Вычисление собственных значений матричного полинома осуществляет процедура polyeig. Обращение [ R, d ] = polyeig(A0, A1,... , Ap ) позволяет решить полную проблему собственных значений для матричного поли- нома степени р вида (A0 + λ *A1 + ... + p λ *Ap)*r = 0. 1.4. Функции прикладной численной математики 62 Входными переменными этой процедуры являются р+1 квадратные матри- цы A0, A1, ... Ap порядка n. Исходными переменными - матрица собственных векторов R размером (n*(n*p)) и вектор d собственных значений размером (n*p). Функция polyvalm предназначена для вычисления матричного полинома вида Y X p X p X p X p X p n n n n ( ) = ⋅ + ⋅ + + ⋅ + ⋅ + − − 1 1 2 2 1 0 по заданному значению матрицы Х и вектора p = [pn, pn-1, ... , p0] коэффициентов полинома. Для этого достаточно обратиться к этой процедуре по схеме: Y = polyvalm(p, X). Пример: p = 1 8 31 80 94 20 » X X = 1 2 3 0 -1 3 2 2 -1 » disp(polyvalm(p,X)) 2196 2214 2880 882 864 1116 1332 1332 1746 Примечание. Следует различать процедуры polyval и polyvalm. Первая вы- числяет значение полинома для любого из элементов матрицы аргумента, а вто- рая при вычислении полинома возводит в соответствующую степень всю мат- рицу аргумента. Процедура subspace(А,В) вычисляет угол между двумя подпространствами, которые "натянуты на столбцы" матриц А и В. Если аргументами являются не матрицы, а векторы A и B, вычисляется угол между этими векторами. 1.4.4. Аппроксимация и интерполяция данных Полиномиальная аппроксимация данных измерений, которые сформиро- ваны как некоторый вектор Y, при некоторых значениях аргумента, которые обра- зуют вектор Х такой же длины, что и вектор Y, осуществляется процедурой polyfit (X, Y, n). Здесь n - порядок аппроксимирующего полинома. Результатом действия этой процедуры является вектор длиной (n +1) из коэффициентов ап- проксимирующего полинома. Пусть массив значений аргумента имеет вид: x = [1 2 3 4 5 6 7 8], а массив соответствующих значений измеренной величины - вид: y = [ -1.1 0.2 0.5 0.8 0.7 0.6 0.4 0.1]. Тогда, применяя указанную функцию при разных значениях порядка ап- проксимирующего полинома, получим: » x = [1 2 3 4 5 6 7 8]; » y = [ -1.1 0.2 0.5 0.8 0.7 0.6 0.4 0.1]; 1.4. Функции прикладной численной математики 63 » polyfit(x,y,1) ans = 0. 1143 -0. 2393 » polyfit(x,y,2) ans = -0. 1024 1. 0357 -1. 7750 » polyfit(x,y,3) ans = 0. 0177 -0. 3410 1. 9461 -2. 6500 » polyfit(x,y,4) ans = -0. 0044 0. 0961 -0. 8146 3. 0326 -3. 3893. Это означает, що заданную зависимость можно аппроксимировать или пря- мой y x x ( ) , ,2393 = − 0 1143 0 , или квадратной параболой y x x x ( ) , , , = − + − 0 1024 1 0357 1 775 2 , или кубической параболой y x x x x ( ) , , , , = − + − 0 0177 0 341 1 9461 2 65 3 2 , или параболой четвертой степени y x x x x x ( ) , , , , , = − + − + − 0 0044 0 0961 0 8146 3 0326 3 3893 4 3 2 Построим в одном графическом поле графики заданной дискретной функ- ции и графики всех полученных при аппроксимации полиномов: x = [1 2 3 4 5 6 7 8]; y = [ -1.1 0.2 0.5 0.8 0.7 0.6 0.4 0.1]; P1=polyfit(x,y,1) ; P2=polyfit(x,y,2); P3=polyfit(x,y,3); P4=polyfit(x,y,4) ; stem(x,y); x1 = 0.5 : 0.05 : 8.5; y1=polyval(P1,x1); y2=polyval(P2,x1); y3=polyval(P3,x1); y4=polyval(P4,x1); hold on plot(x1,y1,x1,y2,x1,y3,x1,y4), grid, set(gca, 'FontName', 'Arial Cyr', 'FontSize', 16), title('Полиномиальная аппроксимация '); xlabel('Аргумент'); ylabel('Функция') Результат представлен на рис. 1.18. Функция spline(X,Y,Xi) осуществляет интерполяцию кубическими сплай- нами . При обращении Yi = spline(X,Y,Xi) 1.4. Функции прикладной численной математики 64 Рис. 1.18 она интерполирует значение вектора Y, заданного при значениях аргумента, пред- ставленных в векторе Х, и выдает значение интерполирующей функции в виде вектора Yi при значениях аргумента, заданных вектором Xi. В случае, если вектор Х не указан, по умолчанию принимается, что он имеет длину вектора Y и любой его элемент равен номеру этого элемента. В качестве примера рассмотрим интерполяцию вектора x = -0.5:0.1:0.2; y = [-1.1 0.2 0.5 0.8 0.7 0.6 0.4 0.1]; x1 = -0.5:0.01:0.2; y2 = spline(x,y,x1); plot (x,y,x1,y2), grid set(gca,'FontName','Arial Cyr','FontSize',16), title('Интерполяция процедурой SPLINE '); xlabel('Аргумент'); ylabel('Функция') Результат приведен на рис. 1.19. Одномерную табличную интерполяцию осуществляет процедура interp1. Обращение к ней в общем случае имеет вид: Yi = interp1(X,Y,Xi,’<метод>’), и позволяет дополнительно указать метод интерполяции в четвертом входном ар- гументе: ' nearest' - ступенчатая интерполяция; 'linear' - линейная; ‘cubic' - кубическая; ‘spline' - кубическими сплайнами. 1.4. Функции прикладной численной математики 65 Рис. 1.19 Если метод не указан, осуществляется по умолчанию линейная интерполя- ция. Например, (для одного и того же вектора): x = -0.5:0.1:0.2; y = [-1.1 0.2 0.5 0.8 0.7 0.6 0.4 0.1]; x1 = -0.5:0.01:0.2; y1 = interp1(x,y,x1); y4 = interp1(x,y,x1,'nearest'); y2 = interp1(x,y,x1,'cubic'); y3 = interp1(x,y,x1,'spline'); plot (x1,y1,x1,y2,x1,y3,x1,y4), grid plot (x1,y1,x1,y2,'.',x1,y3,'--',x1,y4,':'), grid set(gca,'FontName','Arial Cyr','FontSize',8), legend('линейная','кубическая','сплайновая','ступенчатая',0) set(gca,'FontName','Arial Cyr','FontSize',14), set(gcf,'color','white') title('Интерполяция процедурой INTERP1 '); xlabel('Аргумент'); ylabel('Функция') Результаты приведены на рис.1.20. |