№4 интерполяционный многочлен методом кубического сплайна

import numpy as np

# Функция для метода Гаусса
def gauss_method(A, b):
    n = len(b)
    # Прямой ход
    for i in range(n):
        # Поиск максимального элемента для избежания деления на малые числа
        max_row = i + np.argmax(np.abs(A[i:, i]))
        A[[i, max_row]] = A[[max_row, i]]  # Меняем строки
        b[[i, max_row]] = b[[max_row, i]]

        # Нормализуем строку
        div = A[i, i]
        A[i] /= div
        b[i] /= div

        # Вычитаем текущую строку из нижележащих
        for j in range(i + 1, n):
            factor = A[j, i]
            A[j] -= factor * A[i]
            b[j] -= factor * b[i]

    # Печать состояния матрицы после прямого хода
    print(f"\nМатрица после прямого хода:")
    extended_matrix = np.hstack([A, b.reshape(-1, 1)])
    for row in extended_matrix:
        for elem in row:
            print(f"{elem:>7.2f}", end=" ")
        print()

    # обратный ход
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - np.dot(A[i, i + 1:], x[i + 1:])) / A[i, i]

    print("\nОбратный ход:")
    for elem in x:
        print(f"{elem:>7.2f}", end=" ")
    print()

    return x

x = [0, 1, 2, 3]
f = [8, 9, 16, 43]
n = len(x) - 1

# расширенная матрица СЛАУ
A = np.zeros((4 * n, 4 * n))
b = np.zeros(4 * n)


# высчитываем значения, сначала в узлах a_i = f(x_i)
row = 0
for i in range(n):
    A[row, 4 * i] = 1
    b[row] = f[i]
    row += 1

# затем на концах отрезков
for i in range(n):
    h = x[i + 1] - x[i]
    A[row, 4 * i] = 1
    A[row, 4 * i + 1] = h
    A[row, 4 * i + 2] = h ** 2
    A[row, 4 * i + 3] = h ** 3
    b[row] = f[i + 1]
    row += 1

# уравнения для непрерывности первой производной
for i in range(n - 1):
    h = x[i + 1] - x[i]
    A[row, 4 * i + 1] = 1
    A[row, 4 * i + 2] = 2 * h
    A[row, 4 * i + 3] = 3 * h ** 2
    A[row, 4 * (i + 1) + 1] = -1
    row += 1

# уравнения для непрерывности второй производной
for i in range(n - 1):
    h = x[i + 1] - x[i]
    A[row, 4 * i + 2] = 2
    A[row, 4 * i + 3] = 6 * h
    A[row, 4 * (i + 1) + 2] = -2
    row += 1

# граничные условия: вторая производная на концах равна 0
A[row, 2] = 2
row += 1
A[row, 4 * (n - 1) + 2] = 2
A[row, 4 * (n - 1) + 3] = 6 * (x[n] - x[n - 1])

print("Таблица функции:")
for i in range(len(x)):
    print(f"{x[i]:>7.2f}", end=" ")
print()
for i in range(len(f)):
    print(f"{f[i]:>7.2f}", end=" ")
print()

print("\nРасширенная матрица СЛАУ:")
extended_matrix = np.hstack([A, b.reshape(-1, 1)])
for row in extended_matrix:
    for elem in row:
        print(f"{elem:>7.2f}", end=" ")
    print()

A_copy = A.copy()
b_copy = b.copy()
solution = gauss_method(A_copy, b_copy)

print("\nРешение СЛАУ:")
for elem in solution:
    print(f"{elem:>7.2f}", end=" ")
print()


# функция для вычисления значения сплайна в точке x
def spline_value(x_val, coeffs, x_points, index):
    a, b, c, d = coeffs[4 * index: 4 * index + 4]
    h = x_val - x_points[index]
    return a + b * h + c * h ** 2 + d * h ** 3


# точки для интерполяции (узлы, середины, конец отрезков)
interp_points = []
for i in range(n):
    interp_points.append(x[i])
    interp_points.append((x[i] + x[i + 1]) / 2)
interp_points.append(x[n])

# вычисление значений интерполяции
interp_values = []
for point in interp_points:
    index = min(n - 1, max(0, np.searchsorted(x, point) - 1))
    interp_values.append(spline_value(point, solution, x, index))

# Вывод результатов с форматированием
print("\nТаблица интерполянты:")
for point in interp_points:
    print(f"{point:>7.2f}", end=" ")
print()
for value in interp_values:
    print(f"{value:>7.2f}", end=" ")
print()