ХО ; 1 093306841
Приведенные расчетные данные демонстрируют быстрое приближение минимального собственного значения к нулю при увеличении размерности матрицы п.
Упражнение 6.2 Напишите программу для пахооюдения собственных значений и собственных векторов симметричной вещественной матрицы методом вращений (методом Якоби). С ее помощью найдите собственные значение матрицы Лемера (Lehmer) А, для которой
г
min (i,j)
тах(г, j) ’
= 1,2, j = 1,2, ...,п,
при п = 8.
Приведем необходимые расчетные формулы, которые связаны с применением матрицы вращения. Пусть
В = T*(kl)AT(kl).
Для измененных элементов матрицы В имеем
bkk = kk + s2a,t - 2сзаы,
Ьц = s2a,kk + с2 o,a + 2 csakiy
192 Спектральные задачи линейной алгебру
bid = (с2 - s2)akt + cs(akk- аи),
Ьы = caki - sau, гфк, i ф l,
Ьц = saki + can, гфк, i ф l.
В силу этого, для обнуления элемента Ьы необходимо (с2 - s2)akl + cs(akk - аи) = 0.
Пусть с = cos в, s = sin в {в — угол поворота) тогда достаточно положить
I _ 1 _ акк — аи
tan(20) 2 аы
Полагая t = s/c = tan 0, для определения t получим квадратное уравнение
t2 + 2 - 1 = 0,
решение которого есть
ь = -ф± у/ф2ТТ.
Меньший по модулю корень \t\ < 1 соответствует повороту на угол меньший 7г/4 и порождает более устойчивый алгоритм. Поэтому
t = sign(<£)(-|0| + фф2 + 1).
Для более точных вычислений при больших \ф\ положим
= sign(0)
\Ф\ + ФФ2 + г
причем для очень больших \ф\ имеем t « 1/(20).
По значению t находим
1
S = tc.
С = —,
фТ+t2
Расчетные формулы для элементов матрицы В принимают вид
Ькк — &кк tQ'kh
Ьц — an + tau,
Ьы — 0?
Ьы = аы — s(an -Ь тс1ы)) ^ ф к, % ф /,
Ьц = ац + з(аы-тан), гфк, гф1,
где т = s/(l + с).
В модуле jacobi функция elemMaxO находит максимальный по модулю эле- мент матрицы, для обнуления отдельного элемента матрицы используется функция rotate().
Модуль Jacobi
import пшпру as np import math as mt jef 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) показывает хорошее согласие результатов.
Задачи
Задача 6.1 Напишите программу для нахоэюдепия ближайшего к заданному числу собственного значения и соответствующего собственного вектора симметричной веществетюй матрицы при использовании обратных итераций со сдвигом. С ее помощью найдите первые три минимальных по модулю собственных значение матрицы Паскаля А, для которой
aii = [i - l)!(j - 1)! ® = >п» j = l,2, ...,n.
при n — 8.
Задача 6.2 Напишите программу для нахоэюдепия минимального по модулю собственного значения и соответствующего собственного вектора симметричной трехдиагональной матрицы. С ее помощью найдите первое минимальное по модулю собственное значение матрицы п х п, для которой
Оц 2, <2^2 — 1 Qf2t2+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,...,п
|
— функции п переменных
|
|
|
|
Do'stlaringiz bilan baham: |