Метод Якоби, выводящий неправильные собственные значения - PullRequest
Новая
       16

Метод Якоби, выводящий неправильные собственные значения

0 голосов
/ 17 октября 2019

Я работаю над созданием калькулятора собственных значений с использованием метода Якоби, и он работает без ошибок. Однако он не находит правильных собственных значений и не находит правильных собственных векторов. По некоторым причинам, я всегда получаю собственные значения 0. Я думаю, что это может не сохранять матрицу, которую я вводил для MatrixA. (Ссылка на метод Якоби, если вы не знакомы: http://fourier.eng.hmc.edu/e176/lectures/ch1/node1.html)

import numpy as np
import bettertimeit as time
import matplotlib as plt

def Jacobi(A):
    n     = A.shape[0]            # matrix size #columns = #lines
    maxit = 100                   # maximum number of iterations
    eps   = 1.0e-15               # accuracy goal
    pi    = np.pi
    info  = 0                     # return flag
    ev    = np.zeros(n,float)     # initialize eigenvalues
    U     = np.zeros((n,n),float) # initialize eigenvector
    for i in range(0,n): U[i,i] = 1.0

    for t in range(0,maxit):
         s = 0    # compute sum of off-diagonal elements in A(i,j)
         for i in range(0,n): s = s + np.sum(np.abs(A[i,(i+1):n]))
         if (s < eps): # diagonal form reached
              info = t
              for i in range(0,n):ev[i] = A[i,i]
              break
         else:
              limit = s/(n*(n-1)/2.0)       # average value of off-diagonal elements
              for i in range(0,n-1):       # loop over lines of matrix
                   for j in range(i+1,n):  # loop over columns of matrix
                       if (np.abs(A[i,j]) > limit):      # determine (ij) such that |A(i,j)| larger than average
                                                         # value of off-diagonal elements
                           denom = A[i,i] - A[j,j]       # denominator of Eq. (3.61)
                           if (np.abs(denom) < eps): phi = pi/4         # Eq. (3.62)
                           else: phi = 0.5*np.arctan(2.0*A[i,j]/denom)  # Eq. (3.61)
                           si = np.sin(phi)
                           co = np.cos(phi)
                           for k in range(i+1,j):
                               store  = A[i,k]
                               A[i,k] = A[i,k]*co + A[k,j]*si  # Eq. (3.56)
                               A[k,j] = A[k,j]*co - store *si  # Eq. (3.57)
                           for k in range(j+1,n):
                               store  = A[i,k]
                               A[i,k] = A[i,k]*co + A[j,k]*si  # Eq. (3.56)
                               A[j,k] = A[j,k]*co - store *si  # Eq. (3.57)
                           for k in range(0,i):
                               store  = A[k,i]
                               A[k,i] = A[k,i]*co + A[k,j]*si
                               A[k,j] = A[k,j]*co - store *si
                           store = A[i,i]
                           A[i,i] = A[i,i]*co*co + 2.0*A[i,j]*co*si +A[j,j]*si*si  # Eq. (3.58)
                           A[j,j] = A[j,j]*co*co - 2.0*A[i,j]*co*si +store *si*si  # Eq. (3.59)
                           A[i,j] = 0.0                                            # Eq. (3.60)
                           for k in range(0,n):
                                store  = U[k,j]
                                U[k,j] = U[k,j]*co - U[k,i]*si  # Eq. (3.66)
                                U[k,i] = U[k,i]*co + store *si  # Eq. (3.67)
         info = -t # in case no convergence is reached set info to a negative value "-t"
    return ev,U,t

n = int(input("Enter the matrix size: "))

A = np.zeros((n, n))

for i in range(n):
    A[i] = input().split(" ")

MatrixA = np.array(A)
print("A= ")
print(A)
for i in range(A.shape[0]):
    row = ["{}*x{}".format(A[i, j], j + 1) for j in range(A.shape[1])]

# Jacobi-method
ev,U,t = Jacobi(A)
print ("JACOBI METHOD: Number of rotations = ",t)
print ("Eigenvalues  = ",ev)
print ("Eigenvectors = ")
print (U)
...