Stirling and Basel interpolation methods

bassel_n_stirling.py

import numpy as np


def get_q(x, cp_x):
    x0 = x[len(x) // 2]
    h = x[1] - x[0]
    return (cp_x - x0) / h


def get_x0_idx(x, cp_x):
    return (np.abs(x - cp_x)).argmin()


def get_window(x, y, idx):
    min_left_idx = 0
    max_right_idx = idx + min(len(x) - idx, len(x))
    bound_range = min(idx - min_left_idx, max_right_idx - idx)
    left_bound = idx - bound_range
    right_bound = idx + bound_range
    return x[left_bound:right_bound + 1], y[left_bound:right_bound + 1]


def execute(x, y, cp_x):
    nearest_x_idx = get_x0_idx(x, cp_x)
    x, y = get_window(x, y, nearest_x_idx)
    q = get_q(x, cp_x)
    if -0.25 < q < 0.25:
        return stirling(y, q)
    elif 0.25 < abs(q) < 0.75:
        return bessel(y, q)
    else:
        print('|q|>0.75, cannot apply any method')


def build_difference_table(y):
    n = len(y)
    deltas = [[0 for j in range(i + 1)] for i in reversed(range(0, n))]
    deltas[0] = y
    eps = 5 * (10 ** -3)
    for i in range(1, n):
        for j in range(n - i):
            deltas[i][j] = deltas[i - 1][j + 1] - deltas[i - 1][j]
            if abs(deltas[i][j]) < eps:
                deltas[i][j] = 0
                print('Критерий остановки выполнился после %i итераций' % i)
                break
        else:
            continue
        return deltas
    raise Exception('Критирий остановки не выполнился, увеличьте точность апроксимации')


def bessel(y, q):
    n = len(y)
    m_idx = n // 2
    deltas = build_difference_table(y)
    f_cp_x = (deltas[0][m_idx] + deltas[0][m_idx + 1]) / 2
    temp1 = 1 / 2
    temp2 = (q - 0.5)
    fact = 1
    d1 = 0
    d2 = 0
    for i in range(1, n):
        if deltas[i][0] == 0:
            break
        m_idx_i = len(deltas[i]) // 2
        if i % 2 != 0:
            f_cp_x += temp2 * deltas[i][m_idx_i] / fact
            temp2 *= (q + d2) * (q - (d2 + 1))
            d2 += 1
        else:
            temp1 *= (q + d1) * (q - (d1 + 1))
            d1 += 1
            f_cp_x += temp1 * (deltas[i][m_idx_i] + deltas[i][m_idx_i - 1]) / fact
        fact *= (i + 1)
    return f_cp_x


def stirling(y, q):
    n = len(y)
    m_idx = n // 2
    deltas = build_difference_table(y)

    # for arr in deltas:
    #     print(arr)
    f_cp_x = y[m_idx]
    q2 = q ** 2
    temp1 = q / 2
    temp2 = q2
    d1 = 1
    d2 = 1
    fact = 1
    for i in range(1, n):
        if deltas[i][0] == 0:
            break
        m_idx_i = len(deltas[i]) // 2
        if i % 2 != 0:
            f_cp_x += temp1 * (deltas[i][m_idx_i] + deltas[i][m_idx_i - 1]) / fact
            temp1 *= (q2 - d1 ** 2)
            d1 += 1
        else:
            f_cp_x += temp2 * deltas[i][m_idx_i] / fact
            temp2 *= (q2 - d2 ** 2)
            d2 += 1
        fact *= (i + 1)
    return f_cp_x

main.py

import numpy as np
from methods import execute
import math


def f(x):
    return math.sin(x ** 2) + math.exp(x) * math.log(x) + math.cosh(x)


def build_table_xy(f):
    x = np.asarray([i * 0.25 for i in range(1, 100)])
    y = [f(p) for p in x]
    return x, y


if __name__ == "__main__":
    data = np.genfromtxt("resources/24.csv", delimiter=',', names=True)
    cp = np.genfromtxt("resources/24-cp.csv", delimiter=',', names=True)
    P = 3.439
    x, y = build_table_xy(f)
    my = execute(x, y, P)
    print('diff between my impl and real value %f' % (my - f(P)))
    print('my own impl %f' % my)
    print('numpy impl %f' % np.interp(P, x, y))
    print('real value %f' % f(P))