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

ХО ; 1 093306841
Приведенные расчетные данные демонстрируют быстрое приближение мини­мального собственного значения к нулю при увеличении размерности матри­цы п.
Упражнение 6.2 Напишите программу для пахооюдения собственных зна­чений и собственных векторов симметричной вещественной матрицы ме­тодом вращений (методом Якоби). С ее помощью найдите собственные зна­чение матрицы Лемера (Lehmer) А, для которой
г
min (i,j)
тах(г, j) ’
= 1,2, j = 1,2, ...,п,
при п = 8.
Приведем необходимые расчетные формулы, которые связаны с применением матрицы вращения. Пусть
В = T*(kl)AT(kl).
Для измененных элементов матрицы В имеем
bkk = kk + s²a,t - 2сза_ы,
Ьц = s²a,kk + с² o,a + 2 csakiy
192 Спектральные задачи линейной алгебру
bid = (с² - s²)a_kt + cs(a_kk- а_и),
Ьы = ca_ki - sa_u, гфк, i ф l,
Ьц = sa_ki + can, гфк, i ф l.
В силу этого, для обнуления элемента Ьы необходимо (с² - s²)a_kl + cs(a_kk - а_и) = 0.
Пусть с = cos в, s = sin в {в — угол поворота) тогда достаточно положить
I _ 1 _ ^акк — аи
tan(20) 2 а_ы
Полагая t = s/c = tan 0, для определения t получим квадратное уравнение
t² + 2 - 1 = 0,
решение которого есть
ь = -ф± у/ф²ТТ.
Меньший по модулю корень \t\ < 1 соответствует повороту на угол меньший 7г/4 и порождает более устойчивый алгоритм. Поэтому
t = sign(<£)(-|0| + фф² + 1).
Для более точных вычислений при больших \ф\ положим
₌ sign(0)
\Ф\ + ФФ² + г
причем для очень больших \ф\ имеем t « 1/(20).
По значению t находим
1
S = tc.


С = —,
фТ+t²
Расчетные формулы для элементов матрицы В принимают вид
Ькк — &кк tQ'kh
Ьц — an + tau,
Ьы — 0?
Ьы ⁼ аы — s(an -Ь тс1ы)) ^ ф к, % ф /,
Ьц = ац + з(аы-та_н), гфк, гф1,
где т = s/(l + с).
В модуле jacobi функция elemMaxO находит максимальный по модулю эле- мент матрицы, для обнуления отдельного элемента матрицы используется функция rotate().
Модуль Jacobi
import пшпру as np import math as mt j_ef elemMax(A):
к и и
Find the largest (absolute value) off-diagonal element A[k,l] in the upper half of A.
и и и
n = len(A) aMax = 0.
for i in range(n-l):
for j in range(i+1, n):
if abs(A[i,j]) >= aMax: aMax = abs(A[i,j]) к = i
1 = j
return aMax, к, 1 def rotate(A, P, к, 1):
и и и
Rotate of A to make A[k,l] =0.
и и и
n = len(A) d = A [1,1] - A [k, k] if abs(A[k,l]) < abs(d)*l.Oe-36: t = A[k,l] / d else:
phi = d / (2*A[k,l])
t = 1 / (abs(phi) + mt.sqrt(l + phi**2)) if phi < 0.: t = - t
с = 1 / mt.sqrt(l + t**2) s = t*c
tau = s / (1 + c)
# Modify the matrix elements tt = A[k,l]
A[k,l] = 0.0 A[k,k] = A[k,k] - t*tt A [1,1] = A [1,1] + t*tt for i in range(k): tt = A[i,k]
A[i,k] = tt - s*(A[i,l] + tau*tt)
A[i, 1] = A[i, 1] + s*(tt - tau*A[i,l]) for i in range(k+1,1): tt = A[k,i]
А [к, i] = tt - s*(A[i,l] + tau*A[k,i])
A[i, 1] = A[i,l] + s*(tt - tau*A[i,l]) for i in range(l+l,n): tt = A[k,i]
A[k,i] = tt - s*(A[l,i] + tau*tt)
A[l,i] = A[l,i] + s*(tt - tau*A[l,i])

• 
  Update transformation matrix for i in range(n):
  

tt = P[i,k]
P[i,k] = tt - s*(P[i,l] + tau*P[i,k])
P[i,l] = P[i,l] + s*(tt - tau*P[i,l]) def jacobi(A, tol = 1.0e-12):
к и и
Solution of eigenvalue problem by Jacobi’s method.
Returns eigenvalues in vector lam
and the eigenvectors as columns of matrix .
и и и
n = len(A)

• 
  Number of rotations limit rotMax = 5*(n**2)
  

P = np.identity(n)

• 
  Jacobi rotation loop for i in range(rotMax):
  

aMax, к, 1 = elemMax(A) if aMax < tol:
return np.diagonal(A), P rotate(A, P, к, 1)
print ’Jacobi method did not converge’
Решение полной проблемы собственных значений для матрицы Гильбе при использовании метода Якоби дается следующей программой.
import numpy as np from jacobi import jacobi n = 8
A = np.zeros((n, n), ’float’) for i in range(0, n):
for j in range(0, n): if i < j:
A[i, j] = (i+1.)/(j+1) else:
A[i, j] = (j+1.)/(i+1) lam, x = jacobi(A) print ’All eigenvalue:\n’, lam
[ о 74602786 1 46486865 0 2760521 0 439372
О 12603956 0 18325435 0.087072^6 ; 4
Сравнение с результатами вычисления минимального собственного значения методом Якоби и обратными итерациями (см. упражнение 6.1) показывает хорошее согласие результатов.

4. 
   Задачи
   

Задача 6.1 Напишите программу для нахоэюдепия ближайшего к задан­ному числу собственного значения и соответствующего собственного век­тора симметричной веществетюй матрицы при использовании обратных итераций со сдвигом. С ее помощью найдите первые три минимальных по модулю собственных значение матрицы Паскаля А, для которой
^aii ⁼ [i - l)!(j - 1)! ® ⁼ >^п» j = l,2, ...,n.
при n — 8.
Задача 6.2 Напишите программу для нахоэюдепия минимального по моду­лю собственного значения и соответствующего собственного вектора сим­метричной трехдиагональной матрицы. С ее помощью найдите первое ми­нимальное по модулю собственное значение матрицы п х п, для которой
Оц 2, <2^2 — 1 Qf2_t2+1 1
при различных п. Сравните численное решение с точным.
Задача 6.3 Напишите программу для преобразования симметричной веще­ственной матрицы А к симметричной трехдиагональной матрице РАР с помощью ортогонального преобразования Хаусхольдера (Р^-1 = Р). С ее по­мощью проведите преобразование матрицы А, элементы которой есть
а_ц = min(i,j), г = l,2,...,n, j = 1,2
при различных п.
Задача 6.4 Напишите программу для вычисления собственных значений трехдиагональной вещественной матрицы А с использованием QR алго­ритма. Найдите собственные числа трехдиагопальпой матрицы А, для ко­торой
о>а — 2, &г,г+1 — 1 oi, 6X2,1—1 — 1 О!,
при различных значениях п и параметра а (0 < а < 1). Сравните найденные собственные значения с точными.
Нелинейные уравнения и системы
Многие прикладные задачи приводят к необходимости нахождения прибли­женного решения нелинейных уравнений и систем нелинейных уравнений. С этой целью используются итерационные методы. Приведены алгоритмы решения нелинейных уравнений с одним неизвестным и систем нелинейных уравнений. Применяются итерационные метод последовательных приближе­ний (простой итерации) и метод Ньютона в различных модификациях.

	Основные обозначения
	/(*)	— функция одной переменной
	, г = 1,2,...,п	— функции п переменных