Вычислительная математика лекции

соответствии с этим обобщением многочлен Лагранжа – Сильвестра как и прежде можно записать в виде формулы:

s

r( ) [ f ( k ) k1 ( ) f '( k ) k 2 ( ) ... f (mk 1) ( k ) k ,mk ( )].

k 1

Функция φk,j(λ) (j=1,2,…,mk; k = 1,2 …,s) не зависит от f(λ) и представляет собой интерполяционный многочлен Лагранжа – Сильвестра для функции, у которой все значения на спектре матрицы А равны нулю за исключением одного f(j-1)(λk), равного 1.

Остаточный член многочлена Лагранжа – Сильвестра и в этом случае

равен r( )=f ( ) r( ) ( ) f (m) ( ) . . Поэтому в этих условиях

m!

попрежнему f(A) = r(A).

П р и м е р.

Пусть ψ(λ) = (λ-λ1)2 (λ-λ2)3. Тогда

r( ) f ( 1) 11( ) f '( 1) 12 ( ) f ( 2 ) 21( ) f '( 2 ) 22 ( ) f (2) ( 2 ) 23 ( ).

При расчете матричной функции можно исключить этап формирования многочлена Лагранжа – Сильвестра, что демонстрируется следующим примером.

П р и м е р 2 Вычислить f(At) =eAt , f(λ)=eλt для матрицы

2		1	1
	0	1		2
А=			1	с минимальным многочленом ψ(λ) = (λ-1) (λ-2).
	1	1
			1

Здесь t – скалярный коэффициент. Расчетная формула для данного случая f(At)=f(1)Z1+f'(1)Z2+f(2)Z3, где Z1, Z2, Z3 искомые матрицы, не зависящие от f(λ). Для их однозначного определения необходимо сформировать три условия.

(1)Пусть f(λ)= λ0=1. Тогда Z1+ Z3=E.

(2)Пусть f(λ)=λ-1. Тогда Z2+ Z3=A-E

91

(3) Пусть f(λ)= (λ-1)( λ-2) . Тогда - Z2=(A-E)(A-2E).

Решая систему трех линейных алгебраических уравнений,

			(2E A) A
получаем
	Z	( A E)( A 2E)
			( A E)		2
			1	0 0				1		1 1			0 0		0
f ( At) f		(1)		0 0			f '(1)		1	1 1		f (2)		1 1	0
			1
				1 1					0	0 0				1 1	0
			1

	f (1)	f	'(1)
		'(1) f		(2)
f (1) f
	f (1) f		(2)

f '(1)	f '(1)
f '(1) f (2)	f '(1)			t
			e
f (1) f (2)	f (1)

et tet					tet			tet
te	t	e	2t	e	2t	te	t	te	t
									.
et e2t				et e2t				et

9.8.Свойства матричной экспоненты.

1.

deAt

AeAt eAt A.

dt

Доказательство.

d

[E At

A2t2

A3t3

...] A A2t

A3t2

...

dt

2!

3!

2!

A[E At

A2t2

A3t3

...] AeAt .

2!

3!

2.

t eA d

t [E A

A2 2

A3 3

A 2

A2 3

...]d [E

...]

t0

0

0

2!

3!

2!

3!

Et

At2

A2t3

... [E At

A2t2

A3t3

... E]A 1 (eAt E) A 1.

2!

3!

2!

3!

3. Функции от одной матрицы коммутируют.

(eAt E)A 1 A 1 (eAt E) .

9.9.Использование канонической формы Жордана при вычислениях матричных функций.

Будем вычислять матричную функцию eAt.

92

eAt E At A2t2 A3t3 ...

2! 3!

Используя подобное преобразование, получим At=SJtS-1, где J

форма Жордана.

Aktk = (SJtS-1)k = (SJtS-1 SJtS-1… SJtS-1) = S(Jt)k S-1. Таким

образом,

eAt E S Jt S 1

S Jt 2 S 1

S Jt 3 S 1

...

2!

3!

S[E Jt

Jt 2

Jt 3

...]S 1 SeJt S 1.

2!

3!

Jt diag(J

t, J t,..., J

t);

(Jt)k diag((J

t)k , (J t)k ,..., (J

t)k );

0

1

s

0

1

s

eJt E Jt

Jt 2

Jt 3

... diag(eJ0t , eJ1t ,..., eJst ).

2!

3!

Если Ji клетка размером k для λ = λi, то Ji = λiE+J(1), где J(1)

матрица размера "к", у которой элементы первой наддиагонали равны 1, а остальные элементы нулевые. Согласно правилам умножения матриц (J(1))p = J(p) , где J(p) матрица размера "к", у

которой элементы p-ой наддиагонали равны 1, а остальные элементы нулевые. При p≥k матрица J(p) становится нулевой. Матричная экспонента

eJit e( i Ek Jk (1) )t e it E eJk	(1)t e it E (E	J (1)t J (2)t 2		1		... J (k 1)t k 1	1	) .
k	k k	k	k	2!		k	(k 1)!

Таким образом,

93

0x0,

	1	t		t 2		....		tk 1
			2!					(k 1)!
								tk 2
	0	1	t			....
								(k 2)!
eJit e it ...		...	...			...	...
	0	0	0			....	t
	0	0	0			....	1

9.10.Дифференциальные уравнения и матричная экспонента.

Рассматриваются обыкновенные линейные дифференциальные

уравнения с постоянными коэффициентами. Обозначим xi = xi(t), i=

1, n ; fj

= fj(t),

j 1, m . Задана система из "n" дифференциальных

уравнений первого порядка.

dx1

a x

a

x

.... a x

b

f .... b

f

m

dt

11 1

12 2

1n n

11 1

1m

dx2

a x

a

x

.... a x

b f

.... b

f

m

dt

21 1

22 2

2n n

21 1

2m

.............

dxn

a x

a

x

.... a x

b f

.... b

f

m

dt

n1 1

n2 2

nn n

n1 1

nm

Начальные условия xi(t0) = xi0

i=1, n . t

[t0,T].

Эквивалентная запись в матричной форме:

dx

Ax Bf ; A n n , B n m ; x n , f m . Начальные условия x(t0) =

dt

x0.

Начнем с однородной системы

dx

Ax . Нетрудно видеть, что

dt

x(t)=eAtu. Вектор "u" определяется из начальных условий u = e-At x(t)=eA(t-t0)x0.

94

Переходим к решению неоднородной системы. Положим x(t) = eAtv(t), где v(t) подлежащая вычислению вектор функция. После подстановки x(t) = eAtv(t) в дифференциальное уравнение получим:

AeAt v(t) eAt

dv

AeAt v(t) Bf (t);

dv

e At Bf (t);

v(t)=v(t0 )

t e A Bf ( )d .

dt

dt

t0

Полученное выражение подставляем в x(t)

t0

x(t)=eAt v(t

0 )

t eA(t ) Bf ( )d . Из начальных условий имеем

x(t ) eAt0 v(t );

v(t ) e At0 x(t

) . Таким образом, решение

0

0

0

0

неоднородной системы уравнений имеет вид:

A(t t0 )

t

A(t )

(*)

x(t)=e

x(t0 )

t0

e

Bf ( )d .

Частные случаи.

1.		dx	Ax.		x(t ) = x . x(t)=eA(t t0 ) x(t				).
		dt			0		0	0
2. x(t0 ) 0.					x(t)=	t	eA(t ) Bf ( )d .
						t0
3.t0 = 0. x(t)=eAt x(0) 0t eA(t ) Bf ( )d .
4.	dx(t)			Ax(t) b. t0 0.				x(t)=eAt x(0) b		t eA(t )d .Здесь b вектор
		dt								0
размера n.
После замены переменных ρ =t – τ
0t eA(t )d t0 eA d 0t eA d . Если существует А-1, то
0t eA d (eAt E) A 1 и							x(t)= eAt x(0) (eAt E)A 1b.
Выводы.
1.Для нахождения точного решения системы линейных

дифференциальных уравнений необходимо произвести расчет матричной экспоненты eAt.

95

2. Использовать представление eAt в форме степенного ряда неэффективно , т. к. для достижения приемлемой точности придется удерживать недопустимо большое число членов ряда.

4. Использование формулы Лагранжа – Сильвестра или разложения Жордана требует расчета собственных чисел. Эта задача плохо обусловлена при большом размере матрицы.

Перечисленные проблемы решаются заменой исходной системы эквивалентной системой разностных уравнений.

9.11.Решение системы разностных уравнений.

П е р е х о д к с и с т е м е р а з н о с т н ы х у р а в н е н и

й.

Выбираем величину шага H. На оси "t" независимой переменной формируем систему узлов t0, t1, …tk, … с шагом H.

Будем вычислять решение шаг за шагом. Пусть решение вычислено в точке tk и требуется вычислить x(tk+H). Тогда, воспользовавшись (*),

имеем x(tk +H)=eAH x(tk ) ttkk H eA(tk H ) Bf (tk )d . Здесь предполагается,

что на промежутке t [tk,tk+H] f(t)≈f(tk). Произведя замену переменных ρ=tk+H-τ в подынтегральной функции, получим

xk 1 =eAH xk ( 0H eA Bd ) f (tk ); x(tk 1 ) xk 1, x(tk ) xk , f (tk ) fk .

Система разностных уравнений выглядит следующим образом xk+1 = Sxk + Gfk; S = eAH, G= 0H eA Bd .

Пошаговое решение этой системы дает следующий результат

x1 = Sx0 + Gf0

x2 = Sx1 +Gf1 = S2x0 +SGf0+Gf1

x3 == Sx2 +Gf2 = S3x0 +S2 Gf0 +SGf1+ Gf2

k 1

xk S k x0 S k 1 lGfl

l 0

96

Р е к у р р е н т н а я а п п р о к с и м а ц и я м а т р и ч н о

йэ к с п о н е н т ы.

Решение системы разностных уравнений требует вычисления

матриц S = eAH и G= 0H eA Bd , которые не зависят от величины tk.

Их вычисление может быть проведено с использованием рекуррентных формул. Пусть h=H/2N, где N целое число. Вычислив exp(Ah), можно найти exp(AH), применяя последовательно формулы

exp(2Ah)= exp(Ah)exp(Ah)

……………………………….

exp(AH)= exp(2NAh)=exp(2N-1Ah)exp(2N-1Ah)

Таким образом, получаем рекуррентное соотношение φk+1=φk2,

где φk=exp(2кAh), k=0,1,…,N-1, причем φN=exp(AH)=S.

Построим теперь рекуррентное соотношение для вычисления

матрицы

G= 0H eA Bd .

Докажем следующее соотношение

02h exp( A )d 0h exp( A )d (E exp( Ah))

. С этой целью преобразуем правую часть

0h exp( A )d (E exp( Ah))

0h exp( A )d 0h exp( A Ah)d

0h exp( A )d h2h exp( A )d 02h exp( A )d ;0h exp( A Ah)d h h2h exp( A )d

Таким образом, g1=(E+φ0)g0, где

97

h exp( A )d B g ;

2h exp( A )d B g ;

0

0

0

1

gk 1 gk (E k ); k=0,1,...,N-1;

2k h

2k 1 h

0

exp( A )d B

gk

; 0

exp( A )d B gk 1;

gN G.

О в ы б о р е ш а г а h.

Шаг h должен выбираться достаточно малым для достижения

приемлемой точности аппроксимации при разумном числе l

учитываемых слагаемых в разложении экспоненты. Обозначим φ0

сумму l

учитываемых слагаемых

( Ah)l

( Ah)l 1

exp( Ah)

0

...

(l 1)!

l!

R

(|| Ah ||)l

(|| A || h)l 1

|| R ||

...

l!

(l 1)!

(|| Ah ||)

l (1

|| A || h

(|| A || h)2

...)

(h || A ||)

l

,(

h || A ||

1).

l 2

l!

l

l!

l

w

1

1

1 h||A||/l

Таким образом, погрешность аппроксимации

(h || A ||)l

. При

l !

заданной величине δ должно быть h<

l

l!

.

||

A ||

Шаг дискретности H=2Nh выбирается из условия наблюдения процесса, описываемого дифференциальным уравнением. При этом необходимая точность обеспечивается малым значением начального шага h. Использование рекуррентных соотношений дает выигрыш по количеству операций в ≈2N раза в сравнении с применением малого шага h на всем диапазоне решения.

98

9.12.Матричные нормы.

Определение матричной нормы эквивалентно определению векторной нормы. Именно, функция f, отображающая элементы множества Rm×n в элементы множества вещественных чисел R,

является матричной нормой, если

f(A)=||A||≥0	A	Rm×n	(||A||=0 A=0),
||A+B|| ||A||+||B||		A, B	Rm×n ,
||αA||=|α| ||A||		α R,	A Rm×n .
Чаще всего используются нормы Фробениуса (евклидова
норма)
					m n
			|| A ||F	aij2
					i 1 j 1
и p – нормы
		||A||p= sup		|| Ax ||p
				|| x ||p
			x 0

О p – нормах говорят, что они подчинены соответствующим

векторным нормам. Для таких норм выполняется условие

||Ax||p ||A||p ||x||p.

Введенные матричные нормы удовлетворяют

мультипликативному условию

||AB|| ||A|| ||B||.

При выполнении этого условия говорят, что нормы матриц A и

B взаимно согласованы.

Некоторые свойства матричных норм.

Для A Rm×n имеем:
||A||2 ||A||F		n || A ||2 ,
max | aij ||| A ||2
	mn max | aij |,
i, j				i, j
99

m

||A||1 =

max | aij | ,

1 j n i 1

n

||A||∞ = max | aij |,

1 i m i 1

1

|| A ||

|| A ||

m || A || ,

n

2

1

|| A ||

|| A ||

n || A || ,

m

1

2

1

|| A ||2 max j ( AT A).

1 j n

9.13.Устойчивость решений обыкновенных дифференциальных и разностных уравнений.

9.13.1. Основные понятия теории устойчивости.

Факт устойчивости устанавливается по результатам анализа решений исходной и возмущенной систем уравнений. Пусть исходная и возмущенная системы имеют вид:

исх. система dxdt f (x(t),u(t)), нач. усл. x(t0 ) x0

возм. система dx * f (x * (t),u(t)), нач. усл. x*(t0 ) x0 * dt

Обозначим x(t)-x*(t)=v(t).

О п р е д е л е н и е. Решение x(t) исходной системы устойчиво

по Ляпунову, если для любого ε > 0 существует δ > 0, такое что || v(t)

|| ε для всех t > t*, если || x0 – x0* || δ.

100