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

алгоритмов для расчета дополнительных начальных условий. В качестве стартовых методов может быть использован любой другой одношаговый метод, например один из методов РунгеКутты. Порядок сходимости стартового метода должен быть не ниже, чем у основного.

10.4. Свойство жесткости систем дифференциальных уравнений

Проблемы, связанные с решением "жестких" систем

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

Пример 1. Пусть матрица А размера 2×2 имеет 2 собственных

значения 1	2,	2 1 .	Решение	системы содержит две
составляющие	с экспонентами		exp( 2t) и	exp( t) . Интервал решения

определяется самой медленной экспонентой Т=1с. и обычно выбирается равным (3 5)T . Для построения графика решения вполне достаточно

выбрать шаг временной дискретности h 0.1 0.2 и тогда график решения будет содержать десять-тридцать точек. Условие численной устойчивости при использовании явного метода ломаных Эйлера связано с выполнением

неравенства h		i		max 2 , что полностью					соответствует поставленной
задаче.
Пример 2. Пусть					104	,		1	. Решение вновь содержит две
					1		2

экспоненты и, т.к. самая медленная осталась прежней, время наблюдения решения остается сохраняется t [0,3] . Но на этот раз на начальном участке имеется резко меняющаяся экспонента с малой постоянной времени 1 1 1 10 4 с . Для построения графика такого решения необходимо взять десяток точек с шагом h 10 5 , а затем, на оставшейся,

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

121

h 2 10 4 c., и весь процесс решения связан с выполнением десятков тысяч

шагов.

"Жесткими" называют системы дифференциальных уравнений, для решения которых характерны два участка. Первый из них малой продолжительностью nc называют пограничным слоем. Решение здесь обладает большими производными и изменяется очень быстро. Второй участок продолжительностью много больше первого (T nc ) обладает сравнительно малыми производными, и решение носит спокойный

характер.

Для линейных систем

x(t) Ax(t),	x(t	) x0 ,
	0

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

cond ( A)

A

A 1

, cond ( A) 1 или

max i ( A)

	k( A)						,	k( A)		1.
		min		( A)
				i
				i
Для нелинейных систем общего вида
	x(t) f (t, x), x(t								) x0	,
								0
жесткость	системы	часто						предопределяется			плохой

i

обусловленностью матрицы Якоби f в процессе решения.

x

На примере решения линейных систем убеждаемся и том, что время наблюдения решения Т определяется минимальным по модулю собственным значением матрицы А, а шаг интегрирования – максимальным. Тогда число шагов N прямо пропорционально числу обусловленности, что приводит к недопустимым затратам

T

1

, h

1

,

N

T

i

max

1.

h

i

min

i

max

i

min

122

В связи с указанными причинами следует избегать использования явных методов численного интегрирования. Для интегрирования жестких систем рекомендуется использовать неявный метод ломанных Эйлера или метод Адамса-Мултона второго порядка точности (неявный метод трапеций), устойчивость которых обеспечивается при выполнении условия Re(i ( A)) 0, i , вне зависимости от выбора величины шага интегрирования. Существенный недостаток этих методов малая точность. Среди неявных линейных многошаговых методов присутствуют A(α) устойчивые методы высоких порядков точности, пригодные для решения жестких систем. В качестве примера можно назвать семейство алгоритмов Гира. A(α) устойчивость означает, что область устойчивости метода распространяется на большую часть левой полуплоскости.

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

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

Система линейных алгебраических уравнений (СЛАУ), записанная в матричной форме, имеет вид

Ax = b, где A Rn×n ; x,b Rn. Требуется найти её решение с заданной точностью.

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

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

Затем, рассчитав коэффициенты обусловленности, выяснить возможность

123

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

Корректность задачи.

Проанализируем условия существования и единственности решения.

			x1		b1
Ax a1		an
			xn		bn

Из этой формулы следует, что вектор b есть линейная комбинация столбцов матрицы ai c весами xi (i= 1,2,…,n). Если ранг r(A) =n, т. е.

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

В случае вырожденной матрицы r < n, вектор столбцы принадлежат подпространству полного пространства. Это подпространство имеет размерность r. Если вектор b не принадлежит этому подпространству,

система несовместна и решение отсутствует. Если вектор b принадлежит подпространству столбцов, то существует бесконечное множество решений. Действительно, общее решение есть сумма решений однородного уравнения и частного решения. Известно, что однородная система с вырожденной матрицей имеет бесконечное число не нулевых решений. Вектор b в этом случае должен иметь n – r элементов,

являющихся заданной линейной комбинацией остальных r элементов.

Продемонстрируем ситуацию следующим примером. Используя метод исключения Гаусса, получаем

124

	2	2	1	b1			2	2	1	b1
	0	6	5	b			0	6	5	b
				2						2
	1	2	2	b			0	6	5	2b b
				3						3 1
		2		2			1	b1
		0		6			5	b
								2
		0		0			0	2b b b
								3	1 2
В данном случае r(A)=2. Для существования частного решения
необходимо выполнение условия 2b3-b1-b2=0. Неизвестную x3 полагают
равной произвольной константе,							x2 и x1	определяются в ходе обратной
подстановки.
Обусловленность задачи.

Введем понятия погрешностей вектора и матрицы. Пусть компоненты вектора и матрицы А получают малые приращения.

Возмущенная система имеет матрицу А*=А+ΔА и правую часть b*=b+Δb.

Ее решение x*=x+Δx. Назовем абсолютной погрешностью вектора и матрицы величины b=|| b*-b ||=|| Δb ||, ΔA=||A* - A||=|| А ||

соответственно, а величины δb =|| b*-b ||/ || b ||, δA = ||A* - A||/ ||A||

относительными погрешностями. Будет показано, что погрешности решения связаны с погрешностями исходных данных соотношениями

|| Δx|| = ν (|| Δb || +|| ΔA ||) и

|| x*-x ||/ || x || = νδ(|| b*-b ||/ || b || + || A*-A ||/ || A ||). Здесь ν и νδ

коэффициенты абсолютной и относительной обусловленности

(чувствительности ) решения задачи. Их величина зависит от свойств матрицы системы и будет вычислена позднее.

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

Вычислительная неустойчивость приводит к накоплению недопустимо больших погрешностей, а иногда даже к аварийному останову или

125

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

Любая хорошо обусловленная задача может быть испорчена вычислительно неустойчивым алгоритмом. С другой стороны самый совершенный алгоритм не справится с плохо обусловленной задачей.

При решении СЛАУ используют численные методы двух разновидностей: прямые и итерационные. Характерным признаком прямых методов является возможность получить точное решение (в

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

формирующий последовательность векторов x0, x1, …,xk,… таким образом, что xk → x при k → ∞, где x точное решение системы.

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

Прямые методы, как правило, используют при решении задач среднего размера ( n ≈ 100) с заполненной матрицей.

11.3.Обусловленность задачи решения системы линейных алгебраических уравнений.

Предположим сначала, что x решение системы Ax = b и x+Δx

решение системы A(x+Δx ) = b+Δb с измененной правой частью.

Следовательно, AΔx = Δb, Δx = A-1Δb, ||Δx|| ||A-1 || ||Δb|| (*).

Таким образом, в данных условиях абсолютное число обусловленности ν = ||A-1 ||.

С практической точки зрения полезнее использовать относительные погрешности. Из Ax = b следует ||b|| ||A|| ||x||, что вместе с (*) дает

126

|||Δx|| ||b||	A||	||A-1 || ||Δb|| |||x||		или (b≠0)
(|||Δx|| /|||x||)		A||	||A-1 || (||Δb|| /||b||).
Относительное число обусловленности					νδ = ||A||		||A-1 || .
Произведение ||A||			||A-1 || играет очень важную роль и называется
числом обусловленности матрицы A. Стандартное обозначение этого
числа cond(A).
Рассмотрим теперь случай, когда изменяются элементы матрицы, так
что возмущенное уравнение имеет вид
(A + ΔA)(x + Δx) = b или с учетом Ax=b
AΔx = b - (A + ΔA)x – ΔAΔx =- ΔA(x+Δx ) , -Δx							A-1 ΔA(x+Δx ) .
Следовательно,
||Δx	A-1	ΔA	(x+Δx )|| = cond(A)(			ΔA /	A ) ||(x+Δx )||,
так что
Δx /	(x+Δx )||		cond(A)(	ΔA / A	) .
Обычно ||x|| Δx||. Поэтому ||(x+Δx )||						≈ ||(x|| .
Если с погрешностью заданы как правая часть, так и матрица,
причем cond(A)δ(A) 1, можно доказать справедливость неравенства

δx cond(A) (δ(b) + δ(A)).

Во всех рассмотренных случаях относительное число обусловленности решения совпадает с числом обусловленности матрицы .

Полученные результаты нужно правильно интерпретировать. Если cond(A) мало, например, близко к 1, то малые относительные изменения данных обязательно приводят лишь к малому изменению решения. С

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

Отметим следующие свойства числа обусловленности.

1. Для единичной матрицы cond(E) = 1. Это следует из того, что

E-1 =E и ||E||=1.

127

2.	Из условия ||E||=1 и из равенства E=AA-1 следует
1 = ||E|| ||A-1 || ||A|| = cond(A).
3.	Число обусловленности матрицы не меняется при умножении

матрицы на произвольное число α≠0, Это следует из соотношений (αA)-1

=α-1 A-1 и ||αA|| = |α| ||A||.

На практике влияние большого числа обусловленности зависит от

точности данных и длины слова используемой ЭВМ.

Проиллюстрируем наши результаты на примере. Рассмотрим

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

		1.03		0.991				2.51	, A 1	87.4	91.8
A		0.991		0.943		, b=					95.4
								2.41		91.8
cond			( A) || A || ||		A 1 ||		2.021*187.2 378

Решением системы является x1 ≈ 1.981, x2 ≈ 0.4735. Правая часть системы в лучшем случае известна с точностью 0.005, если считать, что ее элементы получены, например, округлением "истинных" значений при вводе в память трехзначной десятичной ЭВМ. Возмутим каждую из компонент вектора b на величину 0.005, положив b*=(2.505, 2.415)T .

Решением системы теперь является x*=(2.877, -0.4629)T . Решение оказалось полностью искаженным. Относительная погрешность задания правой части ||b-b*||∞ /||b||∞ =

=0.005/2.51≈0.2% привела к относительной погрешности решения

||x –x*||∞ / ||x||∞ =0.9364/1.981 ≈ 47.3%. Погрешность увеличилась в

237 раз. Это несколько меньше, чем cond(A) = 378. Однако, видимо существует такой вектор b, при котором коэффициент возрастания достигает значения cond(A) .

11.4.Прямые методы решения систем линейных алгебраических уравнений.

11.4.1.Метод Гаусса и LU разложение.

128

1.Схема единственного деления.

Из числа прямых наиболее широкое применение находит метод исключения Гаусса. Суть этого метода состоит в эквивалентном преобразовании к системе с верхней треугольной матрицей. При этом формируется последовательность (A | b)→(A1 | b1)→…→(Ak | bk)→…→

→(An-1 | bn-1)=(U | bn-1), где U верхняя треугольная матрица. Матрица

Ak должна содержать нули ниже главной диагонали в первых k столбцах.

Пусть уже осуществлено обнуление элементов ниже главной диагонали в предыдущих k-1 столбцах, т. е. сформирована матрица Ak-1 . Опишем аналогичную процедуру обработки k – ого столбца, переход

(Ak-1 | bk-1)→(Ak | bk). В простейшем варианте для обнуления элементов k-

ого столбца матрицы Ak-1 ниже главной диагонали осуществляют следующие операции: рассчитывают коэффициенты lk+i,k = ak+i,k/ak,k

(i=1,2,…,n). Здесь ai,j элементы матрицы Ak. Строку с k-ым номером матрицы (Ak-1 | bk-1) умножают на lk+i,k и результат умножения вычитают из строки с i-ым номером этой матрицы (i=1,2,…,n). Указанные операции должны быть последовательно произведены для всех столбцов с номерами k=1,2,…,n-1.

В результате (n-1) таких операций получаем систему с верхней треугольной матрицей, последнее уравнение которой содержит лишь одну неизвестную. Последовательной подстановкой снизу вверх (обратный ход)

находим остальные неизвестные.

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

неизвестных, чтобы обнулить n-1 элементов первого столбца ниже главной диагонали, нужно вычислить n-1 штук элементов l и для каждого из них провести (n-1) операций умножения и столько же вычитания, всего

Q1=2 n(n-1)+n-1=2(n-1)2+3(n-1) операций. На каждом следующем шаге

129

число обрабатываемых уравнений уменьшается на единицу. Поэтому

Q2=2 n(n-2)+n-2=2(n-2)2+3(n-2), Qk=2(n-k)2+3(n-k) .Общее число операций

	n 1	n 1	n 1
	Q Qk 2 (n k)2 3 (n k)
	k 1	k 1	k 1
	n 1	n 1	2(n 1)n(2n 1)		(n 1)n		2
	2 i2 3 i			3				n3.
			6		2		3
	i 1	i 1

Обратная подстановка идет значительно быстрее. Последнее неизвестное находится с помощью одной операции ( деление на последний ведущий элемент), предпоследнее уравнение требует двух операций

(умножение-вычитание, а потом деление). На k–ом шаге требуется k

n	1		n	2
операций, k=1, …, n. Общее число операций Q k		n(n 1)			.
k 1	2	2

Метод исключения Гаусса эквивалентен разложению матрицы на произведение треугольных матриц A=LU, где U верхняя треугольная матрица, а L нижняя треугольная с единицами на главной диагонали.

На первом шаге метода рассчитывают набор коэффициентов li1, i=2, …, n. Сформируем матрицу

		1	0	0
1		l21	1	0		. Легко убедится, что
L1
		l	0	1
		31

			1	0	0 a11		a12	a13			a11		a12			a13
	1		l21	1	0		a22	a23				0	a22 l21a12						A1
	L1	A				a21										a23 l21a13
			l	0	1	a	a	a				0	a l a			a l a
			31			31	32	33					32	31	12	33 31 13
		Для обнуления элемента ниже главной диагонали во втором
											1	0	0			a32 l31a12
								1
	столбце формируем матрицу L2										0	1	0	,l32			. Тогда
																a22 l21a12
											0	l	1
												32

130