Матлаб. К. Ю. Петрова введение в matlab учебное пособие
 Скачать 2.57 Mb.
|
|
Матричная экспонента В MATLAB имеется несколько команд для вычисления функций от квадратных матриц – это polyvalm, sqrtm, expm, logm. Команда polyvalm предназначена для вычисления значений матричного полинома E a A a A a A F n n 0 1 ) ( . Ее первый аргумент – вектор коэффициентов полинома, второй – матрица А. Функция sqrtm (A) вычисляет квадратный корень из матрицы A. Если матрица A имеет отрицательные собственные числа, результат будет комплексной матрицей. 65 Команда expm (A) обеспечивает вычисление матричной экспоненты A e путем использования аппроксимации Паде, а операция logm (A) – вычисление матричного логарифма A ln . Напомним определение матричной экспоненты. Как известно, обычная экспонента a e задается разложением в ряд Тэйлора ! ! 3 ! 2 1 3 2 n a a a a e n a Матричная экспонента для заданной квадратной матрицы А вводится как сумма аналогичного ряда ! ! 2 2 n A A A E e n A (*) Она представляет собой квадратную матрицу того же размера, что и матрица А, и обладает рядом свойств, характерных для обычной экспоненты. Особенно легко находится матричная экспонента для диагональных матриц – для этого достаточно взять обычные экспоненты от диагональных элементов (заметим, что это единственный случай, когда expm (А)=exp(А)). В общем случае вычисление матричной экспоненты – трудоемкая процедура. Для контроля правильности результата можно использовать проверку равенства , det trA A e e связывающего определитель матричной экспоненты со следом исходной матрицы. Пример. Найдем матричные экспоненты для следующих трех матриц: 0 1 1 0 , 1 0 2 1 , 2 0 0 1 3 2 1 A A A Возводя их в степени и выполняя поэлементное суммирование членов ряда (*), получаем: 1 cos 1 sin 1 sin 1 cos , 0 2 , 0 0 3 2 1 2 A A A e e e e e e e e След матрицы А 1 равен трем, определитель матрицы 1 A e равен е 3 , так что равенство trA A e e det выполняется. Легко проверить его справедливость и для остальных матриц. Вычислим те же матричные экспоненты в MATLAB с помощью команды expm: >>A1=[1 0;0 2]; >>e1=expm(A1) >>A2=[1 2;0 1]; >>e2=expm(A2) >>A3=[0 1;-1 0]; >>e3=expm(A3) 2.7183 0 0 7.3891 2.7183 5.4366 0 2.7183 0.5403 0.8415 -0.8415 0.5403 Элементами полученных матриц служат величины , 4366 , 5 2 , 7,3891 , 7183 , 2 2 e e e а также sin 57,3 0 и cos 57,3 0 У матрицы А 3 собственные числа были чисто мнимые i. В общем случае матрица второго порядка имеет комплексные собственные числа i , преобразованием подобия она может быть приведена к виду A . Матричная экспонента для таких матриц определяется формулой cos sin sin cos a A e e 66 Пример. Для матрицы 1 5 0 5 0 1 A получаем e e A / 5 0 cos 5 0 sin 5 0 sin 5 0 cos Применение функции expm дает тот же результат: >>exp(-1)*[cos(-1/2), sin(-1/2); -sin(-1/2),cos(-1/2)] >> A=[-1 -.5;.5 -1]; >>expm(A) >>A=[-1 -.5;.5 -1]; >>expm(sym(A)) ans = 0.3228 -0.1764 0.1764 0.3228 ans = 0.3228 -0.1764 0.1764 0.3228 ans = [exp(-1)*cos(1/2),-exp(-1)*sin(1/2)] [exp(-1)*sin(1/2),exp(-1)*cos(1/2)] Последний столбец и следующий пример демонстрируют возможность тулбокса SYMBOLIC вычислять матричную экспоненту в символьном виде. Пример. Найти матричную экспоненту e tA , если 5 3 3 2 2 1 2 1 2 A Решение в MATLAB. >> A=[-2 1 -2;1 -2 2;3 -3 5]; eig(a)' ans = -1 3 -1 >> syms t, F=expm(t*A); 4*F 4F = [ -exp(3*t)+5*exp(-t), -exp(-t)+exp(3*t), 2*exp(-t)-2*exp(3*t)] [ -exp(-t)+exp(3*t), -exp(3*t)+5*exp(-t), -2*exp(-t)+2*exp(3*t) ] [ -3*exp(-t)+3*exp(3*t), 3*exp(-t)-3*exp(3*t), 6*exp(3*t)-2*exp(-t)] Это означает, что матричная экспонента имеет вид: t t t t t t t t t t t t t t t t t t e e e e e e e e e e e e e e e e e e F 2 6 3 3 3 3 2 2 5 2 2 5 4 1 3 3 3 3 3 3 3 3 3 Структуру собственных чисел и векторы матрицы A можно найти с помощью команд eig и jordan: >> [V,D]=eig(sym(A)), [v,d]=jordan(sym(A)) V D v d [ -1, 1, -2] [ 1, 1, 0] [ 3, 0, 1] [ 3, 0, 0] [ 0, -1, 0] [ 0, 0, -1] [ -3/4, -1/4, -2] [ -1/4, 1/4, 0] [ 1/4, 3/4, 1] [ -1, 0, 0] [ 0, 3, 0] [ 0, 0, -1] Хотя у матрицы А есть кратное собственное число -1, она диагонализируема, т.е. является матрицей простой структуры. Одно из важных применений матричной экспоненты – аналитический расчет свободного движения линейных динамических объектов. Математически эта задача сводится к решению системы однородных дифференциальных уравнений , , ) 0 ( , 0 CX y X X AX X где А – квадратная матрица, С – прямоугольная матрица, Х – вектор переменных состояния. 67 Стандартное средство для решения этой задачи в MATLAB – команда initial тулбокса CONTROL. Другой способ – использование матричной экспоненты. Он опирается на запись решения в форме ) ( 0 X Ce t X At Такой способ может оказаться более удобным, если надо получить несколько решений одной и той же системы при различных начальных условиях. При этом матричная экспонента вычисляется только один раз, в то время как при использовании команды initial придется многократно повторять процесс моделирования. Пример. Колебания маятника описываются дифференциальным уравнением 0 4 2 , 0 y y y Требуется получить графики y(t) при четырех вариантах начальных условий: 0 ; 1 ; 75 , 0 ; 5 , 0 ; 25 , 0 0 0 y y Решение. Преобразуем уравнение к форме Коши, обозначая : ; 2 1 y x y x 2 , 0 4 , 2 1 2 2 1 x x x x x Выпишем матрицы А, С этой системы и матрицу Х 0 четырех вариантов начальных условий: 0 0 0 0 1 75 , 0 5 , 0 25 , 0 , 0 1 , 2 , 0 4 1 0 0 X C A Составим программу expmpend.m для получения решения согласно формуле : 0 X Ce y At % Program expmpend A=[0 1; -4 -0.2]; c=[1 0]; % матрицы А, С X0=[1/4 1/2 3/4 1; 0 0 0 0]; % матрица начальных условий h=0.1; t=0:h:10; N=length(t); % массив времени Eh=expm(A*h); F=zeros(N,2); F(1,:)=c; for i=2:N F(i,:)=F(i-1,:)*Eh; % цикл вычисления массива At Ce F размера Nх2 end Y=F*X0; % решение для четырех начальных условий plot(t,Y), grid % вывод графиков Графики четырех полученных решений показаны на рис. 4.21. 68 0 2 4 6 8 10 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 T y 1 , y 2 , y 3 , y 4 Y Рис. 4.21 Более короткий путь их получения, не требующий организации цикла, основан на использовании символьных вычислений и вывода графиков командой ezplot. A=[0 1; -4 -0.2]; c=[1 0]; % матрицы А, С X0=[1/4 1/2 3/4 1; 0 0 0 0]; % матрица начальных условий syms T E=expm(A*T);Y=c*E*X0; ezplot(Y(1),[0 10]),hold on,ezplot(Y(2),[0 10]),grid,ezplot(Y(3),[0 10]), ezplot(Y(4),[0 10]),hold off, title('y_1, y_2, y_3, y_4') К сожалению, ezplot не позволяет выводить в одно графическое окно сразу несколько кривых, поэтому приходится дополнительно прибегать к команде hold. Канонические формы матриц Преобразование подобия При решении многих задач линейной алгебры подходящей заменой переменных удается упростить вид матриц, входящих в эти задачи. Одно из наиболее распространенных преобразований описывается формулой , 1 TAT B где А – исходная квадратная матрица, Т – матрица невырожденной замены переменных, В – результат преобразования. Матрицы A и В, связанные таким преобразованием, называются подобными, а само преобразование – подобием. Преобразование подобия часто возникает при решении прикладных задач. Наиболее распространенный пример – решение системы линейных уравнений A X = Y. При замене переменных матрица A этой системы подвергается преобразованию подобия, изменяющему как ее внешний вид, так и многие свойства (например, обусловленность, что необходимо учитывать при решении линейных систем уравнений на ЭВМ). Другой важный пример – моделирование свободного движения линейной системы AX X Переходя к новым переменным X = T Y, получаем систему , 1 ATY T Y матрица которой подобна матрице исходной системы. 69 Одна из команд MATLAB, использующая преобразование подобия, называется balance. Она подбирает диагональную матрицу преобразования Т таким образом, чтобы выровнять нормы строк и столбцов матрицы А по величине. По существу, речь идет о разумном масштабировании переменных, улучшающем обусловленность матрицы (иногда говорят о балансировке или уравновешивании элементов). Пример. Сформируем матрицу А с большим разбросом элементов и выполним ее балансировку: >>A=[1 100 0.1; 1 0 0; 0 100 1], [T, B]=balance(A) В результате получаем: A T B 1 100 0.1 1 0 0 0 100 1 2 0 0 0 0.25 0 0 0 32 1 12.5 1.6 8 0 0 0 0.7813 1 Разброс элементов в матрице В по сравнению с исходной матрицей А заметно уменьшился. На диагонали матрицы Т стоят масштабные коэффициенты, они представляют собой степени двойки (чтобы избежать погрешностей округления). За счет преобразования подобия можно достичь значительного упрощения вида матрицы А, придав ей регулярную структуру и сделав нулевыми или единичными большинство ее элементов. Стандартные формы, к которым может быть приведена матрица А путем преобразования подобия, называются каноническими. Наибольшую известность получили сопровождающая каноническая форма (companion form), диагональная каноническая форма (форма Жордана), треугольная форма Шура, форма Хессенберга, форма Шварца. В пакете MATLAB существуют команды для получения этих форм: compan, jordan, schur, hess. Остановимся на двух первых из них. Фробениусова каноническая форма Параметрами этой канонической формы служат коэффициенты характеристического полинома 0 1 1 1 a p a p a p A pE n n n Если степень минимального полинома матрицы A равна n, то ее фробениусова каноническая форма имеет вид 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 1 2 2 1 0 a E a a a a a B n n где a = (a 0 , ..., a n-1 ). Все варьируемые элементы матрицы расположены в последней строке и совпадают с точностью до знака с коэффициентами ее характеристического полинома. По этой причине ее называют также матрицей, сопровождающей характеристический полином, ее собственные числа равны его корням. Вместо матрицы В можно использовать подобные ей матрицы 70 , 0 , 0 , 0 a E E a E a получающиеся из нее транспонированием и перенумерацией переменных. Таким образом, имеем два строчных и два столбцовых варианта фробениусовых матриц. Если степени минимального и характеристического полинома матрицы не совпадают, то фробениусова каноническая форма становится клеточно-диагональной с клетками вида В. Процедура получения фробениусовой канонической формы довольно проста. Ранее отмечалось, что команда Р=poly(A) дает вектор коэффициентов характеристического полинома матрицы А. Команда compan решает обратную задачу – она строит матрицу по характеристическому полиному. Набрав B=compan( Р), получим матрицу В, первая строка которой образована коэффициентами полинома Р (деленными на старший коэффициент и взятыми с минусом), а все остальные элементы равны нулю, за исключением первой поддиагонали, на которой стоят единицы. Пример. Пусть задан полином Р третьего порядка: 4 3 2 P 2 3 x x x Для построения сопровождающей матрицы набираем команды: Р=[1 2 3 4]; B=compan(Р). Результатом будет матрица 0 1 0 0 0 1 4 3 2 B Находя ее характеристический полином (это можно сделать с помощью команды poly(B) ), убеждаемся, что он совпадает с Р. В соответствии со сказанным выше, сопровождающими полином Р будут также матрицы 2 1 0 3 0 1 4 0 0 , 2 3 4 1 0 0 0 1 0 , 0 0 4 1 0 3 0 1 2 Все они подобны и имеют одинаковые характеристические полиномы. Если дана произвольная квадратная матрица А с разными собственными числами, то для преобразования ее в сопровождающую каноническую форму достаточно двух команд: p=poly(A); B=compan(p). Матрицы А и В будут иметь одинаковые характеристические полиномы и будут подобны. В случае алгебраически кратных собственных чисел характеристические полиномы матриц А и В будут по-прежнему совпадать, однако подобие уже не будет иметь место. Простого способа получения сопровождающей канонической формы в MATLAB для этого случая не существует. Упражнение 1. Выполните описанный выше переход для матрицы А=eye(2) и проверьте, будут ли подобными матрицы А и В. Упражнение 2. Пусть матрицы А и В =compan(poly(A)) подобны. Подумайте, как найти связывающую их матрицу преобразования Т? Проверьте ваш способ в MATLAB, используя формулу 1 AT T B 71 Жорданова каноническая форма Каноническая форма Жордана и разные ее модификации широко используются в приложениях, возможность ее нахождения предусмотрена в большинстве математических пакетов. Параметрами канонической формы Жордана служат корни характеристического уравнения матрицы, поэтому в общем случае это – комплексная матрица. Если все корни различны, то жорданова каноническая форма имеет диагональный вид, при этом диагональными элементами являются собственные числа исходной матрицы. Жорданова форма полезна при решении многих задач. Один из характерных примеров – решение системы дифференциальных уравнений AX X Если все собственные числа n d d , , 1 матрицы А различны, то существует такая невырожденная матрица V, что заменой переменных VY X исходная система приводится к «диагональному» виду, который описывается уравнениями , , , 2 2 2 1 1 1 n n n y d y y d y y d y Решая их, получаем , , , , 2 1 2 2 1 1 t d n n t d t d n e c y e c y e c y где постоянные i c зависят от начальных условий. Матричное описание полученной системы дифференциальных уравнений имеет вид , DY Y 0 0 0 0 0 0 0 0 0 2 1 n d d d D Матрицы A и D связаны соотношением , 1 AV V D где V – матрица, составленная из собственных векторов матрицы А. Матрица D называется диагональной или жордановой канонической формой матрицы А для матриц простой структуры. В MATLAB матрицы V и D можно найти с помощью команды [V, D]=eig(A). |