lab3_3

#LU-разложение
def get_LU(A):
    n = np.shape(A)[0]
    for i in range(n - 1):
        for j in range(i + 1, n):
            c = A[j][i] / A[i][i]
            A[j][i] = c
            for k in range(i+1, n):
                A[j][k] -= c * A[i][k]

    return A

# Решение СЛАУ
def Solve(LU, b):
    n = LU.shape[0]
    x = np.zeros(n)
    y = np.zeros(n)

    # Прямой ход: решение Lz=b
    def sum_y(i):
        res = 0
        for k in range(n):
            res += LU[i, k] * y[k]
        return res

    for i in range(n):
        y[i] = b[i] - sum_y(i)

    # Обратный ход: решение Ux=z
    def sum_x(i):
        res = 0
        for j in range(i, n):
            res += LU[i, j]*x[j]
        return res

    for i in reversed(range(n)):
        x[i] = (y[i] - sum_x(i))/LU[i, i]

    return x

def MNK1(X, f):
    n = np.shape(X)[0]

    A = np.array([[n, np.sum(X)],
                 [np.sum(X), np.sum(X*X)]])
    b = np.array([[np.sum(f)],
                  [np.sum(X*f)]])

    return Solve(get_LU(A), b)

def MNK2(X, f):
    n = np.shape(X)[0]

    A = np.array([[n, np.sum(X), np.sum(X*X)],
                  [np.sum(X), np.sum(X*X), np.sum(X*X*X)],
                  [np.sum(X*X), np.sum(X*X*X), np.sum(X*X*X*X)]])
    b = np.array([[np.sum(f)],
                  [np.sum(X*f)],
                  [np.sum(X*X*f)]])

    return Solve(get_LU(A), b)

def Error(f, polinom):
    return np.sqrt(np.sum((polinom-f)**2))

Xi = np.array([-0.7, -0.4, -0.1, 0.2, 0.5, 0.8])
fi = np.array([-0.7754, -0.41152, -0.10017, 0.20136, 0.5236, 0.9273])

a1_0, a1_1 = MNK1(Xi, fi)
a2_0, a2_1, a2_2 = MNK2(Xi, fi)

x = np.arange(-0.7, 0.8, 0.01)
print("Сумма квадратов ошибок:")
print("Для многочлена 1 степени:", Error(fi, np.array([float(a1_0+a1_1*xk) for xk in Xi])))
print("Для многочлена 2 степени:", Error(fi, np.array([float(a2_0+a2_1*xk+a2_2*xk*xk) for xk in Xi])))

plt.figure(figsize=(12, 7))
plt.scatter(Xi, fi, color="black", s=25)
plt.plot(x, [float(a1_0+a1_1*xk) for xk in x], label='Многочлен 1 степени', linewidth=1, color="blue")
plt.plot(x, [float(a2_0+a2_1*xk+a2_2*xk*xk) for xk in x], label='Многочлен 2 степени', linewidth=1, color="red")
plt.grid()
plt.legend()
plt.show()