add more interpolation

Author: Panda-Lewandowski

lab2.1.py

@@ -0,0 +1,110 @@
+"""Интерполяция методом полинома Ньютона"""
+from  math import cos, pi, fabs, factorial
+import numpy as np
+from prettytable import PrettyTable
+
+
+def get_table(up, down, step):
+    xs = []
+    ys = []
+    for x in np.arange(up, down, step):
+        xs.append(x)
+        ys.append(cos(x))
+    return xs, ys
+
+
+def bin_search(lst, x):
+    a = 0
+    b = len(lst) - 1
+    while a < b:
+        m = int((a + b) / 2)
+        if x > lst[m]:
+            a = m + 1
+        else:
+            b = m
+    return b
+
+
+def divided_difference(xs, ys):
+    l = len(xs)
+    if l == 1:
+        return ys[0]
+    else:
+        return (divided_difference(xs[:-1], ys[:-1]) - divided_difference(xs[1:], ys[1:])) / (xs[0] - xs[l - 1])
+
+
+def newton_interpolation(lst_x, lst_y, x, n):
+    i = bin_search(lst_x, x)  # поиск ближайшего
+    if len(lst_x) < i + (n + 1) / 2:  # если не получается взыть выборку с нужным х посередине
+        sample_x = lst_x[len(lst_x) - (n + 1):]
+        sample_y = lst_y[len(lst_y) - (n + 1):]
+    else:
+        if n % 2 != 0:  
+            sample_x = lst_x[i - int((n + 1) / 2): i + int((n + 1) / 2)]
+            sample_y = lst_y[i - int((n + 1) / 2): i + int((n + 1) / 2)]
+        else:
+            sample_x = lst_x[i - int((n + 1) / 2) - 1: i + int((n + 1) / 2)]
+            sample_y = lst_y[i - int((n + 1) / 2) - 1: i + int((n + 1) / 2)]
+    y_x = sample_y[0]
+    for i in range(1, len(sample_x)):
+        k = 1
+        for j in range(i):
+           k *= (x - sample_x[j])
+        dd = divided_difference(sample_x[:i+1], sample_y[:i+1])
+        y_x += (k * dd)
+    return y_x
+
+
+def gen_error(x, y_inter, n, xs):
+    y_real = cos(x)
+    diff = abs(y_real - y_inter)
+    omega = x - xs[n]
+    for i in range(0, n):
+        omega *= (x - xs[i])
+    omega = abs(omega)
+    d = abs(cos(x + (pi * n)/2)) / factorial(n + 1)
+    return diff <= d * omega
+
+
+a = float(input('Введите нижнюю границу значений Х: '))
+b = float(input('Введите верхнюю границу значений Х: '))
+h = float(input('Введите шаг: '))
+
+xss, yss = get_table(a, b, h)
+
+table = PrettyTable()
+table.add_column("X", xss)
+table.add_column("Y", yss)
+
+try:
+    x = float(input('Введите значение х в пределах [{0}, {1}]: '.format(a, b)))
+    if x > b or x < a:
+        x = float(input("X должен лежать в пределах [{0}, {1}]! Повторите ввод х: ".format(a,b)))
+    n = int(input('Введите степень полинома: '))
+    if n < 0:
+        n = int(input('Степень должна быть больше 0!\nПовторите ввод: '))
+
+except ValueError:
+    print("Неверный ввод! Повторите")
+    x = float(input('Введите значение х в пределах [{0}, {1}]: '.format(a, b)))
+    if x > b or x < a:
+        x = float(input("X должен лежать в пределах [{0}, {1}]! Повторите ввод х: ".format(a, b)))
+    n = int(input('Введите степень полинома: '))
+    if n < 0:
+        n = int(input('Степень должна быть больше 0!\nПовторите ввод: '))
+
+print(table)
+in_res = newton_interpolation(xss, yss, x, n)
+r_res = cos(x)
+absolute_error = fabs(in_res - r_res)
+try:
+    ratio_error = absolute_error / fabs(r_res)
+except ZeroDivisionError:
+    ratio_error = None
+    print('Упс')
+print("Результат при интерполяции:  ", in_res)
+print("Реальный результат:  ", r_res)
+print("Абсолютная погрешность: {0:.5e}".format(absolute_error))
+print("Относительная погрешность: {0:.5e}".format(ratio_error))
+print('Погрешность полинома Ньютона:', gen_error(x, in_res, n, xss))
+

lab2.2.py

@@ -0,0 +1,149 @@
+"""Обратная интерполяция"""
+from math import cos, sin, pi, fabs, ceil, log
+from prettytable import PrettyTable
+
+relative_accuracy = 0.001
+
+left_fist = -5
+right_first = 1
+
+left_second = -3
+right_second = 3
+
+
+def f(x, y):
+    return x ** 3 - 15 * y + 4
+
+
+def g(x, y):
+    return cos(x) - y
+
+
+def bisection(a, b, x, eps, f):
+    if f(x, a) == eps:
+        return a
+    if f(x, b) == eps:
+        return b
+    c = (a + b) / 2
+    while abs(b - a) > eps * abs(c) + eps:
+        c = (a + b) / 2
+        if f(x, b) * f(x, c) < 0:
+            a = c
+        else:
+            b = c
+    return (a + b) / 2
+
+
+def for_first_equation(x):
+    return bisection(left_fist, right_first, x, relative_accuracy, f)
+
+
+def for_second_equation(x):
+    return bisection(left_fist, right_first, x, relative_accuracy, g)
+
+
+def create_table(left, right, step, function):
+    x = []
+    y = []
+
+    while left <= right:
+        x.append(left)
+        y.append(function(left))
+        left += step
+
+    return [x, y]
+
+
+def find_beg(x, table, size, near, deg):
+    deg += 1
+    if near == 0 and table[near] > x:
+        return 0
+
+    if near == size - 1 and table[near] < x:
+        return size - deg - 1
+
+    if x <= table[near]:
+        if near < deg / 2:
+            return 0
+        if (size - 1 - near) < (ceil(deg / 2) - 1):
+            return size - deg - 1
+        return near - deg / 2
+
+    if x > table[near]:
+        if near < (ceil(deg / 2) - 1):
+            return 0
+        if size - 1 - deg < deg / 2:
+            return size - 1 - deg
+        return near - (ceil(deg / 2) - 1)
+
+    return 0
+
+
+def nearest_value(x, table, size):
+    if x < table[0]:
+        return 0
+
+    if x > table[size - 1]:
+        return size - 1\
+
+    diff = fabs(x - table[0])
+    first_y = 0
+
+    for i in range(1, size):
+        if fabs(x - table[i]) < diff:
+            first_y = i
+            diff = fabs(x - table[i])
+
+    return first_y
+
+
+def newton_interpolation(x, degree, beginnig, table_x, table_y):
+    result = table_y[beginnig]
+
+    for i in range((beginnig + 1), beginnig + degree):
+        divided = 0
+        for j in range(beginnig, i + 1):
+            difference = 1
+            for k in range(beginnig, i + 1):
+                if (k != j):
+                    difference *= (table_x[j] - table_x[k])
+            divided += (table_y[j] / difference)
+        for k in range(beginnig, i):
+            divided *= (x - table_x[k])
+        result += divided
+
+    return result
+
+
+a = float(input('Введите левый предел: \n'))
+b = float(input('ВВедите правый предел: \n'))
+step = float(input("Вывдете шаг: \n"))
+
+x_column_first, y_column_first = create_table(a, b, step, for_first_equation)
+
+x_column_second, y_column_second = create_table(a, b, step, for_second_equation)
+
+x_result, y_result = [], []
+
+for index in range(len(y_column_first)):
+    y_result.append(x_column_first[index])
+    x_result.append(y_column_second[index] - y_column_first[index])
+
+table_data = PrettyTable()
+table_data.add_column('X', x_result)
+table_data.add_column('Y(X)', y_result)
+print(table_data)
+
+degree = int(input('Enter degree of Newton\'s polynomial: \n'))
+x = 0
+
+for_first_y = nearest_value(x, x_result, len(x_result))
+beginning = find_beg(x, x_result, len(x_result), for_first_y, degree)
+
+result = newton_interpolation(x, degree, int(beginning), x_result, y_result)
+
+print('\n -------- Results -------- \n')
+print('X = ', result)
+y_first = for_first_equation(result)
+y_second = for_second_equation(result)
+print('Y = ', y_first + (y_second - y_first) / 2)