SHARE
TWEET

Untitled

a guest May 22nd, 2019 79 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import copy
  2. import numpy as np
  3. import sympy as sp
  4. import matplotlib.pyplot as plt
  5. from prettytable import PrettyTable
  6.  
  7. def get_interpolation_lagrange_polynomial(list_x, list_y):
  8.  
  9.     x = sp.symbols('x')
  10.  
  11.     expr = 1
  12.     for x_i in list_x:
  13.         expr = expr * (x - x_i)
  14.  
  15.     func_l = 0
  16.     for x_i, y_i in zip(list_x, list_y):
  17.         func_l = func_l + y_i * (expr / (x - x_i)) / (expr / (x - x_i)).subs(x, x_i)
  18.  
  19.     return sp.lambdify(x, func_l, "numpy")
  20.  
  21.  
  22. def table_finite_difference(x, y):
  23.  
  24.     n = len(y)
  25.     diff = [x, y]
  26.     for diff_num in range(n - 1):
  27.         diff.append(['' for _ in range(n)])
  28.         for i in range(1, n - diff_num):
  29.             diff[diff_num + 2][i - 1] = diff[diff_num + 1][i] - diff[diff_num + 1][i - 1]
  30.  
  31.     return diff
  32.  
  33.  
  34. def divided_difference(X, Y):
  35.  
  36.     # dif_diff(x0..xn) = sum(f(xj) / П(i!=j)(xj - xi), j = 1..n)
  37.     x = sp.symbols('x')
  38.  
  39.     expr = 1
  40.     for xi in X:
  41.         expr = expr * (x - xi)
  42.  
  43.     diff = 0
  44.     for xi, yi in zip(X, Y):
  45.         diff += yi / (expr / (x - xi)).subs(x, xi)
  46.  
  47.     return diff
  48.  
  49.  
  50. def table_divided_difference(x, y):
  51.  
  52.     n = len(y)
  53.     diff = [x, y]
  54.     for order in range(1, n):
  55.         diff.append(['' for _ in range(n)])  # empty list with n elements
  56.         for i in range(n - order):
  57.             diff[order + 1][i] = divided_difference(x[i:i + order + 1], y[i:i + order + 1])  # + 1 last element
  58.  
  59.     return diff
  60.  
  61.  
  62. def newton_polynomial(div_diffs):
  63.  
  64.     x = sp.symbols('x')
  65.  
  66.     expr = 1
  67.     n = len(div_diffs)  # len(x) + 2
  68.     N = 0
  69.     for i in range(1, n):  # div_diffs[0] = x_vector
  70.         N += div_diffs[i][0] * expr
  71.         expr = expr * (x - div_diffs[0][i - 1])
  72.  
  73.     return sp.lambdify(x, N, "numpy")
  74.  
  75.  
  76. ############################################################################################################
  77.  
  78. def linear_interpolation(list_x, list_y):
  79.  
  80.     n = len(list_x)
  81.     a = np.empty(n - 1)
  82.     b = np.empty(n - 1)
  83.  
  84.     for i in range(n - 1):
  85.  
  86.         # F(x) = a_i * x + b_i
  87.         a[i] = (list_y[i + 1] - list_y[i]) / (list_x[i + 1] - list_x[i])
  88.         b[i] = list_y[i] - a[i] * list_x[i]
  89.  
  90.         space_x = np.linspace(list_x[i], list_x[i + 1])
  91.         plt.plot(space_x, a[i] * space_x +           b[i])
  92.  
  93.     plt.scatter(list_x, list_y)
  94.     plt.xlabel("Axis: X")
  95.     plt.ylabel("Axis: Y")
  96.     plt.title("Linear Interpolation")
  97.     plt.show()
  98.  
  99.     return a, b
  100.  
  101.  
  102. def quadratic_interpolation(list_x, list_y):
  103.  
  104.     n = len(list_x)
  105.     a = np.empty(n - 1)
  106.     b = np.empty(n - 1)
  107.     c = np.empty(n - 1)
  108.  
  109.     for i in range(n - 2):
  110.  
  111.         # F(x) = a_i * x + b_i + c_i * x^2
  112.         a[i] = (list_y[i + 2] - list_y[i]) / (list_x[i + 2] - list_x[i]) / (list_x[i + 2] - list_x[i + 1])  \
  113.              - (list_y[i + 1] - list_y[i]) / (list_x[i + 1] - list_x[i]) / (list_x[i + 2] - list_x[i + 1])
  114.  
  115.         b[i] = (list_y[i + 1] - list_y[i]) / (list_x[i + 1] - list_x[i]) - a[i] * (list_x[i + 1] + list_x[i])
  116.  
  117.         c[i] = list_y[i] - b[i] * list_x[i] - a[i] * list_x[i] ** 2
  118.  
  119.     for i in range(0, n - 2, 2):
  120.         space_x = np.linspace(list_x[i], list_x[i + 2])
  121.         plt.plot(space_x, a[i] * space_x ** 2 + b[i] * space_x + c[i])
  122.  
  123.     curr_len = len(list_x)
  124.     if curr_len % 2 == 0:
  125.         i = curr_len - 2
  126.         ax = (list_y[i + 1] - list_y[i]) / (list_x[i + 1] - list_x[i])
  127.         bx = list_y[i] - ax * list_x[i]
  128.  
  129.         space_x = np.linspace(list_x[i], list_x[i + 1])
  130.         plt.plot(space_x, ax * space_x + bx)
  131.  
  132.  
  133.     plt.scatter(list_x, list_y)
  134.     plt.xlabel("Axis: X")
  135.     plt.ylabel("Axis: Y")
  136.     plt.title("Quadratic Interpolation")
  137.     plt.show()
  138.  
  139.     return a, b, c
  140.  
  141.  
  142. def cubic_interpolation(list_x, list_y):
  143.  
  144.     n = len(list_x)
  145.  
  146.     h_i = np.array([list_x[i] - list_x[i - 1] for i in range(1, n)])
  147.     l_i = np.array([(list_y[i] - list_y[i - 1]) / h_i[i - 1] for i in range(1, n)])
  148.     delta_i = np.empty(n - 2, float)
  149.     lambda_i = np.empty(n - 2, float)
  150.  
  151.     delta_i[0] = -0.5 * h_i[1] / (h_i[0] + h_i[1])
  152.     lambda_i[0] = 1.5 * (l_i[1] - l_i[0]) / (h_i[0] + h_i[1])
  153.  
  154.     for i in range(1, n - 2):
  155.         delta_i[i] = - h_i[i + 1] / (2 * h_i[i] + 2 * h_i[i + 1] + h_i[i] * delta_i[i - 1])
  156.         lambda_i[i] = (2 * l_i[i + 1] - 3 * l_i[i] - h_i[i] * lambda_i[i - 1]) / \
  157.                       (2 * h_i[i] + 2 * h_i[i + 1] + h_i[i] * delta_i[i - 1])
  158.  
  159.     a = np.array(copy.copy(list_y)[1:])
  160.     b = np.empty(n - 1)
  161.     c = np.empty(n - 1)
  162.     d = np.empty(n - 1)
  163.     c[n - 2] = 0
  164.  
  165.     for i in range(n - 3, -1, -1):
  166.         c[i] = delta_i[i] * c[i + 1] + lambda_i[i]
  167.     for i in range(n - 2, -1, -1):
  168.         b[i] = l_i[i] + 2 / 3 * c[i] * h_i[i] + 1 / 3 * h_i[i] * c[i - 1]
  169.         d[i] = (c[i] - c[i - 1]) / (3 * h_i[i])
  170.  
  171.     for i in range(n - 1):
  172.         space_x = np.linspace(list_x[i], list_x[i + 1])
  173.         plt.plot(space_x, a[i] + (b[i] + (c[i] + d[i] * (space_x - list_x[i + 1])) * (space_x - list_x[i + 1])) *
  174.             (space_x - list_x[i + 1]))
  175.  
  176.     plt.scatter(list_x, list_y)
  177.     plt.xlabel("Axis: X")
  178.     plt.ylabel("Axis: Y")
  179.     plt.title("Cubic Interpolation")
  180.     plt.show()
  181.  
  182.     return a, b, c, d
  183.  
  184.  
  185. def linear_spline_function(x_value, x_points, a, b):
  186.     interval_index = len(x_points) - 2
  187.     for i in range(1, len(x_points)):
  188.         if x_value < x_points[i]:
  189.             interval_index = i - 1
  190.             break
  191.  
  192.     return a[interval_index] * x_value + b[interval_index]
  193.  
  194.  
  195. def quadratic_spline_function(x_value, x_points, a, b, c):
  196.     interval_index = len(x_points) - 3
  197.     for i in range(1, len(x_points) - 1):
  198.         if x_value < x_points[i]:
  199.             interval_index = i - 1
  200.             break
  201.  
  202.     return c[interval_index] + (b[interval_index] + a[interval_index] * x_value) * x_value
  203.  
  204.  
  205. def cubic_spline_function(x_value, x_points, a, b, c, d):
  206.     ii = len(x_points) - 2  # ii == interval_index
  207.     for i in range(1, len(x_points)):
  208.         if x_value < x_points[i]:
  209.             ii = i - 1
  210.             break
  211.  
  212.     return a[ii] + (b[ii] + (c[ii] + d[ii] * (x_value - x_points[ii + 1])) * (x_value - x_points[ii + 1])) * (
  213.             x_value - x_points[ii + 1])
  214.  
  215. def main():
  216.  
  217.     x = [0.357, 0.871, 1.567, 2.032, 2.628, 3.228]
  218.     y = [0.548, 1.012, 1.159, 0.694, -0.503, 0.703]
  219.  
  220.     table = PrettyTable()
  221.     table.add_column("#", [i for i in range(len(x))])
  222.     table.add_column("x", x)
  223.     table.add_column("y", y)
  224.  
  225. #region Task_1
  226.  
  227.     print("\nTable Values:")
  228.     print(table)
  229.  
  230.     func_l = get_interpolation_lagrange_polynomial(x, y)
  231.     table_for_func_l = PrettyTable()
  232.     table_for_func_l.add_column("#", [i for i in range(len(x))])
  233.     table_for_func_l.add_column("x", x)
  234.     table_for_func_l.add_column("y", [func_l(x_i) for x_i in x])
  235.  
  236.     print("\nLagrangian polynomial:")
  237.     print(table_for_func_l)
  238.  
  239.     print("\nValue of L(x1 + x2):")
  240.     print(func_l(x[0] + x[1]))
  241.  
  242.     func_l_range = np.linspace(0, x[-1])
  243.     plt.plot(func_l_range, func_l(func_l_range))
  244.     plt.scatter(x, y)
  245.     plt.xlabel("Axis: X")
  246.     plt.ylabel("Axis: Y")
  247.     plt.title("Lagrangian polynomial")
  248.     plt.grid()
  249.     plt.show()
  250.  
  251. #endregion
  252.  
  253. #region Task_2
  254.  
  255.     din_diffs = table_finite_difference(x, y)
  256.  
  257.     print("\nFinite difference:")
  258.     finite_table = PrettyTable()
  259.     finite_table.add_column("#", [i for i in range(len(x))])
  260.     finite_table.add_column("x", x)
  261.     finite_table.add_column("y", y)
  262.  
  263.     curr_len = len(x)
  264.     for i in range(len(din_diffs) - 1):
  265.         add = ["None"]
  266.         for j in range(1, curr_len):
  267.             add.append(din_diffs[j][i])
  268.         finite_table.add_column("del_" + str(i + 1), add)
  269.     print(finite_table)
  270.  
  271.     print("\nDivide difference:")
  272.     div_diffs = table_divided_difference(x, y)
  273.     div_table = PrettyTable()
  274.     div_table.add_column("#", [i for i in range(len(x))])
  275.     div_table.add_column("x", x)
  276.     div_table.add_column("y", y)
  277.     for i in range(len(div_diffs) - 1):
  278.         add = ["None"]
  279.         for j in range(1, curr_len):
  280.             add.append(div_diffs[j][i])
  281.         div_table.add_column("del_" + str(i + 1), add)
  282.     print(div_table)
  283.  
  284.  
  285.  
  286.     f_N = newton_polynomial(div_diffs)
  287.     print("\nValue of N(x1 + x2):")
  288.     print(f_N(x[0] + x[1]))
  289.  
  290.     range_x = np.linspace(0, x[-1])
  291.     plt.plot(range_x, f_N(range_x))
  292.     plt.scatter(x, y)
  293.     plt.xlabel("Axis X")
  294.     plt.ylabel("Axis Y")
  295.     plt.title("Newton polynomial")
  296.     plt.show()
  297.  
  298. #endregion
  299.  
  300.     linear_coefficients = linear_interpolation(x, y)
  301.     quadratic_coefficients = quadratic_interpolation(x, y)
  302.     cubic_coefficients = cubic_interpolation(x, y)
  303.  
  304.     space = np.linspace(0, x[-1])
  305.     plt.xlabel("Axis: X")
  306.     plt.ylabel("Axis: Y")
  307.     plt.plot(space, func_l(space), label="Lagrange")
  308.     plt.plot(space, f_N(space), label="Newton")
  309.     plt.plot(space, [linear_spline_function(x_value, x, *linear_coefficients) for x_value in space], label="Linear")
  310.     plt.plot(space, [quadratic_spline_function(x_value, x, *quadratic_coefficients) for x_value in space],label="Quadratic")
  311.     plt.plot(space, [cubic_spline_function(x_value, x, *cubic_coefficients) for x_value in space], label="Cubic")
  312.     plt.scatter(x, y)
  313.     plt.legend()
  314.     plt.show()
  315.  
  316. if __name__ == '__main__':
  317.     main()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top