Untitled

a guest
Dec 3rd, 2019
83
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import matplotlib.pyplot as plt
2.
3. x = [15.5, 12.5, 8, 10, 7, 4, 8, 10, 9.5, 14, 18, 17, 22, 25, 19, 24.5, 23, 17, 16, 12.5, 16.5, 21, 17, 11, 5.5, 7.5, 10, 12]
4. y = [32.5, 28.5, 29, 33, 33, 37, 39.5, 38.5, 42, 43.5, 42, 40, 41.5, 37, 35, 33.5, 29.5, 30.5, 32, 19.5, 24.5, 22, 15, 10.5, 2.5, 8, 14.5, 20]
5. t = [k/27 for k in range(len(x))]
6.
7.
8. def gauss(t_, x_):
9. n = len(t_)
10. h = [0 for k in range(n-1)]
11. b = [0 for k in range(n-1)]
12. for i in range(n-1):
13. h[i] = t_[i+1] - t_[i]
14. b[i] = 6*(x_[i+1]-x_[i])/h[i]
15. u = [0 for k in range(n-1)]
16. v = [0 for k in range(n-1)]
17. u[0] = 2*(h[0] + h[1])
18. v[0] = b[1] - b[0]
19. for i in range(1,n-1):
20. u[i] = 2*(h[i-1] + h[i]) - h[i-1]*h[i-1]/u[i-1]
21. v[i] = b[i] - b[i-1] - h[i-1]*v[i-1]/u[i-1]
22. z = [0 for k in range(n)]
23. for i in range(n-2,0,-1):
24. z[i] = (v[i] - h[i]*z[i+1])/u[i]
25. return z
26.
27.
28. z_x = gauss(t,x)
29. z_y = gauss(t,y)
30.
31.
32. def s(x_,y_,z_):
33. i = 0
34. n = len(t)
35. for k in range(n-1,0,-1):
36. if x_-k >= 0:
37. i = k
38. h_i = t[i+1]-t[i]
39. s_i = z_[i]/(6*h_i)*(t[i+1]-x_)**3 + z_[i+1]/(6*h_i)*(x_-t[i])**3 \
40. + (y_[i+1]/h_i - z_[i+1]*h_i/6)*(x_- t[i])+ (y_[i]/h_i - z_[i]*h_i/6)*(t[i+1] - x_)
41. return s_i
42.
43. M=1000
44. s_u_x = [s(k/M,x,z_x) for k in range(M)]
45. s_u_y = [s(k/M,y,z_y) for k in range(M)]
46. plt.plot(s_u_x, s_u_y)
47. plt.show()
RAW Paste Data