16. МАТРИЧНА ЕКСПОНЕНТА ТА її ЗАСТОСУВАННЯ ДЛЯ РОЗВ'ЯЗАННЯ РІВНЯНЬ СТАНУ

Теорія диференціальних рівнянь настільки глибоко інтегрувала в себе матричні методи, що було введено спеціальне поняття - матрична експонента. Застосування цього поняття дозволяє в загальному аналітичному вигляді про інтегрувати систему лінійних диференціальних рівнянь.*

Якщо А є квадратна матриця, то запис ЄА позначає так звану матричну експоненту (матричний експоненціал). Ця матриця, в свою чергу, має той же порядок, що і матриця А, і для її обчислення можна застосувати формулу розкладання в ряд Тейлора, яка співпадає за формою з випадком скалярної експоненти, представленої рядом Тейлора:

Є = 1 + А + + +.... = ^=о— (16.1)

У цьому виразі І -це одинична матриця , порядок якої співпадає з порядком А.

Існує багато методів обчислення матричного експоненціалу , які випливають з матричної теорії , однак найпростішим методом є саме обчислення за формулою (16.1). При практичних обчисленнях кількість врахованих членів розкладання може визначатися шляхом оцінки того уточнення, яке вносить наступний член у величини елементів матричного експоненціалу.

Для подальшого застосування матричної експоненти слід визначити операції її диференціювання та інтегрування.

Для операції диференціювання здійснимо почленне диференціювання правої частини виразу (16.1):

А2 х /2    А" х /3 _ Акік

(еАІ)=А+(2А¥2!)+ (3А¥/3!)+ +(4А4і3/4!)+...=А+А21+ (А¥/2!)+

+(А4і3/3!)+...=А(і+Аі+(А21/2!)+ (А¥/3!)+...= АеА . (16.2)

Таким чином, щоб продиференціювати матричний експоненціал, слід помножити його на матрицю А. При цьому можна довести, що дозволяється переставити матриці А і е :

— (еАІ)= еАІА (16.3)

СІІ

Якщо показник матричного експоненціалу містить позначку "мінус ", то

— (е - А1)= -Ае - А1= -еАІА. (16.4)

Для   того,   щоб   проінтегрувати   матричну   експоненту, розглянемо рівняння(16.2) і візьмемо інтеграли від обох частин:

Г й ( м \]

Я .1 И)

—х е ;J сіі= |а еАІси (16.5)

Взяття інтегралів приводить до співвідношення

еА1= А ї еА1((1 . (16.6)

Домножимо обидві частини на А-1 і отримаємо :

їеА ((1= А-1 еА1 . (16.7) Тут також дозволяється переставляти матриці А-1 і еА1 :

|еА (1= еА1 А-1 . (16.8)

При наявності позначки "мінус " в показнику експоненти:

|е~м (1= - А-1 е--А1= -е-А1А-1 (16.9)

Тепер розглянемо систему неоднорідних лінійних диференціальних рівнянь, складених за методом змінних стану для електричної системи :

— = АХ+В<2. (16.10)

Будемо мати на увазі, що Х і 0 є функціями часу, тобто Х=Х(1) і 0=0(1). Домножимо обидві частини (16.10) на е-А1:

е-Ах—=е-АхАХ+ е-Ах В<2. (16.11) —х

Перенесемо перший доданок правої частини в ліву:

е-Ах—-е-АхАХ= е-Ах ВЦ. (16.12) —х

Ліва частина (16.12) являє собою похідну від добутку двох функцій и=е-Ах і V = Х(х), тому можна записати :

— \емХ(х)] = е -АхВЦ(Х). (16.13)

Зазначимо, що тепер Х і Ц записані як функції часу. Візьмемо визначений інтеграл від обох частин від 1=0 до 1=Т , де Т представляє крок інтегрування :

Т — т

ї — [є ~АхХ (X)— = | є -АхБ0(Х(16.14)

о —

Оскільки в лівій частині береться інтеграл від похідної, то можна записати:

е-Ах Х(х)

т

т г

-Аг

ї е-АхБд(Х)сіХ (16.15)

о

о

о

Підставляючи межі інтегрування, отримаємо :

еАТх(т) -х(0)= [Vаіб<2(іуь .

Запишемо в лівій частині лише перший доданок :

е-АТ х(т) = х(0)+ [V а'б<2(і уи

(іб.іб)

(16.17)

Тепер для того, аби в лівій частині залишилось лише Х(Т), домножимо обидві частини (16.17) на ЄАТ :

Рівняння (16.19) представляє аналітичне матричне роз'язання системи диференціальних рівнянь за методом змінних стану. Перший доданок відображає вільний рух системи, якщо Х(0) не дорівнює нулю, тобто система має запас електричної енергії. Другий доданок являє вплив на перехідний процес джерел енергії. В обчисленні цієї складової є найбільші складності, оскільки треба брати інтеграл від функції, залежної від (((і), що може бути будь

якою. Задача спрощується, коли вважати, що величина 0(1) протягом кроку

інтегрування є незмінною і дорівнює величині на початку кроку інтегрування :

Тоді при обчисленні інтегралу В( може бути винесено за позначку інтегралу:

АТ

Раціонально також винести і е   за позначку інтегралу, тоді треба визначити

еАТ \1 е~А(ЖБд = еАТ(-А -1) е -Аї| 0*В((= еАТА-1е "АТ|0 *В(( (16.21)

Х(т) = еАТ Х(0)+ еАТ [ е~АБШУЬ

АТ

Можна внести е    під знак інтегралу :

Х(т) = еАТ Х(0)+ ГеА(Г~')Б&і)сіі

(16.18)

(16.19)

0(1)= 0(0)= 0 = СОПБІ

 

(16.20)

Обчислимо окремо матрицю С, яка помножується на В((:

С=еАТА"1е -АТ|0 = еАТ(І-е -АТ)А"1= (е АТ-І) А"1

(16.22)

т

 

(16.23)

Таким чином, можна записати розв'язання у вигляді

Х(Т)= еАТ Х(0)+СВ0,

(16.24)де згідно з (16.22)

С=(еАТ-І )Л-1 , (16.25) а якщо враховувати (16.23), то

А Т2     А2Т3 А3Т4

С=Ц+ АТ-++АТ-+.... (16.26) 2!       3! 4!

Якщо перейти до рекурентних формул, то (16.24) можна подати у вигляді Х(к+1)= ЄАТ Хк+ СВ))к (16.27)

Наведена формула надає дві можливості для обчислення матриць, що входять до рекурентних формул. По-перше, можна використати обчислення матричної експоненти розкладанням в ряд :

еАТ = І+АТ+ А-Т^+— + 2! 3!

При цьому можна заздалегідь вказати кількість врахованих членів ряду або вести розрахунок доти, доки наступний член вже не буде суттєво уточнювати результат. Після того, як обчислено матрицю експоненту, матрицю можна обчислити за (16.25). Однак аналізуючи формулу (16.25), можна зазначити, що вона пов'язана з необхідністю здійснення операцій з одиничною матрицею А-1 . Можна довести, що експоненціал і матрицю С можна обчислювати й уточнювати паралельно. Тому в другому алгоритмі діють наступним чином. Обчисливши перший член експоненціалу, можна помножити його на Т і отримати перший член матриці С. Якщо помножити перший член

матриці С на А, отримаємо другий член експоненціалу. Помноживши цей член

Т

на — , отримаємо другий член матриці С. Таким чином, для отримання

наступного члена матриці С слід помножати останній обчислений член

ТТТ

експоненціалу на Т, —, —, — і т.д., а для отримання наступного члена

експоненціалу слід останній обчислений член матриці С помножати на А.

Серед функцій МЛТЬЛВ взагалі введена функція обчислення матричної експоненти ехрт(Л), якою можна користуватися при обчисленні коефіцієнтів рекурентних формул.

На рис. 16.1 наведено робочий аркуш МашСАЕ), на якому реалізовано

АТ

розрахунок з використанням матричної експоненти. Через е позначено матричний експоненціал, обчислення якого обмежено 5 членами розкладання в ряд. Матриця С обчислюється за формулою (16.25), а рекурентна формула (16.27) забезпечує розрахунок перехідного процесу протягом заданого інтервалу часу.

РОЗВ'ЯЗАННЯ  ЗА ДОПОМОГОЮ МАТРИЧНОЇ ЕКСПОНЕНТИ

Е := 10С

І := 10    Ю := 0.01

лл/

:= 10

Ь := 0.1

АЛЛ

С := 0.01

АЛЛ/

сік := 0.0005

С-И2 С А :=

I   Ь       Ь )

Г 0 -1

С

В :=

1 0

eigenvals (А)

I :=

-5.05 + 31.233ІЛ -5.05 - 31.233і)

А

-10 100 -10 -0.1

В

0 -100 10 0

31.233

10 01

0.201

0.198

еАТ := I + ік • А + ^ • А2 + *

1 2

3

- -А3

сік

еАТ

0.994888

1- 2-3 0.049872

1- 2-3 - 4

-А4

С

V -4.987188х 10 3 0.999825)

к := 0.. 10000

4.987313 х 10

С, = (еАТ - I) - А" 4

1.247872 х 10

5 ^

V-1.247872 х 10 6 4.999667 х 10 4 )

:= еАТ -

1 )

V   к )

+ С-В- ((

200

150

100

50

"50

л

 

 

 

[

 

 

-

і

 

 

0

 

9995

99.8002

 

9996

99.8002

 

9997

99.8002

 

9998

99.8002

 

9999

99.8002

 

10000

99.8002

 

10001

99.8002

 

0.5

к-сШ

 

0

9995

19.98002

9996

19.98002

9997

19.98002

9998

19.98002

9999

19.98002

10000

19.98002

10001

19.98002

Рис.16.1

КОНТРОЛЬНІ ЗАПИТАННЯ І ЗАВДАННЯ ДЛЯ САМОСТІЙНОЇ РОБОТИ

1. Як визначається матрична експонента?

2. Як обчислюються похідна і інтеграл від матричної експоненти?

3. Наведіть викладки до роз'язання рівнянь стану із застосуванням матричної експоненти.

4. Які способи можна використати для обчислення матриць в рекурентних формулах для розрахунків перехідних процесів із застосуванням матричної експоненти?

5. Поясніть структуру робочого аркуша МашСАО для реалізації методу матричної експоненти.