Кубический интерполяционный сплайн

Для проверки работоспособности программы требуется для функции f(x) = (x - 1)⁶, задав отрезок [a, b] из крайних точек, выбирая промежуточные точки с шагом ℎ = (b - a) / N, N задано, построить сплайн.

Теорию рекомендую почитать здесь.

Прежде чем запускать основной код, нужно сформировать код для метода прогонки и проинициализировать его.

Вспомогательный код для решения СЛАУ методом прогонки (подробнее в этой статье):

def Metod_Progonki(A,B):
    #Известные константы
    k1 = -A[0,1]
    m1 = B[0]
    k2 = -A[A.shape[0] - 1, A.shape[1] - 2]
    m2 = B[B.shape[0] - 1]
    alfa = k1
    beta = m1
    #Поиск альф и бет
    c = 2
    a = 0
    b = 1
    alf = [alfa]
    bet = [beta]
    for i in range(1, A.shape[0] - 1):
        beta = (B[i] - A[i,a] * beta) / (A[i,a] * alfa + A[i,b])
        alfa = -A[i,c] / (A[i,a] * alfa + A[i,b])
        a += 1
        b += 1
        c += 1
        alf.append(alfa)
        bet.append(beta)
    #расчет игриков
    y = (k2 * beta + m2) / (1 - k2 * alfa)
    otv = [y]
    for i in range(len(alf) - 1, -1, -1):
        y = alf[i] * y + bet[i]
        otv.append(y)
    #переворачиваем значения в списке
    otvet = []
    for i in reversed(otv): otvet.append(i)
    return otvet

Этот код следует добавить в отдельный файл в ту же папку PyCharm, что и основной код:

Основной код:

#библиотека для параметризации функции и построения графиков
from sympy import *

#для решения трехдиагональной СЛАУ
from Project1 import Progonka
import numpy

x = symbols('x')
y = (x - 1) ** 6 #исходная функция
a, b = 0.5, 1.5 #отрезок из крайних точек
N = 10 #число разбиений

def cube_interp(y, a, b, N):
    #получение списков координат опорных точек сплайна
    global xi, yi
    xi, yi = [], []
    h = (b - a) / N #шаг
    while a < b:
        xi.append(a)
        yi.append(y.evalf(subs={'x':a}))
        a += h
    if round(a - h, 5) != b:
        xi.append(b)
        yi.append(y.evalf(subs={'x': b}))
    #ищем константы ci уравнений
    A = numpy.zeros([N - 1, N - 1])
    A[0, 0] = 2 * (xi[2] - xi[0])
    A[0, 1] = xi[2] - xi[1]
    A[N - 2, N - 3] = xi[N - 1] - xi[N - 2]
    A[N - 2, N - 2] = 2 * (xi[N] - xi[N - 2])
    for i in range(2, N - 1):
        A[i - 1, i - 2] = xi[i] - xi[i - 1]
        A[i - 1, i - 1] = 2 * (xi[i + 1] - xi[i - 1])
        A[i - 1, i] = xi[i + 1] - xi[i]

    a = numpy.zeros([N - 1])
    for i in range(1, N):
        a[i - 1] = 6 * ((yi[i + 1] - yi[i]) / (xi[i + 1] - xi[i]) - (yi[i] - yi[i - 1]) / (xi[i] - xi[i - 1]))

    #с помощью метода прогонки находим константы ci
    ci = [0]
    for j in Progonka.Metod_Progonki(A, a):
        ci.append(float(j))
    ci.append(0)

    #находим константы di
    di = [0]
    for i in range(1, N + 1):
        di.append((ci[i] - ci[i - 1]) / (xi[i] - xi[i - 1]))

    #находим константы bi
    bi = [0]
    for i in range(1, N + 1):
        hi = xi[i] - xi[i - 1]
        bi.append(ci[i] * hi / 2 - di[i] * hi ** 2 / 6 + (yi[i] - yi[i - 1]) / hi)

    #находим константы ai
    ai = []
    for i in yi: ai.append(i)

    #вычисление уравнений кусков кубического сплайна с привязкой к своим отрезкам
    yn = {}
    for i in range(1, N + 1):
        yn[xi[i - 1], xi[i]] = ai[i] + bi[i] * (x - xi[i]) + ci[i] * (x - xi[i]) ** 2 / 2 + di[i] * (x - xi[i]) ** 3 / 6
    return yn

#вызов функции и вывод на экран
print(cube_interp(y, a, b, N))

#построение функции и сплайна в одних осях
p = plot(y, (x, a, b), show=False)
for i, j in cube_interp(y, a, b, N):
    p.extend(plot(cube_interp(y, a, b, N)[i, j], (x, i, j), line_color='r', show=False))
p.show()

Результат: