Posobie_MathCAD_v2

x		a(1)	x		a(1)
	1	12		2	13
		a22(1) x2 a23(1)
...
		(1)			(1)
		am2	x2 am3

x3 ... a1(1m)x3 ... a2(1)m

x3 ... amm(1)

xm f1(1)xm f2(1)

xm fm(1)

Если a22(1) 0 , то, продолжая аналогичное исключение, приходим к системе уравнений с верхней треугольной матрицей

x

a(1)

x

a(1)

x

... a(1) x

f (1)

1

12

2

13

3

1m

m

1

x2 a23(2)

x3 ... a2(2m)

xm f2(2)

x3 ... a3(3)m xm f3(3)

xm fm(m)

Из нее в обратном порядке находим все значения xi:

x

m

f

(m)

m

(m 1)

(m 1)

xm 1

xm

fm 1

am 1m

(1)

(1)

(1)

x3

(1)

xm

x1

f1

a12

x2 a13

... a1m

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

a(0)		a(0)	, i, j = 1, 2, …, m,
ii		ij

31

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

Метод Гаусса – Жордано использует аналогичный прием с исключением элементов в столбцах матрицы таким образом,

чтобы перевести исходную систему			к			, где	A –
	Ax	f		A x	f

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

Пример 2.1. Решить СЛАУ

x1 x2 x3 22x1 x2 x3 3x1 x2 x3 6

методом Гаусса – Жордано и найти обратную матрицу. Записываем расширенную матрицу и с помощью эквивалентных преобразований переводим матрицу системы в единичную:

	A			I		f
1	1	– 1	1	0	0	2
– 2	1	1	0	1	0	3
1	1	1	0	0	1	6

Прибавляем ко второй строке первую, умноженную на 2, а из третьей вычитаем первую:

1	1	– 1	1	0	0	2
0	3	– 1	2	1	0	7
0	0	2	– 1	0	1	4

Делим вторую строку на 3 и вычитаем ее из первой:

32

1	0	1	3		1	3		1	3			0		1	3
0	1	13			2 3			13				0		7 3
0	0	2			– 1									4
								0				1

Делим последнюю строку на 2 и прибавляем ее с коэффициентом 13 к первой строке и с коэффициентом 13 – ко второй:

1	0		0			0				1	3			1	3			1
0	1		0			12				13				16				3
0	0		1			1	2			0				1	2			2
		I							A 1									x
				A A 1
Выполняя		проверку					I , Ax f						,	убеждаемся,				что

решение найдено правильно.

Часто возникает необходимость в решении СЛАУ, матрицы которых являются слабо заполненными, т.е. содержат много нулевых элементов. В то же время эти матрицы имеют определенную структуру. Среди таких систем выделим системы с матрицами ленточной структуры, в которых ненулевые элементы располагаются на главной диагонали и на нескольких побочных диагоналях. Для решения систем с ленточными матрицами коэффициентов вместо метода Гаусса можно использовать более эффективные методы, например, метод прогонки.

Рассмотрим наиболее простой случай: систему с трехдиагональной матрицей коэффициентов, к которой сводится решение ряда численных задач (сплайн-интерполяция таблично заданной функции, дискретизация краевых задач для дифференциальных уравнений методами конечных разностей и др.). В этом случае СЛАУ имеет вид:

33

c1 x1 b1 x2 f1

ai xi 1 ci xi bi xi 1 fi

am xm 1 cm xm

fm

Матрица системы (2.6) – (2.8) имеет трехдиагональную структуру, что хорошо видно из следующего эквивалентного векторноматричного представления:

c1		b1		0 0 0				0
a	2	c	2	b	0		0	0
				2
0		a3		c3	b3	0	0
0		0 a4			c4	b4	0
		0 0 0 am 1						cm 1
0
0		0		0	0		0	am

	0		x1					f1
	0		x					f
					2				2
	0
				x3				f3
	...									.
	...
	...
								fm 1
	bm 1		xm 1
			x					f
	c				m				m
		m

Если при этом выполняется условие ci bi ai , то говорят, что

матрица данной системы имеет диагональное преобладание. Предположим, что существуют такие наборы чисел αi

i = 2, 3, …, m, при которых

xi i 1xi 1 i 1 .

Записав уравнение (1.6) в виде (2.9):

и βi,

(2.9)

x1 b1 x2 f1 , c1 c1

получим формулы для определения α2, β2:

2		b1	,	2		f1	.	(2.10)
		c1
						c1

Уменьшим в (2.9) индекс на единицу: xi 1 i xi i , и подставим полученное выражение в (2.7):

ai i xi ai i ci xi bi xi 1 fi ,

откуда

34

x

bi

x

ai i fi

.

i

i 1

ci ai i

ci ai i

Данное равенство

совпадает

с (2.9), если при

всех

i = 1, 2, …, m – 1 выполняются рекуррентные соотношения

i 1

bi

, i 1

ai i fi

.

(2.11)

ci ai i

ci ai i

Они позволяют получить все остальные коэффициенты αi, βi. При i = m – 1 из (2.9) получим xm 1 mxm m . Подставляя

это выражение в (2.8) и разрешая полученное выражение относительно xm, записываем:

x	am m fm	,	(2.12)
m	cm am m

где αm и βm известны. Далее по формулам (2.9) последовательно находятся xm–1, xm–2, …, x1.

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

Будем называть прогонку корректной, если знаменатели прогоночных коэффициентов в формуле (2.11) не обращаются в нуль, и устойчивой, если |αi| < 1 при i = 1, 2, 3, …, m.

Теорема. Пусть коэффициенты ai, bi уравнения (2.7) отличны от нуля и пусть ci bi ai при i = 1, 2, 3, …, m. Тогда

прогонка корректна и устойчива.

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

2.6. Итерационные методы решения СЛАУ

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

35

векторов
	x	(0)	, x	(1)	, x	(2) , , x(k ) , , пределом которой должно
									(k ) .
быть точное решение							x	* lim x
								k

На практике построение последовательности обрывается, как только достигается желаемая точность. Чаще всего для достаточно малого значения 0 контролируется выполнение оцен-

ки	*	(k )	.
	x	x

Метод последовательных приближений может быть построен по следующей схеме. Эквивалентными преобразованиями при-

ведем систему (2.5) к виду
				(2.13)
x	C x		d ,
				– некоторые новые мат-
где x – тот же самый вектор, а		C и d

рица и вектор соответственно.

При решении методом последовательных приближений не-

обходимо выбрать начальное (нулевое) приближение. За нуле-

вое приближение можно принять столбцы правых частей f , d

системы (2.13) или нулевой вектор. Следующее приближение

1	определяется рекуррентным равенством
x
			1	0
	2		x	C x	d .
			:
Далее находим x
			2	1
			x	C x	d ,
и т.д. Для k-й итерации получаем
			k 1	k		0, 1, 2, ... (2.14)
			x	C x	d , k

Такой итерационный процесс будем называть одношаговым итерационным методом.

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

C				k при k имела пределом
	, чтобы последовательность x
*	, т.е. была решением системы (2.13), эквивалентной исходной
x
	k	*	.
системе (2.5): lim x		x
	k

36

Достаточным условием сходимости итерационного метода (2.13) к решению системы (2.5) при любом начальном векто-

0		является требование		C		1, где		C		– норма матрицы С.
ре x

В общем виде одношаговый итерационный метод можно записать в так называемой канонической форме:

			(k 1)		(k )	(k )
	(k 1) x			x
B						Ax	f .
			(k 1)

Здесь B(k 1) – матрица, задающая итерационный метод, k 1 –

итерационный параметр. Если B(k 1) E (где E – единичная матрица), то метод называют явным, а в противном случае –

неявным. Если матрица B(k 1) и итерационный параметр(k 1) не зависят от номера итерации ( B(k 1) B, (k 1) ), то

метод называют стационарным, и нестационарным – в про-

тивном случае.

Использование неявных методов сопровождается обращени-

ем матрицы Bk 1 , поэтому для сохранения эффективности алгоритма эта матрица должна быть легко обратима (например,

Bk 1 – диагональная, треугольная, трехдиагональная или ортогональная).

Приведем некоторые примеры. Для этого представим мат-

рицу	A в виде A A1 D A2 , где D – диагональная матрица,
A1 –	левая треугольная, A2 – правая треугольная. A1 и A2 име-

ют нулевую главную диагональ. Ненулевые элементы всех трех матриц совпадают с соответствующими элементами матрицы A.

Метод релаксации (простой итерации). Здесь B(k 1) E , а

(k 1) .

Итерационный метод Ричардсона. Здесь B(k 1) E , а (k 1)

–переменный параметр.

Метод Якоби. B(k 1) D , а (k 1) 1.

37

Метод верхней релаксации.	B(k 1) D A ,	k 1 ,
	1

0 2 , где – заданный числовой параметр. Для 1 как частный случай получается метод Гаусса –Зейделя. Для симметричных положительно определенных матриц A условие 0 2 является условием сходимости метода.

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

Предположим, что диагональные элементы матрицы A ис-

ходной системы (2.5) не равны 0 (aii 0, i = 1, 2, …, m). Разрешим первое уравнение системы (2.5) относительно x1, второе

относительно x2 и т.д. Получим следующую эквивалентную систему, записанную в скалярном виде:

	x1 C12x2 C13x3 ... C1n xn d1
	x2	C21x1	C23x3	... C2n xn d2	,
	....................................................
	xn	Cn1x1	Cn2 x2	... Cn,n 1xn 1 dn

что совпадает с формулой (2.13), в которой матрица C и век-

тор d определены по следующим формулам:

		aij	, i j
						fi
Cij		a		,	di		, i 1, 2, ..., n.	(2.15)
		ii	i j			aii
	0,

Метод, основанный на таком приведении системы (2.5) к виду (2.13), называют методом Якоби. Теперь, задав нулевое приближение, по рекуррентным соотношениям (2.14) можем выполнять итерационный процесс.

Условие сходимости C 1 в методе Якоби равносильно условию диагонального преобладания для исходной матрицы A:

			m
	aii			aij	,	i 1, 2, ..., m.
			j 1
			i j

38

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

венства на aii , получим неравенство

m				a
1				ij				,	i 1, 2, ..., m.
				a
j 1
				ii
i j

С учетом (2.15) можно перейти к следующему неравенству:

m
1	Cij	,	i 1, 2, ..., m.
j 1

То есть сумма элементов любой строки матрицы C меньше 1. Выполнение такого условия равносильно выполнению условия

				m
сходимости метода	C		max		Cij	1 .
			1 i m
				j 1
Под методом Гаусса – Зейделя обычно понимается такое ви-

доизменение одношагового итерационного метода (2.15) решения СЛАУ, в котором для подсчета i-й компоненты (k + 1)-го

приближения	x	k 1	*	используются уже
			к искомому вектору x
	i

вычисленные на этом, т.е. (k + 1)-м шаге, значения первых i – 1 компонент xi 1 k 1 . Это означает, что если система (2.5) тем или

иным способом сведена (например, с помощью метода Якоби) к системе (2.13) с матрицей коэффициентов С и вектором свобод-

ных членов d ,

то ее приближение к решению по методу Зейде-

ля определяется системой равенств

x

k 1 C x k C x

2

k

C x

k ...

C

x

k

d

1

1

11 1

12

13

3

1m

m

x

k 1 C

x k 1 C

22

x

k C

x k

... C

2m

x

m

k d

2

2

21 1

2

23 3

.

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

x

m

k 1 C

x

k 1 C

m2

x

k 1

... C

m,m 1

x

k 1

C

m,m

x

m

k d

3

m1 1

2

m 1

С точки зрения компьютерной реализации одношагового

итерационного метода использование метода Гаусса – Зейделя

означает, что элементы массива x будут постепенно замещать-

39

ся новыми элементами. В связи с такой интерпретацией метод Гаусса – Зейделя иногда называют методом последовательных смещений.

Иногда исходную систему (2.5) не удается привести к виду (2.13), выполнив при этом условие сходимости метода. В этом случае можно воспользоваться методом релаксации. Этот метод основывается на соотношении

k 1

k

k

x

x

,

Ax

f

откуда

k

, где τ

– итерационный параметр.

x

k 1 x

k Ax

f

Скалярные формулы метода релаксации имеют следующий вид:

x

k 1 x

k a x k a x

2

k ...

a

x

n

k f

1

1

1

11 1

12

1n

x

k 1 x

k a

x

k a

x

k ...

a

x

k f

. (2.16)

2

2

21 1

22

2

2n

n

2

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

x

n

k 1 x

k a

x k a

n2

x

k ...

a

nn

x

k f

n

1

n1 1

2

n

Раскрыв скобки, можно привести (1.16) к виду (1.14), где коэф-
фициенты	матрицы		C и			вектор		свободных членов d		будут
	Cij	1 aij				,i j		di fi ,	i 1, 2, ..., n.
иметь вид:					,	i j	,			Подбо-
		a		ij

ром параметра τ можно добиться сходимости метода релаксации. При использовании итерационных методов мы можем найти

решение исходной СЛАУ (2.5) лишь приближенно с заданной точностью. Поэтому важной проблемой является вопрос о способе остановки итерационного процесса при достижении точности. Наиболее простой способ – это сравнение между собой соответствующих неизвестных с двух соседних итераций: (k + 1) и (k). Если максимальная из всех разностей становится меньше заданной точности , то итерационный процесс останавливается:

max xik xik 1 .

1 i m

2.7. Решение СЛАУ в пакете MathCAD

40