Вычислительная математика лекции. Конспект лекций. Санкт петербург 2011 2 Оглавление
 Скачать 3.86 Mb.
|
|
9.6.2. Определение матричных функций с помощью рядов. Рассмотрим функции, представимые степенными рядами 2 3 3 5 2 4 1 2! 3! sin( ) 3! 5! cos( ) 1 2! 4! x x x e x x x x x x x x 89 Функцией f(A) от матрицы А будем называть тот предел, который получается при подстановке матрицы А в соответствующее степенное разложение. 2 3 2! 3! A A A e E A При таком способе определения матричных функций сохраняется свойство, ранее доказанное для степенных многочленов : λ(f(А))=f(λ(A)). 9.7. Многочлен Лагранжа – Сильвестра. Оказывается матричную функцию можно определить, используя не бесконечный степенной ряд, но алгебраический полином конечной степени. Как и ранее задача ставится следующим образом: задана функция f(λ), требуется вычислить f(A). Многочлен Лагранжа – Сильвестра это интерполяционный многочлен, построенный на спектре матрицы А для заданной функции f(λ). Пусть минимальный аннулирующий многочлен имеет вид ψ(λ) = (λ-λ 1 ) m 1 (λ-λ 2 ) m 2 …(λ-λ s ) m s , m 1 + m 2 +…+ m s =m. Набор из m чисел f(λ k ), f’(λ k ), …, f (m k -1) (λ k ), k=1,2,…,s будем называть значениями f(λ) на спектре матрицы А. Среди всех многочленов с комплексными коэффициентами, принимающих те же значения на спектре, что и f(λ) имеется один и только один многочлен r(λ) степени не выше m-1. Этот многочлен однозначно определяется интерполяционными условиями r(λ k )= f(λ k ), r’(λ k )= f’(λ k ),…, r (m k -1) (λ k )= f (m k -1) (λ k ), (k=1,2,…,s). Иными словами это интерполяционный полином, построенный на указанных узлах. 90 З а м е ч а н и е. В общем случае набор интерполяционных условий содержит кратные узлы ( в одном и том же узле заданы значения интерполируемой функции и ряда её производных). Существует обобщение интерполяционной формулы Лагранжа на этот случай (это общение называют полиномом Эрмита). Многочлен r(λ) называют интерполяционным многочленом Лагранжа – Сильвестра для функции f(λ) на спектре матрицы А. Т е о р е м а Л а г р а н ж а – С и л ь в е с т р а. Пусть f(λ) функция, определенная на спектре матрицы А, а r(λ) – соответствующий интерполяционный многочлен Лагранжа – Сильвестра. Тогда f(A) = r(A). З а м е ч а н и е. Если минимальный многочлен ψ(λ) не имеет кратных корней, то для того чтобы f(A) имело смысл, достаточно, чтобы f(λ) была определена в характеристических точках λ i (i = 1,2,…,m). Если же кратные корни имеются, то в некоторых характеристических точках должны быть определены и производные от f(λ) до известного порядка. Для матриц простой структуры всегда m=n. В первом случае интерполяционный многочлен Лагранжа – Сильвестра строится традиционным образом: ( ) 1 ( ) ( ) ( ) ( ) , r( )=f ( ) ( ) ( ) ( ) '( ) ! m m k k k k f f r r m . Подставляя в формулах матрицу А вместо λ, получим ∆r(А) = 0, так как ω(А)=ψ(А)=0.Таким образом 1 ( ) ( ) ( ) ( ) '( ) m k k k k A f f A A E E Во втором случае способ построения интерполяционного полинома усложняется. В интерполяционных условиях присутствуют кратные узлы. Как уже указывалось, существует обобщение способа построения интерполяционной формулы Лагранжа на этот случай. В 91 соответствии с этим обобщением многочлен Лагранжа – Сильвестра как и прежде можно записать в виде формулы: ( 1) 1 2 , 1 ( ) [ ( ) ( ) '( ) ( ) ... ( ) ( )]. k k s m k k k k k k m k r f f f Функция φ k,j (λ) (j=1,2,…,m k ; k = 1,2 …,s) не зависит от f(λ) и представляет собой интерполяционный многочлен Лагранжа – Сильвестра для функции, у которой все значения на спектре матрицы А равны нулю за исключением одного f (j-1) (λ k ), равного 1. Остаточный член многочлена Лагранжа – Сильвестра и в этом случае равен ( ) ( ) r( )=f ( ) ( ) ( ) ! m f r m . Поэтому в этих условиях попрежнему f(A) = r(A). П р и м е р. Пусть ψ(λ) = (λ-λ 1 ) 2 (λ-λ 2 ) 3 . Тогда (2) 1 11 1 12 2 21 2 22 2 23 ( ) ( ) ( ) '( ) ( ) ( ) ( ) '( ) ( ) ( ) ( ). r f f f f f При расчете матричной функции можно исключить этап формирования многочлена Лагранжа – Сильвестра, что демонстрируется следующим примером. П р и м е р 2 Вычислить f(At) =e At , f(λ)=e λt для матрицы А= 2 1 1 0 1 1 1 1 1 с минимальным многочленом ψ(λ) = (λ-1) 2 (λ-2). Здесь t – скалярный коэффициент. Расчетная формула для данного случая f(At)=f(1)Z 1 +f'(1)Z 2 +f(2)Z 3 , где Z 1 , Z 2 , Z 3 искомые матрицы, не зависящие от f(λ). Для их однозначного определения необходимо сформировать три условия. (1) Пусть f(λ)= λ 0 =1. Тогда Z 1 + Z 3 =E. (2) Пусть f(λ)=λ-1. Тогда Z 2 + Z 3 =A-E 92 (3) Пусть f(λ)= (λ-1)( λ-2) . Тогда - Z 2 =(A-E)(A-2E). Решая систему трех линейных алгебраических уравнений, получаем 2 (2 ) ( )( 2 ) ( ) E A A Z A E A E A E 2 2 1 0 0 1 1 1 0 0 0 ( ) (1) 1 0 0 '(1) 1 1 1 (2) 1 1 0 1 1 1 0 0 0 1 1 0 (1) '(1) '(1) '(1) (1) '(1) (2) '(1) (2) '(1) (1) (2) (1) (2) (1) t t t t t t t t f At f f f f f f f e te te te f f f f f f e te e e te f f f f f 2 2 t t t t t t t te e e e e e 9.8. Свойства матричной экспоненты. 1. At At At de Ae e A dt Доказательство. 2 2 3 3 3 2 2 2 2 3 3 [ ...] 2! 3! 2! [ ...] 2! 3! At d A t A t A t E At A A t dt A t A t A E At Ae 2. 2 2 3 3 2 2 3 0 0 0 2 2 3 2 2 3 3 1 1 [ ...] [ ...] 2! 3! 2! 3! ... [ ] ( ) 2! 3! 2! 3! t t A t At A A A A e d E A d E At A t A t A t Et E At E A e E A 3. Функции от одной матрицы коммутируют. 1 1 ( ) ( ) At At e E A A e E 9.9. Использование канонической формы Жордана при вычислениях матричных функций. Будем вычислять матричную функцию e At 93 2 2 3 3 2! 3! At A t A t e E At Используя подобное преобразование, получим At=SJtS -1 , где J форма Жордана. A k t k = (SJtS -1 ) k = (SJtS -1 SJtS -1 … SJtS -1 ) = S(Jt) k S -1 . Таким образом, 0 1 2 3 1 1 1 2 3 1 1 0 1 0 1 2 3 S Jt S S Jt S S Jt S 2! 3! Jt Jt [ Jt ...] 2! 3! ( , ,..., ); ( ) (( ) , ( ) ,..., ( ) ); Jt Jt Jt ( , ,..., ). 2! 3! s At Jt k k k k s s J t J t J t Jt e E S E S Se S Jt diag J t J t J t Jt diag J t J t J t e E diag e e e Если J i клетка размером k для λ = λ i , то J i = λ i E+J (1) , где J (1) матрица размера "к", у которой элементы первой наддиагонали равны 1, а остальные элементы нулевые. Согласно правилам умножения матриц (J (1) ) p = J (p) , где J (p) матрица размера "к", у которой элементы p-ой наддиагонали равны 1, а остальные элементы нулевые. При p ≥k матрица J (p) становится нулевой. Матричная экспонента (1) (1) ( ) (1) (2) 2 ( 1) 1 1 1 ( ) 2! ( 1)! i i k k i k i J t E J t t J t t k k k k k k k k e e e E e e E E J t J t J t k Таким образом, 94 2 1 2 1 2! ( 1)! 0 1 ( 2)! 0 0 0 0 0 0 1 i i k k J t t t t t k t t k e e t 9.10. Дифференциальные уравнения и матричная экспонента. Рассматриваются обыкновенные линейные дифференциальные уравнения с постоянными коэффициентами. Обозначим x i = x i (t), i= 1, n ; f j = f j (t), 1, j m . Задана система из "n" дифференциальных уравнений первого порядка. 1 11 1 12 2 1 11 1 1 2 21 1 22 2 2 21 1 2 1 1 2 2 1 1 n n m m n n m m n n n nn n n nm m dx a x a x a x b f b f dt dx a x a x a x b f b f dt dx a x a x a x b f b f dt Начальные условия x i (t 0 ) = x i0 i= 1, n . t [t 0 ,T]. Эквивалентная запись в матричной форме: ; A , B ; x , f n n n m n m dx Ax Bf dt Начальные условия x(t 0 ) = x 0 Начнем с однородной системы dx Ax dt . Нетрудно видеть, что x(t)=e At u. Вектор "u" определяется из начальных условий u = e -At 0 x 0 , x(t)=e A(t-t 0 ) x 0 95 Переходим к решению неоднородной системы. Положим x(t) = e At v(t), где v(t) подлежащая вычислению вектор функция. После подстановки x(t) = e At v(t) в дифференциальное уравнение получим: 0 0 ( ) ( ) ( ); ( ); v(t)=v(t ) ( ) t At At At At A t dv dv Ae v t e Ae v t Bf t e Bf t e Bf d dt dt Полученное выражение подставляем в x(t) 0 ( ) 0 x(t)=e v(t ) ( ) t At A t t e Bf d Из начальных условий имеем 0 0 0 0 0 0 ( ) ( ); v( ) ( ) At At x t e v t t e x t . Таким образом, решение неоднородной системы уравнений имеет вид: (*) Частные случаи. 1. dx Ax dt x(t 0 ) = x 0 0 ( ) 0 x(t)=e x(t ). A t t 2. 0 ( ) 0. x t 0 ( ) x(t)= ( ) t A t t e Bf d 3.t 0 = 0. ( ) 0 x(t)=e x(0) ( ) t At A t e Bf d 4. 0 ( ) ( ) . t 0. dx t Ax t b dt ( ) 0 x(t)=e x(0) t At A t b e d Здесь b вектор размера n. После замены переменных ρ =t – τ 0 ( ) 0 0 t t A t A A t e d e d e d . Если существует А -1 , то 1 0 ( ) t A At e d e E A и 1 x(t)= (0) ( ) At At e x e E A b Выводы. 1.Для нахождения точного решения системы линейных дифференциальных уравнений необходимо произвести расчет матричной экспоненты e At 0 0 ( ) ( ) 0 x(t)=e x(t ) ( ) t A t t A t t e Bf d 96 2. Использовать представление e At в форме степенного ряда неэффективно , т. к. для достижения приемлемой точности придется удерживать недопустимо большое число членов ряда. 4. Использование формулы Лагранжа – Сильвестра или разложения Жордана требует расчета собственных чисел. Эта задача плохо обусловлена при большом размере матрицы. Перечисленные проблемы решаются заменой исходной системы эквивалентной системой разностных уравнений. 9.11. Решение системы разностных уравнений. П е р е х о д к с и с т е м е р а з н о с т н ы х у р а в н е н и й. Выбираем величину шага H. На оси "t" независимой переменной формируем систему узлов t 0 , t 1 , …t k , … с шагом H. Будем вычислять решение шаг за шагом. Пусть решение вычислено в точке t k и требуется вычислить x(t k +H). Тогда, воспользовавшись (*), имеем ( ) x(t +H)=e x(t ) ( ) k k k t H A t H AH k k k t e Bf t d Здесь предполагается, что на промежутке t [t k ,t k +H] f(t) ≈f(t k ). Произведя замену переменных ρ=t k +H-τ в подынтегральной функции, получим 1 1 1 0 x =e x ( ) ( ); x(t ) , x(t ) , ( ) H AH A k k k k k k k k k e Bd f t x x f t f Система разностных уравнений выглядит следующим образом x k+1 = Sx k + Gf k ; S = e AH , G= 0 H A e Bd . Пошаговое решение этой системы дает следующий результат x 1 = Sx 0 + Gf 0 x 2 = Sx 1 +Gf 1 = S 2 x 0 +SGf 0 +Gf 1 x 3 == Sx 2 +Gf 2 = S 3 x 0 +S 2 Gf 0 +SGf 1 + Gf 2 1 1 0 0 k k k l k l l x S x S Gf |