LU_decomposition

# import matplotlib
import matplotlib.pyplot as plt
from numpy.linalg import inv
import numpy as np
from numpy import linspace
import sys
from math import fabs

'''
Функция, которая производит LU-разложение матрицы А и возвращает матрицы
L, U и P. Если permute = True, то разложение реализуется с частичным выбором
ведущего элемента (т.е. LUP-разложение).
'''
def lu(A, permute):
    n = len(A)
    U = np.copy(A)
    L = np.eye(n)
    P = np.eye(n)
    for k in range(n - 1):
        P_0 = np.eye(n)
        U_0 = np.eye(n)
        L_0 = np.eye(n)
        if (permute):
            max_elem = U[k][k].copy()
            i_max = k
            for i in range (k + 1, n):
                if abs(max_elem) < abs(U[i][k]):
                    max_elem = U[i][k].copy()
                    i_max = i
            P_0[i_max], P_0[k] = P_0[k].copy(), P_0[i_max].copy()
            U = P_0.dot(U)
        for i in range(k + 1, n):
            c = U[i][k] / U[k][k]
            U_0[i][k] = -c
            L_0[i][k] = c
        U = U_0.dot(U)
        P = P_0.dot(P)
        L = L.dot(P_0.dot(L_0))
    L = P.dot(L)
    return L, U, P

'''
Функция, которая возвращает решение СЛАУ Ах = b, где матрица
A представлена в виде LU-разложения с матрицей перестановок P,
т.е. PA = LU.
'''
def solve(L, U, P, b):
    n = len(L)
    x = [0] * n
    y = [0] * n
    for i in range(n):
        s = sum(L[i][j] * y[j] for j in range(i))
        y[i] = b[np.where(P[i] == 1)[0][0]] - s
    for i in reversed(range(n)):
        s = sum(U[i][j] * x[j] for j in range(i + 1, n))
        x[i] = (y[i] - s) / U[i][i]
    return x

'''
Функция, которая выводит на печать матрицу M.
'''
def print_matrix(M):
    for i in range(len(M)):
        print([" %.3f" % M[i][j] for j in range(len(M[0]))])
    print("\n")

'''
Функция, которая возвращает значение относительной погрешности
вычисления.
'''
def relative_error_E(x_circum, x_tilde):
    a = [abs(x_circum[i] - x_tilde[i]) for i in range(len(x))]
    b = [abs(x_circum[i]) for i in range(len(x))]
    return max(a) / max(b)


def get_A_b(a_1_1, b_1):
    A = [[a_1_1, 1., -3.],
         [6., 2., 5.],
         [1., 4., -3.]]
    b = [b_1, 12., -39.]
    return A, b

'''
Входные данные
'''
a_1_1 = 3
b_1 = -16
low_border = -12
high_border = 0
p = np.logspace(low_border, high_border, 80)

'''
Рассчёт
'''
A, b = get_A_b(a_1_1, b_1)
L, U, P = lu(A.copy(), True)
x = solve(L, U, P, b)
#print(x)

x_elim = []
x_part = []
E_elim = [0] * len(p)
E_part = [0] * len(p)

for p_i in p:
    A, b = get_A_b(a_1_1 + p_i, b_1 + p_i)
    L, U, P = lu(A.copy(), False)
    x_elim.append(solve(L, U, P, b))
#print(x_elim)

for p_i in p:
    A, b = get_A_b(a_1_1 + p_i, b_1 + p_i)
    L, U, P = lu(A.copy(), True)
    x_part.append(solve(L, U, P, b))
#print(x_part)

E_elim = [relative_error_E(x, x_elim[i]) for i in range(len(p))]
E_part = [relative_error_E(x, x_part[i]) for i in range(len(p))]

'''
Вывод результатов
'''
fig, axes = plt.subplots(nrows = 1, ncols = 1, figsize = (8, 4), dpi = 240)
axes.loglog(p, E_elim, "*", color="blue", label="$E{_{elim}}(p)$")
axes.loglog(p, E_part, '*', color="#551A8B", label="$E{_{part}}(p)$")
axes.grid()
axes.legend(loc='upper right')
fig.tight_layout()
plt.show()
fig.savefig("LU_decomposition.png")