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

ν(t)

δ

ε

t*t

-ε

Если решение x(t)

-устойчиво по Ляпунову

-норма || v(t) || → 0 при t→ ∞,

то это решение асимптотически устойчиво.

Аналогичная ситуация имеет место при анализе устойчивости решения систем разностных уравнений

Д о п о л н и т е л ь н ы е з а м е ч а н и я. 1.Устойчивость по Ляпунову эквивалентна выполнению

условия

max || x(t) x *(t) || K || x0 x0 * ||

t* t

Скаляр K играет роль числа обусловленности.

Это условие получается, если положить ε=Кδ.

2.Возмущенная система может возникать при внесении погрешностей как в начальные условия, так и в правую часть.

Можно показать, что в этом случае для устойчивой по Ляпунову системы выполняется условие

101

max || x(t) x *(t) || K[|| x0 x0* || max || ( ) ||] . Здесь ψ(t) малая

t* t

t0

вектор-функция для дифференциальной системы или малая сеточная вектор-функция для разностной системы, добавленная в виде слагаемого к правой части.

9.13.2..Анализ устойчивости решения систем линейных

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

коэффициентами.

Исходное уравнение имеет вид

dxdt Ax Bf (t); A n n , B n m ; x n , f m .

	x(t0) = x0, t [t0, T]. Полагая как и прежде x(t)-x*(t)=v(t), x(t0)-
x*(t )=v(t ) = v ,получаем			dv	Av, v(t		) v . Здесь x*(t) решение
					0
0	0	0	dt			0

возмущенной системы, полученное из исходной при новых

начальных условиях.

Таким образом, для установления факта устойчивости нужно

проанализировать решение однородной системы дифференциальных уравнений. Это решение вычислено ранее v(t)=eA(t t0 ) v(t0 ).

Для вычисления матричной экпоненты воспользуемся

полиномом Лагранжа – Сильвестра.

Пусть минимальный аннулирующий многочлен матрицы А

имеет вид

ψ(λ) = (λ-λ1)m1 (λ-λ2)m2…(λ-λs)ms, m1+ m2+…+ ms=m.В данном

случае	f(λ) = e (t t0 ) .
Тогда

102

s

eA(t t0 ) [Zk1 k 1

s

[Zk1e k (t t0 ) k 1

s

eA(t t0 ) [Zk1 k 1

f ( k ) Zk 2 f '( k ) ... Zkmk f (mk 1 1) ( k )]

Zk 2 (e k (t t0 ) ) ' ... Zkmk (e k (t t0 ) )(mk 1 1) ]

Zk 2 (t t0 ) ... Zkmk (t t0 )mk 1 ]e k (t t0 ) .

Здесь Zki матрицы с постоянными коэффициентами.

Рассмотрим три случая:

1.Reλk 0 (k = 1,2,…,s), причем для тех λk, для которых Reλk=0, mk = 1 (т. е. чисто мнимые характеристические числа являются простыми корнями минимального многочлена).

2. Reλk<0 (k = 1,2,…,s).

3.Хотя бы для одного k имеем Reλk>0 либо Reλk=0, но mk >

1.Тогда из формулы

s

eA(t t0 ) [Zk1 Zk 2 (t t0 ) ... Zkmk (t t0 )mk 1]e k (t t0 ) следует, что в

k 1

случае (1) элементы матрицы Q(t)= eA(t t0 ) ограничены на интервале

[t0,∞), в случае (2) Q(t)→0 при t→∞, а в случае (3) элементы матрицы

Q(t) не ограничены на интервале (t0,∞).

Таким образом, в случаях (1) и (2) решение устойчиво по Ляпунову причем в случае (2) устойчиво асимптотически. В третьем случае неустойчиво.

9.13.3.Анализ устойчивости решения систем линейных

разностных уравнений с постоянными коэффициентами.

Исходное уравнение имеет вид

xk 1 Sxk Gfk ; S n n , G n m ; x(tk )=xk n , f(tk )=fk m. k [0, kmax ].

Начальные условия: x(t0) = x0 . Если возмущать начальные условия, обозначив, как и прежде, xk-x*k=vk, x0-x*0= v0,получаем vk 1 Svk .

103

Таким образом, для установления факта устойчивости нужно проанализировать решение однородной системы разностных уравнений. Это решение получено ранее vk = Skv0

Пусть минимальный аннулирующий многочлен матрицы S

имеет вид

ψ(λ) = (λ-λ1)m1 (λ-λ2)m2…(λ-λs)ms, m1+ m2+…+ ms=m. В данном случае f(λ) = λk.

Тогда

s

S k [Zi1 i k Zi 2k ik 1 ... Zimi (k(k 1)...(k mi 1)) ik mi ].

i 1

Рассмотрим три случая:

1.| λi | 1 (i = 1,2,…,s), причем для тех λi, для которых | λ i |=1, mi = 1 (т. е. соответствующие характеристические числа являются простыми корнями минимального многочлена).

2. | λi |<1 (i = 1,2,…,s).

3.Хотя бы для одного i имеем | λi |>1 либо | λi |=1, но mi > 1.Тогда из формулы

s

S k [Zi1 i k Zi 2k ik 1 ... Zimi (k(k 1)...(k mi 1)) ik mi ] следует, что

i 1

в случае (1) матрица Sk ограничена на интервале k [0,∞), в случае

(2) Sk →0 при k→∞, а в случае (3) матрица Sk не ограничена на интервале k [0,∞).

Таким образом, в случаях (1) и (2) решение устойчиво по Ляпунову, причем в случае (2) устойчиво асимптотически. В третьем случае неустойчиво.

104

10.Численные методы решения обыкновенных дифференциальных уравнений

10.1. Постановка задачи.

Рассмотрим численные методы решения известной задачи Коши:

dx	f (t, x),	t [t		,t	],	x(t		) x ,	x n.
			0				0
dt				k				0

Прежде чем начинать вычисления, необходимо установить их принципиальную возможность, т. е. факт корректности задачи.

Иными словами нужно установить 1) существование решения, 2)

единственность решения, 3)непрерывную зависимость погрешности решения от погрешности исходных данных (в данном случае начальные условия и правая часть). Ответ на эти вопросы дает теорема Коши существования и единственности.

Если f(t,x) непрерывна в некоторой области вокруг точки

(t0,x0), тогда существует по крайней мере одно решение уравнения x'=f(t,x), принимающее при t=t0 значение x0, определенное и непрерывное в некотором интервале около t0.

Если, кроме того, в этом интервале выполнено условие Липшица

|| f(t,x1) – f(t,x2) || L|| x1 – x2 || причем L не зависит от t, x1 и x2,

то это решение единственно и является непрерывной функцией от x0.

Для оценки значения постоянной Липшица L разложим в ряд Тэйлора функцию f(t,x2) относительно точки x1(x1, x2 Rn).

		f			1
f (t, x ) f (t, x )			(t, )(x x );			; [x , x ].
2	1	x	2 1				i	1i 2i
					n

105

			f1		f1
			x		x
	f		1		n
Здесь		(t, )				матрица Якоби.
	x		fn		fn
			x		x
			1		n

Таким образом, || f (t, x2 ) f (t, x1) || || fx (t, ) || *|| (x2 x1) || .

Все численные методы связаны с переходом от дифференциальной системы к разностной схеме и вычислением приближенных решений в форме сеточной функции. В дальнейшем будем обозначать x(tk) решение дифференциальной системы в узле tk,

а xk решение разностной схемы в том же узле, tk = t0 + kh, где h шаг сетки, k = 0, 1, 2, … .

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

10.2.Метод аналитического продолжения ( метод рядов).

Предполагаем, что x, f(t,x) Rn.

Предположим, что установлен факт корректности задачи и функция f(t,x) непрерывно дифференцируема требуемое число раз.

Тогда решение x(t) можно разложить в ряд Тэйлора около центральной точки tk.

106

x(t

) x(t

) hx '(t

)

h2

x ''(t

) ...

hs

x

( s) (t

)

hs 1

x( s 1) (

)

k 1

k

k

k

k

(s 1)!

k

2!

s!

x(t

k

) h[ (t

k

, x(t

k

), h) r (t

k

, x( s 1) (

k

), h)],

k

r (t

, x( s 1) (

), h)

hs

x( s 1)

(

),

[t

,t

];

k

k

k

k

k

k 1

k

(s 1)!

(t

, x(t

), h) x '(t

)

h

x ''(t

) ...

hs 1

x( s) (t

);

k

k

k

k

k

2!

s!

Следовательно, эквивалентное разностное уравнение имеет вид

x(tk 1 ) x(tk )				(tk		, x(tk ), h) rk (x				( s 1)	( k ), h) (*)
h
						hs
r (x( s 1)	(	k	), h)				x( s 1) (	k	).
k					(s 1)!

Разностное уравнение метода (формулу метода) мы получим,

если пренебрежем погрешностью rk .

xk 1 xk	(t	, x , h) (**)
h	k	k

(tk , xk , h) x 'k h x ''k ... hs 1 x( s)k . 2! s!

Применение метода связано с использованием функции φ,

определяемой на каждом шаге через производные до s – ого порядка включительно. Это предполагает использование соотношений:

для производной первого порядка x'(tk)=f(tk,xk);

для производной второго порядка x''(tk)=

f					f			x '(t			) и т. д. с последующим усложнением
t					x			k
		t t	,x x(t	)			t tk ,x x (t k )
		k	k
выражений. Здесь под								f			понимается матрица размером n×n,
								x		t tk , x x ( tk )

называемая матрицей Якоби

107

	f1		f1
	x		x
	1		n
A
	fn		fn
	x		x
	1		n

Остается открытым вопрос о правомерности пренебрежения остаточным членом rk.

Назовем глобальной погрешностью величину E(h)= max|| x( tk)- xk|| (k=1,2,..). Численный метод решения задачи Коши называют сходящимся, если для него E(h)→ 0 при h→ 0. Принято говорить, что метод сходится с p-ым порядком точности (или имеет p-ый порядок точности), если для погрешности справедлива оценка E(h) Chp.

Очевидно, чтобы быть работоспособным, метод должен быть сходящимся.

Малая сеточная функция rk, добавленная к правой части уравнения (**), преобразует его в разностную систему (*). Её решение совпадает в узлах tk с решением исходного дифференциального уравнения. Величина rk выполняет роль погрешности аппроксимации. Говорят, что разностное уравнение метода (**) аппроксимирует исходное дифференциальное уравнение,

если при h→0 rk→0. В нашем случае в результате выполнения этой операции уравнение (*), а значит и уравнение (**) совпадут с исходным дифференциальным уравнением, т. е. аппроксимация имеет место. Погрешность аппроксимации

				hs			). пропорциональна hs. В таком случае
r (x( s 1)	(	k	), h)		x( s 1) (	k
k				(s 1)!

говорят, что уравнение метода (**) аппроксимирует дифференциальное уравнение с порядком точности s.

108

Для того чтобы разностная схема метода была сходящейся,

остается доказать её устойчивость, т. е. ограниченную чувствительность к изменениям исходных параметров (в данном случае правой части).

Действительно, пусть имеет место факт аппроксимации с

порядком точности s : max || rk || C1hs, C1	1		x(s 1)	( k ). .
	(s
k		1)!

Добавляя к правой части разностной системы (**) rk, получим возмущенную систему(*), решение которой совпадает с решением исходного дифференциального уравнения (x(tk), k=0, 1, 2, …) Если разностная система устойчива, то

||x(tk) – xk ||		K max || rk || KC1 hs = C hs , k=1,2,… То есть
			k
E max || x(t		) x	|| Chs , C не зависит от k и h.
k	k	k

Иногда удобнее использовать л о к а л ь н у ю погрешность L.

Локальной погрешностью называется погрешность на очередном k-

ом шаге в предположении, что значение xk рассчитывается по формуле метода подстановкой в неё точных решений на предыдущих шагах. Из формул (*) и (**) следует Lk+1(tk,x(s+1)(ξk)) = h rk = αk hs+1 (в (**) заменяем xk на x(tk) и затем вычитаем из(*); в результате получаем Lk+1=x(tk+1)-x(tk)=hrk=αkhs+1).

10.2.1.Методы Эйлера.

На практике находят применение следующие две модификации метода первого порядка ( в разложении учитываются слагаемые до первой производной включительно).

Явный метод Эйлера xk+1 = xk + hf(tk ,xk). Неявный метод Эйлера xk = xk+1 - hf(tk+1 ,xk+1).

xk+1 = xk + hf(tk+1 ,xk+1).

109

Для определения x(tk) в последней формуле нужно использовать разложение x(t) около центральной точки tk+1 .

Оценим устойчивость решений, получаемых этими методами.

Как и прежде возмущенная система получается в результате изменения начальных условий.

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

Обозначим v =x – x*. Вычитая из исходной системы

возмущенную, получим

dvdt f (x(t),t) f (x * (t),t), нач. усл. v0 x0 x0 *

Воспользуемся формулой разложения функции нескольких переменных в ряд Тэйлора с остаточным членом:

f(x(t),t) - f(x*(t),t) = = f ((x* v),t) (x(t) – x*(t)). Здесьx

	f ((x* v),t)	матрица Якоби; ϴ Rn×n, ϴ=diag( ϴi ), 0 ϴi
	x
i= 1, 2 ,…, n.

Тогда	dv			f ((x* v),t)		v(t); v		x		x *. После
							0
	dt			x					0	0
подстановки в		f ((x* v),t)			= A Rn×n				начальных условий
				x

1,

это

уравнение становится однородным дифференциальным уравнением с постоянными коэффициентами. Его используют в качестве тестового при анализе устойчивости.

Таким образом, система тестовых уравнений в этом случае

имеет вид v'(t) = Av(t). Пусть исходная задача асимптотически

устойчива, то есть Re(λi(A)) < 0; i=1,2,…,n. Естественно потребовать,

чтобы система разностных уравнений метода также была

асимптотически устойчивой.

110