Практикум j практическое примщенше численных методов

^тк+1

^fc+l
w^k,r^k)
(.Aw^k,w^k)'
т_к+1 (w^k,r^k) 1 \
т_к (w^k~¹,r^k-^l)a_kJ
к = 1,2..., ct\ — 1.
Этот метод наиболее широко используется в вычислительной практике при решении задач с симметричной положительно определенной матрицей.

3. 
   Упражнения
   

Упражнение 5.1 Напишите программу, реализующую приблиэюеииое ре­шение системы линейных алгебраических уравнений итерационным мето­дом Зейделя. С ее помощью найдите решение системы уравнений с трех­диагональной матрицей
Ах = /,
в которой
вц 2, &г,г+1 1 QJ, {2>г,г—1 — 1 “Ь О!,
а правая часть
fo 1 fi б) i 2,3, ...,П 1, f_n — 1 "h qj,
определяет точное решение х^ = 1, г = 1,2, ...,п. Исследуйте зависимость числа итераций от п и параметра а при 0 < а < 1.
В модуле seidel функция seidelO обеспечивает итерационное решение си­стемы линейных уравнений методом Зейделя, результатом выступает при­ближенное решение и число итераций для достижения необходимой точности (нормы невязки).
:.У^::' л_::::":^;;|;::;^: ■
import numpy as np
import math as mt
def seidel(A, f, tol = 1.0e-9):
и и и
Solve the linear system Ax = b by Seidel method
и и и
n = len(f)
x = np.zeros((n), ’float’) r = np.copy(f)
# the maximum number of iterations is 10000 for k in ranged, 10001) : for i in range(n):
x[i] = x[i] + (f [i] - np.dot(A[i,0:n] , x[0:n])) / A[i,i] for i in range(n):
r[i] = f[i] - np.dot(A[i,0:n], x[0:n]) if np.dot(r,r) < tol**2: return x, k print ’Seidel method failed to converge:’
print ’after 10000 iteration residual = ’, mt.sqrt(np.dot(r,r)) return x, к
Решение модельной задачи дается следующей программой.
%хёг -5 1 РУ
import numpy as np from seidel import seidel n = 10 al = 0.5
A = np.zeros((n, n), ’float’) for i in range(0,n):
A [i, i] = 2 if i > 0:
A[i,i-1] = - 1 + al if i < n-1:
A[i,i+1] * - 1 -.al print ’A:\n’, A f = np.zeros((n), ’float’) f CO] = 1 - al f[n-l] = 1 + al print ’f:\n’, f x, iter = seidel(A, f) print ’iter =’, iter print ’x:\n’, x
1

iiilil	-1 5	0.	0.		^: '-M	0 ; !	0.	0 i	■Hi
piss	1 iiiiiii	-1 5	0			0 ;	0 .	, 0	HiH
Iiilil	—0.5	2	-15	0 ■		0	йО;;й]	0 ■. P	■ill!
(siiiis	issliii	-0 5	•Mh;	-1 5	: oSp	o. :S	ШИ	Illlll	illlll
illlill	1ЙК		-0.5	I; 2 pi	-1.5	-■oil	■ill	: ip	0
0	0	iliis	!l:!l	-0.5	2.	lilil	Illlill	SliSlII	0
pill	illlll	mm	.■lllll:	0 :	-0 5	’ 2 ■	Pil 51	siiiis	; 0
iiilil	iiiSlll	0	■lli	0	' '0 S »■	unis	■liii	■Ill	■III
piii	iltilllil	■Illlll	iiilil	■ill	iiiiiii	Illlll	Illlll	iiisti	lill
Illlll	m§m	Illlll	Iiilil	ills	111	Illlll	111!	iilli	lill
lill	illlill	IlSlll	Illlll	Illlll	11111	Illlll	Illlll	111	lill
llllp:	iiiiiss
1111!	Illlill	ISSlll	IIIIISS	iiilil	Illlll	Illlll	lill

11
Данные по числу итераций для рассматриваемой задачи с диагональным пре­обладанием суммированы в табл. 5.1

Таблица 5.1 Метод Зейделя

(\	n = 25	о Ю II e	n= 100
0	1233	4482	> 1000
0.25	289	410	599
0.5	107	159	254
0.75	58	93	160
1	25	50	100

При а ф () матрица системы уравнений несимметрична. Наблюдается вполне ожидаемое увеличение числа итераций при увеличении размерности задачи (параметр п) и сильная зависимость от параметра несимметричности а.
Упражнение 5.2 Напишите программу, реализующую приближенное ре­шение системы линейных алгебраических уравнений с симметричной поло- Э1сительио определенной матрицей методом сопряженных градиентов. С ее помощью найдите решение системы
Ах = f
с матрицей Гильберта
Q>ij ⁼ ~ j I 7 5 ^ 2,..., 71, j 1,2,п
и правой частью
п
fi ⁼ ^ ^ Q>ij ? * 1, 2, ..., 77,
3=1
для которой точное решение есть ж» = 1, i = 1,2, Исследуйте зави­симость числа итераций от п.
Расчетные формулы метода сопряженных градиентов возьмем (сравни с (5.24) при В = Е) в виде:
X^м = x^k+1 - a_k+ls^k, 5^fc+1 = r^k+1 + (3_k+ls^k, k = 0,1,...
при задании итерационных параметров a_k+i,/3_k+1 в виде
_ (r^k,r^k) _ (_rfe+i_>rfc+i) (r^k+1,As^k)
°^k+l ~ (s^k,As^k)’ ^{Л+1 =} (r^k,r^k) ⁼ ~ (s^k,As^k) '
В модуле eg функция eg О реализует описанный вариант метода сопряжен­ных градиентов для системы с самосопряженной положительно определенной матрицей А. На выходе мы имеем приближенное решение и число итераций для достижения необходимой точности по невязке.
import mimpy as np
def cg(A, f, tol = 1.0e-9):
к и и
Solve the linear system Ax = b by conjugate gradient method
и и и
n = len(f)
x = np.zeros((n), ’float¹) r = np.copy(f) for i in range(n):
r[i] = np.dot(A[i,0:n] , x[0:n]) - f[i] s = np.copy(r)
As = np.zeros((n), ’float’) for к in range(n):
for i in range(n):
As[i] = np.dot(A[i,0:n], s[0:n]) alpha = np.dot(r, r) / np.dot(s, As) x = x - alpha*s for i in range(n):
r[i] = np.dot(A[i,0:n], x[0:n]) - f[i] if np.dot(r, r) < tol**2: break else:
beta = - np.dot(r, As) / np.dot(s, As) s = r + beta*s return x, к
Решение тестовой задачи дается следующей программой.
import numpy as np from eg import eg n = 4
A = np.zeros((n, n), ’float’) for i in range(0,n):
for j in range(0,n):
A[i,j] = l./(i+j+l) print ’A:\n’, A f = np.ones((n), ’float’) for i in range(0,n): f[i] = 0.
for j in range(0,n):
f [i] = f[i] + A [i, j] print ’f:\n’, f x, iter = cg(A, f) print ’iter =’, iter print ’x:\n’, x
| [ 1 ^: 0.33333333 0 25 ]
[ 0.33333333 0,25 : Q 16666667]





Зависимость числа итераций от размерности матрицы прослеживается по следующим данным: при п = 25 необходимо 13 итераций, при п — 50 — 17, а при п = 100 — 20.

4. 
   Задачи
   

Задача 5.1 Напишите программу, реализующую приближенное решение системы линейных алгебраических уравнений итерационным методом Яко­би. С ее помощью найдите решение системы уравнений, рассмотренной в упразднении 5.1. Сравните скорости сходимости итерационных методов Якоби и Зейделя при различных параметрах п и а.
Задача 5.2 Напишите программу, реализующую приближенное решение системы линейных алгебраических уравнений методов релаксации (5.9). Ис­следуйте зависимость скорости сходимости этого итерационного метода от итерационного парсиметра Tk+\ = г при численном решении системы уравнений из упразднения 5.1 при различных п и а.
Задача 5.3 Напишите программу для решения системы линейных алгебра­ических уравнений явным (В = Е) итерационным методом минимальных поправок (итерационные параметры определяются согласно (5.23)). Сравни­те скорость сходимости этого метода со скоростью сходимости метода сопряженных градиентов на примере задачи из упражнения 5.2.
Задача 5.4 Пусть
А = А₁ + А₂ = А*> о, А\ = А, = ^l-D + L,
В попеременно-треугольном методе переобуславливатель В зададим в виде В = (Е + изА\)(Е +
Напишите программу для решения системы линейных алгебраических урав­нений с симметричной положительно определенной матрицей поперемен­но-треугольным методом с выбором итерационных параметров по (5.23). Исследуйте зависимость скорости сходимости этого итерационного мето­да от параметра из в задаче из упрсю1спепия 5.2.
Спектральные задачи линейной алгебры
Важной и трудной задачей линейной алгебры является нахождение собствен­ных значений и собственных векторов матриц. Рассматривается проблема устойчивости собственных значений по отношению к малым возмущениям элементов матрицы. Для приближенного нахождения отдельных собствен­ных значений широко используется степенной метод в различных модифика­циях. Для решения полной проблемы для симметричных матриц применяется итерационный метод Якоби и QR-алгоритм.

Основные обозначения
X \Х{ } {^1) 3-2 > • • • > Хп }	— n-мерный вектор
А={ац)	— матрица с вещественными элементами ац
Е	— единичная матрица
D = diag{di,d2,.. .,dn}	— диагональная матрица
Aj, 2 — 1,2, . • • ,72	— собственные значения
<р\ г = 1,2,... ,n n	— собственные вектора
(x, y) = ^2XiVi i= 1	— скалярное произведение

1. 
   Собственные значения и собственные вектора матриц
   

Рассматриваются проблемы нахождения собственных значений и собствен' ных векторов квадратной вещественной матрицы А. Собственным числом называется число Л такое, что для некоторого ненулевого вектора (собствен' ного вектора) р имеет место равенство
Собственные вектора определены с точностью до числового множителя. Мно­жество всех собственных значений матрицы А называется спектром матри­цы А.
С учетом того, что ищется нетривиальное решение уравнения (6.1), то
det(i4 - АЕ) = 0. (6.2)
Тем самым собственные значения А матрицы А являются корнями характе­ристического многочлена n-й степени (6.2). Задача отыскания собственных значений и собственных векторов матрицы сводится к построению характери­стического многочлена, отысканию его корней и последующему нахождению нетривиальных решений уравнения (6.1) для найденных собственных значе­ний.
В вычислительной практике рассматривается как полная проблема собствен­ных значений, когда необходимо найти все собственные значения матрицы А, так и частичная проблема собственных значений, когда ищутся только некоторые собственные значения, например, максимальные (минимальные) по модулю.

2. 
   Численные методы решения задач на собственные значения
   

Начнем с приведением некоторых полезных фактов о свойствах собственных значений и собственных векторов квадратной матрицы. Далее рассматрива­ются методы решения частичной и полной проблемы собственных значений.
Свойства собственных значений и собственных векторов
Квадратная вещественная матрица порядка п имеет п собственных значе­ний, при этом каждое собственное значение считается столько раз, какова его кратность как корня характеристического многочлена. Для симметрич­ной матрицы А собственные значения вещественны, а собственные векто­ра, соответствующие различным собственным значениям, ортогональны, т.е. (<р^г, ip*) =0, если г ф j.
Отметим также некоторые свойства собственных значений и собственных векторов для сопряженной матрицы А*:
А
(6.3)

*ф = рф.
Для спектральных задач (6.1), (6.3) имеем
Аг Pi) Ъ 1, 2, . . . , 72,
(<р^,,ф¹)= 0, i ф j.
Две квадратные матрицы А и В одинаковых размеров называются по ными, если существует такая невырожденная матрица Р, что А = Р~¹ Подобные матрицы имеют одни и те же собственные значения, так ка] (6.1) непосредственно следует
Вф = \ф, ф = Р^р.
Поэтому вычислительные алгоритмы решения спектральных задач баз] ются на подобном преобразовании матрицы к матрице В, для которой с тральная задача решается проще. В качестве В естественно выбирать ди нальную матрицу, причем в данном случае это будет
A = diag{A_bA₂,...,A_n}.
Упорядочим собственные значения симметричной матрицы А по возра нию, т.е. Ai < Л₂ < • • • < Л_п. Свойства собственных значений и собствен
(Ах х)
функций связаны с отношением Релея ~Отметим, например, что
(ж, х)
любого ненулевого вектора х справедливы неравенства
(х,х)
В