book1989

где ф определено согласно (30). Полагая ct= —0,5/, получим

			yW = sin nkj		k, j — 1,2, .... N — 1.	(32)
Собственные функции				(32)	определены с точностью до произволь­
ного	постоянного (не			за­
висящего от /) множителя.
Интересно сопоставить
решения дифференциаль­
ной	(24)	и разностной
(25)	задач	на	собствен­
ные	значения.		Значения
собственных функций (32)
разностной задачи совпа­
дают в точках сетки со
значениями		собственных			Рис. 3. Собственные значения дифференциаль­
функций дифференциаль­					ной задачи (сплошная черта) и разностной
ной задачи. Спектр диф­					схемы
ференциальной			задачи не		в то время как спектр разностной
ограничен, т. е. %к-уоо при

задачи ограничен сверху при каждом фиксированном шаге h числом 4/г 2. Для каждого фиксированного номера k ^ k a, где ka не зависит

от h, собственные значения ?4Л) разностной задачи сходятся при h-+0 к соответствующему собственному значению Хк дифференци­ альной задачи, т. е.

	n k h	= Яь
lim — sin2
й-»о №	2 (ft — а)
При этом собственные	значения	разностной задачи (25) всегда

меньше соответствующих собственных значений дифференциаль­

ной задачи (24). Погрешность l h—'k<k) минимальна для малых но­ меров k и сильно возрастает с ростом k. На рис. 3 изображены гра­

фики %к (сплошная черта) и № в зависимости от номера k для значений а = 0, b= 1, N = 25 и N = 50.

5. Свойства собственных значений и собственных функций.

Перечислим свойства собственных значений и собственных функ­

ций разностной задачи (25).				Прежде всего из (31) видно, что
0 <		< ... <		4 г) < A $i< ... <		MU < 4 - .
						h*
Последнее неравенство неулучшаемо, так как
			=	— cos^	n h
		A / V - i			2 ( 6 - а)
n h	.	при	/i->-0. Оценку снизу			для наименьшего
И COS2 ------------ >■1

2 ( b — а)

41

собственного	значения Я, можно уточнить. Обозначая а =
= nh/(2(b—а)), получим
	л (М_л
	—•Aj
где	2
	— наименьшее собственное значение дифферен­

циальной задачи. Не ограничивая общности, можно предположить, что h ^ .(b —а)/3. Тогда а ^ я /б , и поскольку функция sinct/ct мо­ нотонно убывает при а е [0 , я/6], получим

sin «	\ 2 ^	М	6 \ 2 ___9_
а	/ ^	\ 2	л j	л2	'
;.[Л) 5>9/(Ь — а)2.					(33)

Таким образом, наименьшее собственное значение задачи (25) отделено от нуля константой = 9/(b—а)2, не зависящей от h.

Покажем теперь, что собственные функции (32) задачи (25), отвечающие различным собственным значениям, ортогональны в смысле скалярного произведения

w-i
(и, v) = ^ uivih-	(34)
/=i

Запишем уравнение	(25) для функций у т и у{1) в виде
	„(*>	I	__п	(35)
	Ухх.,+ХьУ1 = U>
	„(О		(Л)„(О _п	(36)
	Ухх,,- +		у/ =
Умножим уравнение	(35)	скалярно на у {1), уравнение		(36) —

на у(к) и вычтем из первого полученного равенства второе. Тогда будем иметь

(yfx, у{1)) - (У?х, У{к)) = Q-ih)- W ) {У{к\ У{1))•

(37)

Из разностного аналога формулы интегрирования по частям (23), учитывая условия yW = yW = 0, получим

и точно так же

,(0 „(«1	■ f t M l
в & л = - < е » п

Следовательно, левая часть равенства (37) обращается в нуль, и поскольку ЦМфХМ при кф1, получаем

(у т , у (1)) = 0, если кф1.

Множество функций

У= (Уо, Уи Уг , - . . , Ух-и Ух), yj = y(Xj),

заданных на сетке (2) и удовлетворяющих нулевым граничным

42

условиям Уа = Уп— 0, образует (N—1)-мерное линейное пространст­

во Н относительно	покоординатного сложения и умножения на
число. Собственные	функции уш, k= \, 2,	N—1,	задачи (25)
ортогональны и, следовательно, линейно независимы			в Н. Тем са­
мым множество собственных функций задачи		(25) образует орто­
гональный базис в Н. Нетрудно показать, что

\yM f = 2 h ( y f) * = 0 ,b (b -a )

/=i

для всех k= \, 2, ..., N— 1. Следовательно, множество собственных функций p(ft), k = 1, 2 ,... , с координатами

P f =	У —— sin ~ , /'= 1 ,2 , .... N — 1,
1	1 b —a N

образует в Н ортонормированный базис. Любой элемент у ^ Н можно единственным образом представить в виде разложения

N- 1

у = 2 Ск^ к)-

k = \

6.Разрешимость и сходимость разностной задачи. Обратимся

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

yi-i—2y{+yi+1= 0, i = l , 2, ..., N— 1, Уо~У;я=0

имеет только тривиальное решение и, следовательно, неоднород­ ная система (16) имеет единственное решение.

Исследуем сходимость при h-*0 решения разностной задачи (16) к решению исходной дифференциальной задачи (8)— (9). Обозначим через

2j= t/i—tl[xi) , Х;£=(Ол,

погрешность в точке хи т. е. разность между решениями задач (25) и (8)— (9). Подставляя в (16) вместо у{ сумму г{+и(х{), i= = 1, 2, ..., N—1, получим, что погрешность удовлетворяет разно­ стному уравнению

"г<~1~ « + г<+1	.... Л - 1 ,	20 = гл =0,(38)
И?
г д е
ф; =	U~x i + /ь	(39)

Сеточная функция ф; называется погрешностью аппроксимации или невязкой разностной схемы (16) на решении задачи (8)— (9).

43

Записывая	в виде
	'А = ( u x x , i — и '' ( * ‘))			+( и '' (х	‘) +	А )
и учитывая, что согласно		(8)	выполняется		равенство u"(xi)+fi= О,
получим	b	=	u7xtl — u"(xi).
Разложение по формуле Тейлора показывает, что если uIV(*)
ограничена,	то ф,= 0(/12)	при	0. По этой причине говорят,				что
разностная схема (16) имеет второй порядок						аппроксимации	на
решении исходной задачи		(8) — (9). Наша ближайшая цель — по­
казать, что схема (16) сходится, т. е.				0 при h-*-0 и, более того,
имеет второй порядок точности, т. е. z,==0(/i2).							(38)
Воспользуемся возможностью получить					решение задачи
в явном виде. Уравнение		(38) отличается			от	изученного ранее

уравнения (18) только обозначениями. Поэтому согласно (20) ре­

шение задачи	(38) представляется в виде
* =	2	h 2	h^ i ~ 2	h 2	h4t'	i = 2 ’ 3........ (4°)
	*=i	/=i	k = i	/=i
		*i = h (b — a)-12			h 2	Лф/.

*=i /=i

Из разложения по формуле Тейлора

+ ( J K ,v &)

иограниченности Hiv(*) следует, что существует постоянная Mt>

>0 , не зависящая от А и от / и такая, что

/= 1, 2, ..., У -1 . Поэтому из формулы (40) следует оценка

|zi|< (M 1A*) \	№(Л/~ 1} N + A2	(t~~-'-А 1 ,
L й — а	2	2	J
Т. е.
I* i:< (Mth2) [ j (N - 1)N + t (t-		1)] ^ .
Выражение в квадратных скобках равно i(N+i—2) и оценива­
ется сверху числом 2АД Таким образом,
|z,| (Mth2)N2h2 = Mt (b—a)2h2,
т. e. | 2(| = 0 (h z) при A->-0, i= l,	2........У—1. В этом случае говорят,

что схема имеет второй порядок точности.

Отметим, что приведенный здесь способ исследования сходи­ мости, основанный на явном представлении решения разностной задачи, непригоден для более сложных задач. Другие методы ис­ следования сходимости разностных схем излагаются в части III.

44

7.	Метод прогонки. Система уравнений	(16) представляет со­
бой частный случай систем линейных алгебраических уравнений
	Ay = f
с трехдиагональной матрицей А= [ай], т. е. с		матрицей, все эле­

менты которой, не лежащие на главной и двух побочных диагона­ лях, равны нулю (а^= 0 при />Н -1 и /C t —1).

Вобщем случае системы линейных алгебраических уравнений

стрехдиагональной матрицей имеют вид

с}у}+Ьдш = —/,, /= 1, 2, .... Л —1,		(41)
=	Ук=	(42)

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

Приведем вывод расчетных формул метода прогонки. Будем искать решение системы (41) в виде

	t/j=aJ+iyJ+I+Pj+1,			/ —0,	1,		..., N	1,	(43)
где a i+i, Ря-1 — неизвестные			пока	коэффициенты. Отсюда найдем
Уt-i = > /* //1 +	Р/ =Щ(a /+i*//+i +		Р/+0 +		Р/	=	Р/),	/ =	.1,. .2,,N — 1.
	=	a / a /+i*//+i +		(а/Р/+1		+
Подставляя полученные выражения для yh								yj^i	в уравнение
(41), приходим при /= 1, 2 ,..., N—1 к уравнению
[ C L j+ i	О )	^ j ] УJ+ I~П [ Pj+i ( £ijCXj					Cj) ^	' / j ]	0 .

Последнее уравнение будет выполнено, если коэффициенты aj+i, pm выбрать такими, чтобы выражения в квадратных скоб­ ках обращались в нуль. А именно, достаточно положить

*/-1-1 ■	й/	Р/э	aj$i+ fj	/= 1 ,2 ,	, N — 1. (44)

Соотношения (44) представляют собой нелинейные разностные уравнения первого порядка. Для их решения необходимо задать начальные значения alt р,. Эти начальные значения находим из требования эквивалентности условия (43) при / = 0, т. е. условия ya= alyl+ $i, первому из уравнений (42). Таким образом, получаем

CCl = « i. Pl = H-f

(45)

Нахождение коэффициентов aj+u pj+1 по формулам (44), (45) называется прямой прогонкой. После того как прогоночные коэф­ фициенты а,-+1, Pi+i, / = 0, 1, ..., N—1, найдены, решение системы (41), (42) находится по рекуррентной формуле (43), начиная с

45

}=N—1. Для начала счета по этой формуле требуется знать yN, которое определяется из уравнений

Уя—X2y.v_i+ p2, Ук-i —

и равно (х2рjv+р,2)/(1—х2«№).	Нахождение у, по формулам
У/ = сty+i^+i + Р/+1,	j = N ~ l ,	N — 2, . . О,
=	ИгРл/ "Ь Р2	(46)
	—:---------
	1— Щ<У-м

называется обратной прогонкой. Алгоритм решения системы (41), (42), определяемый по формулам (44) — (46), называется методом прогонки. Применяются и другие варианты метода прогонки (см. [32]).

Метод прогонки можно применять, если знаменатели выраже­ ний (44), (46) не обращаются в нуль. Покажем, что для возмож­ ности применения метода прогонки достаточно потребовать, что­ бы коэффициенты системы (41), (42) удовлетворяли условиям

а,фО, Ь3ф 0,	|Cj| 5г |aj| + |bj|, /= 1,	2, .... N—1,	(47)
	M s S l, М < 1 .		(48)
Заметим, что числа а}, Ьи си хь х2 могут быть комплексными.			(48)
Сначала докажем	по индукции, что при	условиях (4 7 ),

модули прогоночных коэффициентов <Zj, j= 1, ..., N—1, не превос­

ходят единицы. Согласно (4 5 ), (48)				имеем	| oti | = |	^ 1. Предпо­
ложим, что	|а ,|^ 1	для некоторого j		и докажем, что		|a j+1| ^ l . Из
оценок
	|С—а д | > | |сл-1 —| ctj| | а,| | > | |Cj| —Jал-| |
и условий (47) получаем
		| c—	afii| > | bj| > 0 ,
т. e. знаменатели выражений (44) не обращаются						в нуль. Более
того,
Следовательно,		|сс5|^ 1 ,	/= 1, 2, ..., N. Далее, учитывая второе
из условий	(48) и	только	что доказанное		неравенство |а Л-]<;1,
имеем	]1—x2a jv |^ l — |х 2| | оьлг| ^ 1 —|х 2| >0,

т.е. не обращается в нуль и знаменатель в выражении для уЛ-

Каналогичному выводу можно прийти и в том случае, когда условия (47), (48) заменяются условиями

Щ-ф0, Ь,ф0,	| C j | > | f l i | +	|6j|,	у =	1, 2,	N—1,	(49)
	Ы < 1 >	|х * |< 1 .				(50)
В этом случае из предположения \а3\ ^ 1			следует
I с—ajfl, | ^	11 с, |— |а, 11 > |bj |,			|а,+1 1<	1,
46

т. е. все прогоночные коэффициенты, начиная со второго, по моду­ лю строго меньше единицы. При этом |1—y.2aN| ^ 1— |x 2||a w| ^

1 —|a w| >0.

Таким образом при выполнении условий (47), (48) (так же как

иусловий (49), (50)) система (41)—(42) эквивалентна системе

(44)—(46). Поэтому условия (47), (48) (или условия (49), (50)) гарантируют существование и единственность решения системы

(41), (42) и возможность нахождения этого решения методом про­ гонки. Кроме того, доказанные неравенства | щ | ^ 1, j = 1, 2, ..., N, обеспечивают устойчивость счета по рекуррентным формулам (46). Последнее означает, что погрешность, внесенная на каком-либо шаге вычислений, не будет возрастать при переходе к следующим шагам. Действительно, пусть в формуле (46) при / = /0+1 вместо yk+i вычислена величина yh+i=yk+i4-6ji)+1. Тогда на следующем шаге вычислений, т. е. при / = /„, вместо у^--=ак+1ук+1+^к+1 получим ве­

личину	yJrj = ajo+1 (i/j0-H + 6j0+i) + Pj0+I и погрешность	окажется равной
	®/«= у/а уit, — T.fA'iH'
Отсюда	получим, что |6jo[ ^ |a jo+,| 16jo+i | sg 16J0+11, т. е. погрешность
не возрастает.		(16), записанной
Отметим, что для разностной краевой задачи
в виде	ys-i—2«/;+ «/гц= —h2fjt /= 1, 2, ...,	N—1,

имеем ctj= bj= 1, с, = 2, xi = x2 = 0. Поэтому выполнены условия устой­ чивости (47), (48) и решение задачи (16) можно отыскивать ме­ тодом прогонки.

Ч А С Т Ь II

ЧИСЛЕННЫЕ МЕТОДЫ АЛГЕБРЫ И АНАЛИЗА

Г Л А В А 1

ПРЯМЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ

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

Ax = f,	(1)
где А — матрица тУ.т, х= {хи х2, ...,	хт)г— искомый вектор, / =

= (f,, f2, . . . , fm)T— заданный вектор. Предполагается, что опреде­ литель матрицы А отличен от нуля, так что решение х существует и единственно. Для большинства вычислительных задач характер­ ным является большой порядок матрицы А. Из курса алгебры из­ вестно, что систему (1) можно решить по крайней мере двумя спо­ собами: либо по формулам Крамера, либо методом последователь­ ного исключения неизвестных (методом Гаусса). При больших т первый способ, основанный на вычислении определителей, требует порядка т\ арифметических действий, в то время как метод Гаус­ са— только 0{т 3) действий. Поэтому метод Гаусса в различных вариантах широко используется при решении па ЭВМ. задач линей­ ной алгебры.

Методы численного решения системы (1) делятся на две груп­ пы: прямые методы и итерационные методы. В пряШых (или точ­ ных) методах решение х системы (1) находится за конечное число арифметических действий. Примером прямого метода является ме­ тод Гаусса. Отметим, что вследствие погрешностей округления при решении задач на ЭВМ прямые методы на самом деле не приводят к точному решению системы (1) и называть их точными можно лишь отвлекаясь от погрешностей округления. Сопоставление раз­ личных прямых методов проводится обычно по числу арифметиче­ ских действий (а еще чаще— по асимптотике при больших m чис­ ла арифметических действий), необходимых для получения реше­ ния. При прочих равных условиях предпочтение отдается методу с меньшим числом действий.

Итерационные методы (их называют такжметодами последо­

вательных приближений) состоят в том, что решение х системы (1) находится как предел при п->-оо последовательных приближений х1п), где п — номер итерации. Как правило, за конечное число ите­ раций этот предел не достигается. Обычно задается некоторое ма­

48

лое число е>0 (точность) и вычисления проводятся до тех пор, пока не будет выполнена оценка

\\х(п)—х!|<е. (2)

Число итераций я = я(е), которое необходимо провести для по­ лучения заданной точности е (т. е. для выполнения оценки (2)), для многих методов можно найти из теоретических рассмотрений. Качество различных итерационных процессов можно сравнивать по необходимому числу итераций п(е).

К решению систем линейных алгебраических уравнений сводит­ ся подавляющее большинство задач вычислительной математики. В настоящее время предложено колоссальное количество алгорит­ мов решения задач линейной алгебры (см. [8, 35]), большинство из которых рассчитано на матрицы А специального вида (трехдиа­ гональные, симметричные, ленточные, большие разреженные мат­ рицы) .

Прямые методы, которые рассматриваются в гл. 1, не предпо­ лагают, что матрица А имеет какой-либо специальный вид. На прак­ тике они применяются для матриц умеренного порядка (порядка ста). Итерационные методы, рассмотренные в гл. 2, можно приме­ нять и для матриц высокого порядка, однако их сходимость не очень быстрая. Более совершенные прямые и итерационные мето­ ды, учитывающие структуру матрицы, излагаются в части III.

§ 1. Метод Гаусса численного решения систем линейных алгебраических уравнений

1. Основная идея метода. В ближайших двух главах рассмат­ риваются численные методы решения системы линейных алгебраи­ ческих уравнений

Ax = f,

(1)

где А — вещественная квадратная матрица порядка т, a f — за­ данный и х — искомый векторы. Будем предполагать, что опреде­ литель матрицы А отличен от нуля. Тогда для каждого вектора f система (1) имеет единственное решение.

Запишем систему (1) в развернутом виде

а 1\х 1 + ^ 12*2 +	. . . +	OimXni	=	f l ,
@2lXl ”1“ @22Х2 "	• • •	Щ.-ЛД;	=	/2,	(2)
@mixi “Б ®тгх2		@тт т—[:п■

Метод Гаусса решения системы (2) состоит в последовательном исключении неизвестных хи хг, ..., хт из этой системы. Предполо­ жим, что апфО. Поделив первое уравнение на ап, получим

Xi + cizxz-1---- bCimXm= i/,,	(3)
где
Cii = — , / = 2, ... , m, i/i =	— .
« И	«11
	4 9

Рассмотрим теперь оставшиеся уравнения системы (2):

ciiiXi + a,-2хг+ ...

+a{mxm=fi, г = 2, 3, ..., т.

(4)

Умножим (3) на а,-, и вычтем полученное уравнение из i-ro урав­ нения системы (4), i = 2, т. В результате получим следующую ^систему уравнений:

Х 1 +	с \ г х г +	■. .	+	С ц -Х / + . . .	+	С1тХ т =	У и
	+	• • •	+	<*№ + ...	+	а‘% =	Z'1',
							(5)
аых		“Ь йт/Х/ +					и»
«-т2Л2+						О'тт.Хп = /!
Здесь обозначено
aij —Q-Ц	Ci/&iu	/I1'		У 1^4у			(6)

Матрица системы (5) имеет вид

'1	С12				С1 т
0	{1)				а(1)
	а22				w2т
0	а{1)	•	.	*	а(1)
	u/n2				ц/лт

.Матрицы такой структуры принято обозначать так:

1	X ....	X
0	X ....	X

0 X ..•. X

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

flM*2 + . . . + d^Xj + . . . + а^тХт. =

(7)

ail)x	+ Omi*	.fW
ит2л2+		~f" Q/nniXm ’ - m •

Тем самым мы осуществили первый шаг метода Гаусса. Если ,аМфО, то из системы (7) совершенно аналогично можно исклю­

чить неизвестное х2 и прийти к системе, эквивалентной (2) и имею­ щей матрицу следующей структуры:

1	X	X	...	X "
0	1	X	...	X
0	0	X	. . .	X
0	0	X	.. .	х _

При этом первое уравнение системы (5) остается без изменения.

50