daily pastebin goal
24%
SHARE
TWEET

LU_decomposition

a guest Jun 10th, 2018 124 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # import matplotlib
  2. import matplotlib.pyplot as plt
  3. from numpy.linalg import inv
  4. import numpy as np
  5. from numpy import linspace
  6. import sys
  7. from math import fabs
  8.  
  9. '''
  10. Функция, которая производит LU-разложение матрицы А и возвращает матрицы
  11. L, U и P. Если permute = True, то разложение реализуется с частичным выбором
  12. ведущего элемента (т.е. LUP-разложение).
  13. '''
  14. def lu(A, permute):
  15.     n = len(A)
  16.     U = np.copy(A)
  17.     L = np.eye(n)
  18.     P = np.eye(n)  
  19.     for k in range(n - 1):
  20.         P_0 = np.eye(n)
  21.         U_0 = np.eye(n)
  22.         L_0 = np.eye(n)
  23.         if (permute):
  24.             max_elem = U[k][k].copy()
  25.             i_max = k
  26.             for i in range (k + 1, n):
  27.                 if abs(max_elem) < abs(U[i][k]):
  28.                     max_elem = U[i][k].copy()
  29.                     i_max = i
  30.             P_0[i_max], P_0[k] = P_0[k].copy(), P_0[i_max].copy()
  31.             U = P_0.dot(U)
  32.         for i in range(k + 1, n):
  33.             c = U[i][k] / U[k][k]
  34.             U_0[i][k] = -c
  35.             L_0[i][k] = c
  36.         U = U_0.dot(U)
  37.         P = P_0.dot(P)
  38.         L = L.dot(P_0.dot(L_0))
  39.     L = P.dot(L)    
  40.     return L, U, P
  41.  
  42. '''
  43. Функция, которая возвращает решение СЛАУ Ах = b, где матрица
  44. A представлена в виде LU-разложения с матрицей перестановок P,
  45. т.е. PA = LU.
  46. '''
  47. def solve(L, U, P, b):
  48.     n = len(L)
  49.     x = [0] * n
  50.     y = [0] * n
  51.     for i in range(n):
  52.         s = sum(L[i][j] * y[j] for j in range(i))
  53.         y[i] = b[np.where(P[i] == 1)[0][0]] - s
  54.     for i in reversed(range(n)):
  55.         s = sum(U[i][j] * x[j] for j in range(i + 1, n))
  56.         x[i] = (y[i] - s) / U[i][i]
  57.     return x
  58.  
  59. '''
  60. Функция, которая выводит на печать матрицу M.
  61. '''
  62. def print_matrix(M):
  63.     for i in range(len(M)):
  64.         print([" %.3f" % M[i][j] for j in range(len(M[0]))])
  65.     print("\n")
  66.  
  67. '''
  68. Функция, которая возвращает значение относительной погрешности
  69. вычисления.
  70. '''
  71. def relative_error_E(x_circum, x_tilde):
  72.     a = [abs(x_circum[i] - x_tilde[i]) for i in range(len(x))]
  73.     b = [abs(x_circum[i]) for i in range(len(x))]
  74.     return max(a) / max(b)
  75.  
  76.  
  77. def get_A_b(a_1_1, b_1):
  78.     A = [[a_1_1, 1., -3.],
  79.          [6., 2., 5.],
  80.          [1., 4., -3.]]
  81.     b = [b_1, 12., -39.]
  82.     return A, b
  83.  
  84. '''
  85. Входные данные
  86. '''  
  87. a_1_1 = 3
  88. b_1 = -16
  89. low_border = -12
  90. high_border = 0
  91. p = np.logspace(low_border, high_border, 80)
  92.  
  93. '''
  94. Рассчёт
  95. '''
  96. A, b = get_A_b(a_1_1, b_1)
  97. L, U, P = lu(A.copy(), True)
  98. x = solve(L, U, P, b)
  99. #print(x)
  100.  
  101. x_elim = []
  102. x_part = []
  103. E_elim = [0] * len(p)
  104. E_part = [0] * len(p)
  105.  
  106. for p_i in p:
  107.     A, b = get_A_b(a_1_1 + p_i, b_1 + p_i)
  108.     L, U, P = lu(A.copy(), False)
  109.     x_elim.append(solve(L, U, P, b))
  110. #print(x_elim)
  111.  
  112. for p_i in p:
  113.     A, b = get_A_b(a_1_1 + p_i, b_1 + p_i)
  114.     L, U, P = lu(A.copy(), True)
  115.     x_part.append(solve(L, U, P, b))
  116. #print(x_part)
  117.  
  118. E_elim = [relative_error_E(x, x_elim[i]) for i in range(len(p))]
  119. E_part = [relative_error_E(x, x_part[i]) for i in range(len(p))]
  120.  
  121. '''
  122. Вывод результатов
  123. '''
  124. fig, axes = plt.subplots(nrows = 1, ncols = 1, figsize = (8, 4), dpi = 240)
  125. axes.loglog(p, E_elim, "*", color="blue", label="$E{_{elim}}(p)$")
  126. axes.loglog(p, E_part, '*', color="#551A8B", label="$E{_{part}}(p)$")
  127. axes.grid()
  128. axes.legend(loc='upper right')
  129. fig.tight_layout()
  130. plt.show()
  131. fig.savefig("LU_decomposition.png")
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top