Mathcad - ЛР2

Министерство образования Российской Федерации Санкт - Петербургский государственный политехнический университет

Институт менеджмента и информационных технологий

Кафедра ПО ВТ и АС

Лабораторный практикум по курсу вычислительной математики

ОТЧЕТ

по лабораторной работе №2

Тема:

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

Выполнил: _______________

Группа: _____

Вариант: 14

Проверил: Царев В.А.

Отметка о зачете __________________

" ___ " ______________ 2006 г.

Череповец

2006

ORIGIN := 1

TOL := 10− 15

Введем функции вычисления норм матриц(векторов):

Norm1(A) :=

Norm2(A) :=

NormInf(A) :=

for j 1 .. cols(A)

AbsSumInColumnj ←

max(AbsSumInColumn)

cols(A) rows(A)		(Ai , j)2
å	å
j = 1	i = 1
for	i 1 .. rows(A)

rows(A)

åAi, j

i= 1

cols(A)
AbsSumInRowi ← å		Ai, j

j = 1

max(AbsSumInRow)

Задача 2.1. Написать программу для решения СЛАУ, используя простейшую схему Гаусса. Решить

заданную систему линейных уравнений AХ=В при помощи составленной программы. Выполнить проверку, найти вектор невязки и его норму. Проанализировать полученные результаты.

Исходные данные:

N := 14	æ -0.25		0.67	1.00		0.36 ö			æ 3.45		ö
	ç	6.96×N	0.47	-0.11		0.55		÷	ç	0.02	÷
	ç							÷	ç		÷
A :=	ç	0.42	1	0.35		0.17		÷	B := ç	10.86 ÷
	ç							÷	ç		÷
	è	0.54	-0.32 -0.74			1.00 ø			è-5.92 ø
	æ-0.25		0.67	1	0.36		ö		æ	3.45	ö
	ç		0.47	-0.11	0.55		÷		ç	0.02	÷
	ç97.44						÷		ç		÷
A =	ç	0.42	1	0.35	0.17		÷		B = ç10.86		÷
	ç						÷		ç		÷
	è	0.54	-0.32	-0.74		1	ø		è-5.92		ø
GaussForward(A , B) :=			AB ¬ augment(A , B)

m¬ rows(AB)

n¬ cols(AB)

for i Î 1 .. m - 1 if ABi , i = 0

pos ¬ i + 1

while (ABpos, i = 0) Ù (pos < m + 1) pos ¬ pos + 1

return 0 if pos = (m + 1) for k Î i .. n

temp ¬ ABi , k

ABi, k ¬ ABpos , k

ABpos, k ¬ temp divisor ¬ ABi , i

for j Î i .. n

ABi , j ABi , j ¬ divisor

for i1 Î i + 1 .. m mul ¬ ABi1 , i for j Î i .. n

ABi1 , j ¬ ABi1, j - ABi , j ×mul

ABm , n

ABm , n ¬ ABm , n−1

ABm , n−1 ¬ 1

AB

GaussReverse(AB) :=		m ¬ rows(AB)
		n ¬ cols(AB)
		X ¬ ABánñ
		for i Î m - 1 .. 1
			n−1
		Xi ¬ Xi - å ABi , k×Xk
			k = i+1
		X
Решение задачи:
				æ-0.040998002313417		ö
				ç	12.8916089164482	÷
				ç		÷
X := GaussReverse(GaussForward(A , B))				X = ç	-3.60037178886858	÷
				ç	-4.4368213492501	÷
				è		ø
Вычисляем вектор невязки
V := A×X - B
æ	0		ö
ç	-7.4032446839567 ´ 10− 13		÷
ç			÷
V = ç		− 15	÷
ç			÷
ç	-1.77635683940025	´ 10	÷
è	0		ø

Norm1(V) = 7.4210082523507 ´ 10− 13

Norm2(V) = 7.40326599514668 ´ 10− 13

NormInf(V) = 7.4032446839567 ´ 10− 13

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

Задача 2.2. Написать программу для решения СЛАУ, используя схему Гаусса с частичным выбором

ведущего элемента по столбцу. Решить систему линейных уравнений из задачи 2.1. AХ=В при помощи составленной программы. Выполнить проверку, найти вектор невязки и его норму. Проанализировать полученные результаты.

GaussForwardPlus(A , B) := AB ¬ augment(A , B)

m¬ rows(AB)

n¬ cols(AB)

for i Î 1 .. m - 1
	maxpos ¬ i
	for	hÎ i + 1 .. m
	maxpos ¬ h if				ABh, i	>	ABmaxpos, i
	for	k Î i .. n				if maxpos ¹ i
		temp ¬ ABi , k
		ABi , k ¬ ABmaxpos, k
		ABmaxpos , k ¬ temp
	divisor ¬ ABi , i
	for	j Î i .. n
	ABi, j ¬		ABi, j
			divisor
	for	i1 Î i + 1 .. m
		mul ¬ ABi1 , i
		for j Î i .. n
		ABi1, j ¬ ABi1, j - ABi , j×mul

ABm , n

ABm , n ¬ ABm , n−1

ABm , n−1 ¬ 1

AB

			æ-0.040998002313409 ö
			ç	12.8916089164481		÷
			ç			÷
X2 := GaussReverse(GaussForwardPlus(A , B))			X2 = ç	-3.60037178886858		÷
			ç			÷
			è	-4.4368213492501		ø
Вычисляем вектор невязки:
	æ		0		ö
	ç		0		÷
V2 := A×X2 - B	V2 = ç				÷
				− 15
	ç				÷
	ç	-1.77635683940025 ´ 10			÷
	è		0		ø

Сравниваем векторы невязки результатов расчета обычным и расширенным методами Гаусса:

	æ	0	ö	æ	0	ö
	ç	− 13	÷	ç	0	÷
				V2 = ç		÷
	V = ç	-7.4032446839567 ´ 10	÷
	ç-1.77635683940025 ´ 10− 15		÷	ç	-1.77635683940025 ´ 10− 15 ÷
	ç		÷	ç		÷
				è	0	ø
	è	0	ø

Находим нормы векторов невязки:

Norm1(V) = 7.4210082523507 ´ 10− 13	Norm1(V2) = 2.61943244872498 ´ 10− 15
Norm2(V) = 7.40326599514668 ´ 10− 13	Norm2(V2) = 1.87399279458462 ´ 10− 15
NormInf(V) = 7.4032446839567 ´ 10− 13	NormInf(V2) = 1.77635683940025 ´ 10− 15

Вывод: Очевидно, что существует ощутимая разница в точности использованных

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

Задача 2.3. Задана СЛАУ AХ=F с трехдиагональной матрицей системы А. Написать программу для решения подобных систем методом прогонки и решить предложенную систему. Выполнить проверку, найти вектор невязки и его норму. Проанализировать полученные результаты. Компактное хранение элементов матрицы А в ЭВМ организовать с использованием одномерных массивов.

Формируем "трехдиагональную матрицу" - три массива коэффициентов a,b,c:

N := 14				n := 100 + N		n = 114
a :=		for		i Î 2 .. n			b :=					for	m Î 1 .. n		c :=			for		k Î 1 .. n - 1
				æ i		ö									2
												bm ¬ m + N							ck ¬ k + N ×sin(k)
			ai ¬ ln(i + N) ×cosç			÷
				è N ø								b						c
		a
				1									1							1
	1					0						1	197			1				3.25900311047774
	2			2.74434508303767								2	198			2				3.63718970730273
	3			2.76841360701211								3	199			3				0.581852699118848
	4			2.77319776646053								4	200			4				-3.21084105870682
	5			2.75864340775455								5	201			5				-4.17985400776466
	6			2.72479930025542								6	202			6				-1.24958409587954
	7			2.67181980063004								7	203			7				3.01069081920043
a =	8			2.59996610742638			b =					8	204		c =	8				4.64050151207359
	9			2.50960634343452								9	205			9				1.9764508228623
	10			2.401214631945								10	206			10				-2.66514826196204
	11			2.27536928734133								11	207			11				-4.99995103275352
	12			2.13275021102332								12	208			12				-2.73599577934964
	13			1.97413556464829								13	209			13				2.18325196634822
	14			1.80039778047188								14	210			14				5.24180142015986
	15			1.61249896086119								15	211			15				3.50190719132161
	16			1.41148571433883								16	212			16				-1.57691140918009

Формируем элементы столбца свободных членов:

f :=		for m Î 1 .. n		1
		fm ¬ N - m	1	13
			2	12
		f	3	11
			4	10
			5	9
			6	8
			7	7
		f =	8	6
			9	5
			10	4
			11	3
			12	2
			13	1
			14	0
			15	-1
			16	-2

Решение системы:

	S(a, b, c , f , n) := p ¬ b
	q ¬ b
	p	¬	c1
	1		b1
	q1	¬	f1
			b1

for i Î 2 .. n - 1

ci

pi ¬ bi - pi−1 ai

fi - qi−1×ai qi ¬ bi - pi−1×ai

fn - qn−1×an qn ¬ bn - pn−1×an

xn ¬ qn

for i Î n - 1 .. 1 xi ¬ qi - pi ×xi+1

x

x := S(a , b, c , f , n)

		1
	1	0.065018645601875
	2	0.058707159810799
	3	0.054313569815126
	4	0.04996781278962
	5	0.04490562621678
	6	0.039205734607081
	7	0.033544928797916
x =	8	0.028441577413249
	9	0.023855725390588
	10	0.019327134980565
	11	0.014510458353063
	12	9.52834300157261·10 -3
	13	4.69391346966794·10 -3
	14	7.41298074259234·10 -5
	15	-4.58204518385373·10 -3
	16	-9.51138868752708·10 -3

Проверка решения системы (фактически вычисляем A*X):

f2 :=	f21 ¬ b1×x1 + c1×x2
	for	i Î 2 .. n - 1
	f2i ¬ ai ×xi−1 + bi ×xi + ci ×xi+1
	f2n ¬ an×xn−1 + bn×xn
	f2
		1				1
	1	13			1	13
	2	12			2	12
	3	11			3	11
	4	10			4	10
	5	9			5	9
	6	8			6	8
	7	7			7	7
f2 =	8	6		f =	8	6
	9	5			9	5
	10	4			10	4
	11	3			11	3
	12	2			12	2
	13	1			13	1
	14	0			14	0
	15	-1			15	-1
	16	-2			16	-2

Находим вектор невязки:

v := f2 - f

		1
	1	0
	2	0
	3	0
	4	-1.77635683940025·10 -15
	5	0
	6	0
v =	7	0
	8	0
	9	0
	10	0
	11	0
	12	0
	13	0
	14	0
	15	0
	16	0

Находим нормы вектора невязки:

Norm1(v) = 5.175894435272 × 10− 13

Norm2(v) = 7.99838448334824 × 10− 14

NormInf(v) = 2.8421709430404 × 10− 14

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

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

Задача 2.4. Написать программы решения систем линейных уравнений AХ=В с диагональным преобладанием в матрице А методом Якоби и методом Зейделя, с модулем предварительной проверки достаточного условия сходимости итерационного процесса. Решить предложенную систему AХ=В при

помощи составленных программ с заданной точностью ε = 10−6 . Указать количество выполненных итераций. Составить векторы невязки и вычислить их норму.

Исходные данные:

ε := 10− 6

æ 5.94		0.8	0.6	−3.96	0.2	0.3	ö	æ 11.44		ö
ç				−2.97			÷	ç	−54.75	÷
ç 2.97		6.4	0		0.2	0.2	÷	ç		÷
ç	2.97	3.5	8.7	1.98	0.2	0	÷	ç	−4.64	÷
A := ç	4.95	1.6	1.2	−8.91	0.8	0.3	÷	B := ç	20.47	÷
ç							÷	ç		÷
ç	−0.99	2.5	1.1	−3.96	9	0.4	÷	ç	−95.86	÷
ç							÷	ç		÷
è 5.94		1.4	2.4	0	3.2	13 ø		è 26.92		ø