Чисельне диференцiювання

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd


def f(x):
    return 1 / (x + 13)

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

def diff2(x):
    return 2 / ((x + 13)**3)

def factorial(n):
    if n == 0:
        return 1
    else:
        return n * factorial(n-1)

def polinom_newton(x, y, new_x, h):

    N = len(x)


    xy_df = []
    for i in range(N):
        xy_df.append([x[i], y[i]])
    df1 = pd.DataFrame(data=xy_df, columns=['x', 'y'])

    diff = np.zeros((N, N - 1), dtype=float)
    for i in range(N - 1):
        diff[i][0] = y[i + 1] - y[i]
    for k in range(1, N-1):
        for i in range(N - (k + 1)):
            diff[i][k] = diff[i + 1][k - 1] - diff[i][k - 1]

    df2 = pd.DataFrame(data=diff, columns=[f'delta{i}_y' for i in range(1, N)])
    df = df1.join(df2)
    df.index = ['' for i in range(N)]

    if abs(new_x - x.min()) < abs(new_x - x.max()):
        t = (new_x - x[0]) / h
        u = 1
        s = df.iloc[0, 1]
        for i in range(N-1):
            u = u * (t - i) / factorial(i + 1)
            s = s + u * df.iloc[0, i+2]

        diff_1 = ( df.iloc[0, 2] + df.iloc[0, 3]*(2*t - 1)/2 + df.iloc[0, 4]*(3*t**2 - 6*t + 2)/6 + df.iloc[0, 5]*(4*t**3 - 18*t**2 + 22*t - 6)/24 + df.iloc[0, 6]*(5*t**4 - 40*t**3 + 105*t**2 - 100*t + 24)/120 ) / h
        diff_2 = ( df.iloc[0, 3] + df.iloc[0, 4]*(t - 1) + df.iloc[0, 5]*(12*t**2 - 36*t + 22)/24 + df.iloc[0, 6]*(20*t**3 - 120*t**2 + 210*t - 100)/120 ) / h**2

    else:
        t = (new_x - x[-1]) / h
        u = 1
        s = df.iloc[-1, 1]
        for i in range(N-1):
            u = u * (t + i) / factorial(i + 1)
            s = s + u * df.iloc[N-1-(i+1), i+2]

        diff_1 = (df.iloc[N-2, 2] + df.iloc[N-3, 3]*(2 * t + 1)/2 + df.iloc[N-4, 4]*(3 * t ** 2 + 6 * t + 2)/6 + df.iloc[N-5, 5]*(4*t** 3 + 18*t**2 + 22*t + 6)/12 + df.iloc[N-6, 6] * (5*t**4 + 40*t**3 + 105*t**2 + 100*t + 24) / 120) / h
        diff_2 = (df.iloc[N-3, 3] + df.iloc[N-4, 4]*(t + 1) + df.iloc[N-5, 5]*(6*t**2 + 18*t + 11)/12 + df.iloc[N-6, 6]*(20*t**3 + 120*t**2 + 210*t + 100)/120 ) / (h ** 2)
    # print(f'new_x = {new_x}')
    # print(f'new_y = {s}')


    return (df, s, diff_1, diff_2)

def grafic():
    X = np.linspace(x.min() - 1, x.max() + 1, 100)
    Y = np.zeros(len(X))
    for i in range(len(X)):
        Y[i] = polinom_newton(x, y, X[i], h=h)[1]
    plt.title("f(x)")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.grid()
    plt.plot(X, Y, label="Iнтерполяцiйний полiном Ньютона")
    plt.plot(x, y, 'ro', label='Початкові точки')
    plt.legend()
    plt.show()

def df_with_new_x(X):

    df = pd.DataFrame(columns=['x','y', 'diif1', 'diff2'])
    for i in range(0, len(X)):
        tabl, y1, diff1, diff2 = polinom_newton(x, y, X[i], h=h)
        df.loc[i+1] = [X[i], y1, diff1, diff2]
    return df


def abs_err(X, h):
    return abs(polinom_newton(x, y, new_x=X, h=h)[1] - f(X))

def rel_err(X, h):
    return abs(abs_err(X, h) / f(X))

def prava_rizn(x, h):
    return ( f(x+h) - f(x) ) / h

def center_rizn_diff1(x, h):
    return ( f(x+h) - f(x-h) ) / (2*h)

def center_rizn_diff2(x, h):
    return (f(x + h) - 2*f(x) + f(x-h) ) / (h**2)

def abs_err_diff1(x, h):
    return abs(diff1(x) - center_rizn_diff2(x, h))

def rel_err_diff1(x, h):
    return abs_err_diff1(x, h) / abs(diff1(x))

def abs_err_diff2(x, h):
    return abs(diff2(x) - center_rizn_diff2(x, h))

def rel_err_diff2(x, h):
    return abs_err_diff2(x, h) / abs(diff2(x))

def err(X, h):

    df = pd.DataFrame(columns=[''])
    df.loc['Абсолютна похибка'] = [abs_err(X, h)]
    df.loc['Відносна похибка'] = [rel_err(X, h)]
    df.loc['Права різниця'] = [prava_rizn(X, h)]
    df.loc['Центральна різниця першої похідної'] = [center_rizn_diff1(X, h)]
    df.loc['Центральна різниця другої похідної'] = [center_rizn_diff2(X, h)]
    df.loc['Абсолютна похибка першої похідної'] = [abs_err_diff1(X, h)]
    df.loc['Відносна похибка першої похідної'] = [rel_err_diff1(X, h)]
    df.loc['Абсолютна похибка другої похідної'] = [abs_err_diff2(X, h)]
    df.loc['Відносна похибка другої похідної'] = [rel_err_diff2(X ,h)]
    df.columns = [f'x = {X}; h = {h}']

    return df


h = 0.5
x = np.arange(0, 10.5, h)
y = f(x)

table_x_y = pd.DataFrame(data={'x': x, 'y': y}, index=[i+1 for i in range(len(x))])
print(table_x_y)
x1 = 2
x2 = 8
x3 = 1.05
# grafic()

X = [x1, x2, x3]

new_x = df_with_new_x(X)
print(new_x, '\n')

X1 = err(x1, h)
X2 = err(x2, h)
print(X1, '\n')
print(X2, '\n')

zavd3_1 = err(x1, 1)
zavd3_2 = err(x1, 2)
zavd3_3 = err(x2, 1)
zavd3_4 = err(x2, 2)
print(zavd3_1, '\n')
print(zavd3_2, '\n')
print(zavd3_3, '\n')
print(zavd3_4, '\n')

x = np.arange(0, 5.5, 0.5)
y = f(x)

table_x_y = pd.DataFrame(data={'x': x, 'y': y}, index=[i+1 for i in range(len(x))])
print(table_x_y)
x1 = 2
x2 = 8
x3 = 1.05
# grafic()

X = [x1, x2, x3]

new_x = df_with_new_x(X)
new_x['y_exact'] = [f(x1), f(x2), f(x3)]
print(new_x, '\n')

X1 = err(x1, h)
X2 = err(x2, h)
X3 = err(x3, h)
print(X1, '\n')
print(X2, '\n')