Квадратичный интерполяционный сплайн

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

Сплайн строится путем соединения всех точек в массиве между собой кусками параболы:

Обратите внимание, что квадратичный сплайн очень капризный и при выборе неподходящего участка основной функции никакого сглаживания не произойдет!!!

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

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

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

import numpy

def Metod_Gaussa(arr, brr):
    for k in range(arr.shape[0] - 1):
        #поиск строки с максимальным элементом
        max_elem = 0
        str = 0
        for i in range (k, arr.shape[0]):
            if abs(arr[i,k]) > abs(max_elem):
                max_elem = arr[i,k]
                str = i
        #меняем местами строки квадратной матрицы
        change = numpy.repeat(arr[k], 1)
        arr[k], arr[str] = arr[str], change
        #меняем местами элементы вектора-столбца
        change = numpy.repeat(brr[k], 1)
        brr[k], brr[str] = brr[str], change
        #делим полученную строку на max_elem
        arr[k] = arr[k] / max_elem
        brr[k] = brr[k] / max_elem
        #домножаем строку на коэффициенты и вычитаем ее из остальных строк
        for i in range (k + 1, arr.shape[0]):
            factor = arr[i,k]
            arr[i] = arr[i] - arr[k] * factor
            brr[i] = brr[i] - brr[k] * factor

    #находим неизвестные
    arg = [brr[brr.shape[0] - 1] / (arr[arr.shape[0] - 1, arr.shape[0] - 1])]
    for i in range(arr.shape[0] - 2, -1, -1):
        n = brr[i]
        for j in range(len(arg)):
            n = n - arg[j] * arr[i, arr.shape[0] - 1 - j]
        arg.append(n)

    #переворачиваем значения в списке
    otv = []
    for i in reversed(arg): otv.append(i)
    return otv

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

Основной код для левого сплайна:

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

#для решения СЛАУ
from Project1 import Gaussa
import numpy

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

def left_square_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}))
    #ищем константы уравнений
    const = [[] for i in range(N)]
    for i in range(N):
        A = numpy.array([[1, xi[i], xi[i] ** 2], [1, xi[i + 1], xi[i + 1] ** 2], [0, 1, 2 * xi[i]]])
        if i == 0:
            a = numpy.array([yi[i], yi[i + 1], 0])
        else:
            a = numpy.array([yi[i], yi[i + 1], const[i - 1][1] + 2 * const[i - 1][2] * xi[i]])
        for j in Gaussa.Metod_Gaussa(A, a): const[i].append(float(j))
    #вычисление уравнений кусков определённого слева квадратичного сплайна с привязкой к своим отрезкам
    yn = {}
    for i in range(N):
        yn[xi[i], xi[i + 1]] = const[i][0] + const[i][1] * x + const[i][2] * x ** 2
    return yn

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

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

Результат:

Основной код для правого сплайна:

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

#для решения СЛАУ
from Project1 import Gaussa
import numpy

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

def right_square_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}))
    #ищем константы уравнений
    const = [[] for i in range(N)]
    for i in range(N, 0, -1):
        A = numpy.array([[1, xi[i], xi[i] ** 2], [1, xi[i - 1], xi[i - 1] ** 2], [0, 1, 2 * xi[i]]])
        if i == N:
            a = numpy.array([yi[i], yi[i - 1], 0])
        else:
            a = numpy.array([yi[i], yi[i - 1], const[N - 1 - i][1] + 2 * const[N - 1 - i][2] * xi[i]])
        for j in Gaussa.Metod_Gaussa(A, a):
            const[N - i].append(float(j))
    #переворачиваем значения в списке
    const_new = []
    for i in reversed(const): const_new.append(i)
    #вычисление уравнений кусков определённого справа квадратичного сплайна с привязкой к своим отрезкам
    yn = {}
    for i in range(N):
        yn[xi[i], xi[i + 1]] = const_new[i][0] + const_new[i][1] * x + const_new[i][2] * x ** 2
    return yn

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

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

Результат: