• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

a guest May 22nd, 2019 81 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.
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")
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))])
224.
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))])
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.
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))])
262.
263.     curr_len = len(x)
264.     for i in range(len(din_diffs) - 1):
266.         for j in range(1, curr_len):
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))])
277.     for i in range(len(div_diffs) - 1):
279.         for j in range(1, curr_len):
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)
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")