SHARE
TWEET

Untitled

a guest Dec 3rd, 2019 78 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
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
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top