Iнтерполяцiя кубiчним сплайном

import numpy as np
import matplotlib.pyplot as plt

def tridiagonal_solve(a, b, c, d):
    n = len(d)
    c_ = np.zeros(n)
    d_ = np.zeros(n)
    x = np.zeros(n)

    c_[0] = c[0] / b[0]
    d_[0] = d[0] / b[0]

    for i in range(1, n):
        c_[i] = c[i] / (b[i] - a[i] * c_[i - 1])
        d_[i] = (d[i] - a[i] * d_[i - 1]) / (b[i] - a[i] * c_[i - 1])

    x[n - 1] = d_[n - 1]

    for i in range(n - 2, -1, -1):
        x[i] = d_[i] - c_[i] * x[i + 1]

    return x

def cubic_spline(x, y):
    n = len(x)
    h = np.diff(x)
    alpha = np.zeros(n)
    for i in range(1, n - 1):
        alpha[i] = (3 / h[i]) * (y[i + 1] - y[i]) - (3 / h[i - 1]) * (y[i] - y[i - 1])
    L = np.zeros(n)
    mu = np.zeros(n)
    z = np.zeros(n)
    c = np.zeros(n)

    L[0] = 1.0
    mu[0] = 0.0
    z[0] = 0.0

    for i in range(1, n - 1):
        L[i] = 2.0 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1]
        mu[i] = h[i] / L[i]
        z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / L[i]

    L[n - 1] = 1.0
    z[n - 1] = 0.0
    c[n - 1] = 0.0

    for j in range(n - 2, -1, -1):
        c[j] = z[j] - mu[j] * c[j + 1]
        b = np.zeros(n - 1)
        d = np.zeros(n - 1)
        for i in range(n - 1):
            b[i] = (y[i + 1] - y[i]) / h[i] - (h[i] * (c[i + 1] + 2 * c[i])) / 3
            d[i] = (c[i + 1] - c[i]) / (3 * h[i])

    return b, c, d

def evaluate_cubic_spline(x, y, b, c, d, query_points):
    results = []
    for query in query_points:
        for i in range(len(x) - 1):
            if x[i] <= query <= x[i + 1]:
                dx = query - x[i]
                result = y[i] + b[i] * dx + c[i] * dx**2 + d[i] * dx**3
                results.append(result)
                break

    return results

def func(x):
    return 1 / (13 + 25 * x**2) - 13

A = -1
B = 1
N = 7

x = np.linspace(A, B, N)
y = func(x)

b, c, d = cubic_spline(x, y)

X = np.linspace(A, B, 1000)

Y = evaluate_cubic_spline(x, y, b, c, d, X)

plt.plot(X, func(X), label='y=1/(13+25x^2) - 13', lw=3.5)
plt.plot(X, Y, 'r--', label='Кубічний сплайн')
plt.plot(x, y, 'yo', label='Точки даних')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()