北航数值分析作业二

import cmath
from math import *
sgn = lambda x: 1 if x > 0 else -1


# 坑爹的赋值问题!注意下标从1开始
def init():
   A = [[0 for i in range(10)] for i in range(10)]
   for i in range(10):
      for j in range(10):
         ii, jj = i + 1, j + 1
         A[i][j] = sin(0.5 * ii + 0.2 * jj) if ii != jj else 1.52 * cos(ii + 1.2 * jj)
   return A


# 矩阵转置
def transposition(A):
   ans = [[0 for i in range(len(A))] for j in range(len(A[0]))]
   for i in range(len(A)):
      for j in range(len(A[i])):
         ans[j][i] = A[i][j]
   return ans


# 矩阵与矩阵相乘
def matrix_mul_matrix(A, B):
   assert len(A[0]) == len(B)
   ans = [[0] * len(B[0]) for i in range(len(A))]
   for i in range(len(ans)):
      for j in range(len(ans[i])):
         ans[i][j] = sum([A[i][k] * B[k][j] for k in range(len(B))])
   return ans


# 矩阵与向量相乘
def matrix_mul_vector(A, b):
   assert len(A[0]) == len(b)
   return [sum([A[i][j] * b[j] for j in range(len(b))]) for i in range(len(A))]


# 向量与向量相乘
def vector_mul_vector(A, B):
   assert len(A) == len(B)
   return sum([A[i] * B[i] for i in range(len(A))])


# 向量除以常数k
def div(v, k):
   return [v[i] / k for i in range(len(v))]


# 矩阵的拟上三角化
def quasi_upper_triangular(A):
   for i in range(0, len(A) - 2):  # i表示列数
      if not any([A[j][i] for j in range(i + 2, len(A))]): continue
      d = sqrt(sum([A[j][i] ** 2 for j in range(i + 1, len(A))]))
      c = -sgn(A[i + 1][i]) * d
      h = c * (c - A[i + 1][i])
      u = [0] * (i + 1) + [A[i + 1][i] - c] + [A[j][i] for j in range(i + 2, len(A))]
      p = div(matrix_mul_vector(transposition(A), u), h)
      q = div(matrix_mul_vector(A, u), h)
      t = vector_mul_vector(p, u) / h
      w = [q[i] - t * u[i] for i in range(len(q))]

      for i in range(len(A)):
         for j in range(len(A[i])):
            A[i][j] = A[i][j] - w[i] * u[j] - u[i] * p[j]


# 打印矩阵
def printMatrix(A,s):
   print(s)
   for i in range(len(A)):
      for j in range(len(A[i])):
         print('%.3f' % A[i][j], ' ', end=',')
      print()
   print("=========")


def qr(B, C):
   for i in range(len(B) - 1):
      if not any([B[j][i] for j in range(i + 1, len(B))]): continue
      d = sqrt(sum([B[j][i] ** 2 for j in range(i, len(B))]))
      c = -sgn(B[i][i]) * d
      h = c * (c - B[i][i])
      u = [0] * (i) + [B[i][i] - c] + [B[j][i] for j in range(i + 1, len(B))]
      v = div(matrix_mul_vector(transposition(B), u), h)
      for i in range(len(B)):
         for j in range(len(B[i])):
            B[i][j] -= u[i] * v[j]
      p = div(matrix_mul_vector(transposition(C), u), h)
      q = div(matrix_mul_vector(C, u), h)
      t = vector_mul_vector(p, u) / h
      w = [q[i] - t * u[i] for i in range(len(q))]
      for i in range(len(C)):
         for j in range(len(C[i])):
            C[i][j] = C[i][j] - w[i] * u[j] - u[i] * p[j]


# 带双步位移的QR分解求特征根
def qr_with_double_shift(A):
   root = [0] * len(A)
   i = len(A) - 1
   k = 0  # 步数
   while i >= 0:
      if i == 0:
         root[0] = A[0][0]
         break
      elif abs(A[i][i - 1]) < epsilon:
         root[i] = A[i][i]
         i -= 1
         continue
      d1, d2, d3, d4 = A[i - 1][i - 1], A[i - 1][i], A[i][i - 1], A[i][i]
      delta = cmath.sqrt((d1 + d4) ** 2 - 4 * (d1 * d4 - d2 * d3))
      if i == 1 or abs(A[i - 1][i - 2]) < epsilon:
         (root[i], root[i - 1]) = (((d1 + d4) + delta) / 2, ((d1 + d4) - delta) / 2)
         i -= 2
         continue
      if k == L:
         print("didn't get all of the root after {} steps".format(L))
         break
      s, t = d1 + d4, d1 * d4 - d2 * d3
      M = matrix_mul_matrix(A, A)
      for i in range(len(M)):
         for j in range(len(M[0])):
            M[i][j] -= s * A[i][j]
         M[i][i] += t
      qr(M, A)
      k += 1
   return root


"""
高斯消去法适用范围:各阶主子式大于0
列主元高斯消去法适用范围:行列式值大于0
全主元高斯消去法适用范围:求解任意方程,可以求出一个解空间来
在本问题中,根据特征值求特征向量有两种方法:
   1.零点平移反幂法迭代求最接近lamda的特征值,也能求出特征向量,但它只能求出一个特征向量
   2.全主元高斯消去法,这种方法最完善
"""


# 全主元高斯消去法求方阵A关于特征根root的特征向量
def solve(A, root):
   n = len(A)
   for i in range(n):
      A[i][i] -= root
   ind = [i for i in range(n)]
   rank = n
   for i in range(n):
      x, y = i, i
      for row in range(i, n):
         for col in range(i, n):
            if abs(A[x][y]) < abs(A[row][col]):
               x, y = row, col
      if abs(A[x][y]) < epsilon:  # 最大值也为0,开始回代
         rank = i
         break
      # 将最大行与当前行进行交换
      A[x], A[i] = A[i], A[x]
      # 换列,第y列和第i列
      for row in range(n):
         A[row][y], A[row][i] = A[row][i], A[row][y]
      ind[i], ind[y] = ind[y], ind[i]
      # 单位化一行
      t=A[i][i]
      for j in range(i, n):
         A[i][j] /= t
      for row in range(i + 1, n):  # 第j行-A[j][i]倍的第i行
         t = A[row][i]
         if t==0:continue
         for col in range(i,n):
            A[row][col]-=t*A[i][col]
   # 回代过程
   for row in range(rank - 1, -1, -1):
      for j in range(row + 1, n):  # 第i行-A[i][j]倍的第j行
         # 这里一定要注意必须用t把A[i][j]存起来,否则A[i][j]就变成0了
         t = A[row][j]
         for k in range(j, n):
            A[row][k] -= t * A[j][k]
   # 构造特征向量空间,一定要注意把各个列给换回去!
   ans = [[0] * len(A) for i in range(n - rank)]
   for i in range(rank, n):
      for j in range(n):
         ans[i - rank][ind[j]] = -A[j][i]
      ans[i - rank][ind[i]] = 1
   return ans


epsilon = 1e-12
L = 0xffffff

A = init()
printMatrix(A,'最初的矩阵')
quasi_upper_triangular(A)
printMatrix(A,'拟对角化之后的矩阵')
roots = qr_with_double_shift(A)
print('各个特征根',roots)
printMatrix(A,'qr双步位移法结束之后的矩阵')
for i in roots:
   if type(i) == float:
      A = init()
      x = solve(A, i)
      print('特征根', i, '对应的特征向量为
', x)