#! /usr/bin/env python3 # -*- coding: UTF-8 -*- # TFY4104/TFY4115 Fysikk 2019. # Cubic spline approximation to roller coaster on the lab # Length unit used: mm import matplotlib.pyplot as plt import numpy as np # (x,y)-data in file splinebane.txt (1 header row): # y = measured height of fixed points of roller coaster # shown in splinebane.jpg data=np.loadtxt('splinebane.txt',skiprows=1) print(data) xfast=data[:,0] yfast=data[:,1] h = 200 b = np.zeros(6) print(b) for i in range(0,5): b[i] = (6/h**2)*(yfast[i+2]-2*yfast[i+1]+yfast[i]) #b[0] = (6/h**2)*(yfast[2]-2*yfast[1]+yfast[0]) #b[1] = (6/h**2)*(yfast[3]-2*yfast[2]+yfast[1]) #b[2] = (6/h**2)*(yfast[4]-2*yfast[3]+yfast[2]) #b[3] = (6/h**2)*(yfast[5]-2*yfast[4]+yfast[3]) #b[4] = (6/h**2)*(yfast[6]-2*yfast[5]+yfast[4]) #b[5] = (6/h**2)*(yfast[7]-2*yfast[6]+yfast[5]) print(b) A = np.zeros((6,6)) A[0,:] = (4,1,0,0,0,0) A[1,:] = (1,4,1,0,0,0) A[2,:] = (0,1,4,1,0,0) A[3,:] = (0,0,1,4,1,0) A[4,:] = (0,0,0,1,4,1) A[5,:] = (0,0,0,0,1,4) z = np.linalg.solve(A,b) print(z) print(A) ddy = np.zeros(8) ddy[0] = 0 ddy[7] = 0 ddy[1] = z[0] ddy[2] = z[1] ddy[3] = z[2] ddy[4] = z[3] ddy[5] = z[4] ddy[6] = z[5] print(ddy) x0=np.linspace(0,200,201) x1=np.linspace(201,400,200) x2=np.linspace(401,600,200) x3=np.linspace(601,800,200) x4=np.linspace(801,1000,200) x5=np.linspace(1001,1200,200) x6=np.linspace(1201,1400,200) x = np.concatenate([x0,x1,x2,x3,x4,x5,x6]) # x = x0+x1+x2+x3+x4+x5+x6 S0 = ddy[0]*(xfast[1]-x0)**3/(6*h) + ddy[1]*(x0-xfast[0])**3/(6*h) \ + (yfast[1]/h-ddy[1]*h/6)*(x0-xfast[0]) \ + (yfast[0]/h-ddy[0]*h/6)*(xfast[1]-x0) S1 = ddy[1]*(xfast[2]-x1)**3/(6*h) + ddy[2]*(x1-xfast[1])**3/(6*h) \ + (yfast[2]/h-ddy[2]*h/6)*(x1-xfast[1]) \ + (yfast[1]/h-ddy[1]*h/6)*(xfast[2]-x1) S2 = ddy[2]*(xfast[3]-x2)**3/(6*h) + ddy[3]*(x2-xfast[2])**3/(6*h) \ + (yfast[3]/h-ddy[3]*h/6)*(x2-xfast[2]) \ + (yfast[2]/h-ddy[2]*h/6)*(xfast[3]-x2) S3 = ddy[3]*(xfast[4]-x3)**3/(6*h) + ddy[4]*(x3-xfast[3])**3/(6*h) \ + (yfast[4]/h-ddy[4]*h/6)*(x3-xfast[3]) \ + (yfast[3]/h-ddy[3]*h/6)*(xfast[4]-x3) S4 = ddy[4]*(xfast[5]-x4)**3/(6*h) + ddy[5]*(x4-xfast[4])**3/(6*h) \ + (yfast[5]/h-ddy[5]*h/6)*(x4-xfast[4]) \ + (yfast[4]/h-ddy[4]*h/6)*(xfast[5]-x4) S5 = ddy[5]*(xfast[6]-x5)**3/(6*h) + ddy[6]*(x5-xfast[5])**3/(6*h) \ + (yfast[6]/h-ddy[6]*h/6)*(x5-xfast[5]) \ + (yfast[5]/h-ddy[5]*h/6)*(xfast[6]-x5) S6 = ddy[6]*(xfast[7]-x6)**3/(6*h) + ddy[7]*(x6-xfast[6])**3/(6*h) \ + (yfast[7]/h-ddy[7]*h/6)*(x6-xfast[6]) \ + (yfast[6]/h-ddy[6]*h/6)*(xfast[7]-x6) S = np.concatenate([S0,S1,S2,S3,S4,S5,S6]) # S = S0+S1+S2+S3+S4+S5+S6 # Her kan analytiske uttrykk for dS0/dx og d2S0/dx2 legges inn, # og tilsvarende for S1, ..., S6 dS0 = -ddy[0]*(xfast[1]-x0)**2/(2*h) \ + ddy[1]*(x0-xfast[0])**2/(2*h) \ + (yfast[1]/h-ddy[1]*h/6) \ - (yfast[0]/h-ddy[0]*h/6) dS1 = -ddy[1]*(xfast[2]-x1)**2/(2*h) \ + ddy[2]*(x1-xfast[1])**2/(2*h) \ + (yfast[2]/h-ddy[2]*h/6) \ - (yfast[1]/h-ddy[1]*h/6) dS2 = -ddy[2]*(xfast[3]-x2)**2/(2*h) \ + ddy[3]*(x2-xfast[2])**2/(2*h) \ + (yfast[3]/h-ddy[3]*h/6) \ - (yfast[2]/h-ddy[2]*h/6) dS3 = -ddy[3]*(xfast[4]-x3)**2/(2*h) \ + ddy[4]*(x3-xfast[3])**2/(2*h) \ + (yfast[4]/h-ddy[4]*h/6) \ - (yfast[3]/h-ddy[3]*h/6) dS4 = -ddy[4]*(xfast[5]-x4)**2/(2*h) \ + ddy[5]*(x4-xfast[4])**2/(2*h) \ + (yfast[5]/h-ddy[5]*h/6) \ - (yfast[4]/h-ddy[4]*h/6) dS5 = -ddy[5]*(xfast[6]-x5)**2/(2*h) \ + ddy[6]*(x5-xfast[5])**2/(2*h) \ + (yfast[6]/h-ddy[6]*h/6) \ - (yfast[5]/h-ddy[5]*h/6) dS6 = -ddy[6]*(xfast[7]-x6)**2/(2*h) \ + ddy[7]*(x6-xfast[6])**2/(2*h) \ + (yfast[7]/h-ddy[7]*h/6) \ - (yfast[6]/h-ddy[6]*h/6) dSanal = np.concatenate([dS0,dS1,dS2,dS3,dS4,dS5,dS6]) d2S0 = ddy[0]*(xfast[1]-x0)/(h) \ + ddy[1]*(x0-xfast[0])/(h) d2S1 = ddy[1]*(xfast[2]-x1)/(h) \ + ddy[2]*(x1-xfast[1])/(h) d2S2 = ddy[2]*(xfast[3]-x2)/(h) \ + ddy[3]*(x2-xfast[2])/(h) d2S3 = ddy[3]*(xfast[4]-x3)/(h) \ + ddy[4]*(x3-xfast[3])/(h) d2S4 = ddy[4]*(xfast[5]-x4)/(h) \ + ddy[5]*(x4-xfast[4])/(h) d2S5 = ddy[5]*(xfast[6]-x5)/(h) \ + ddy[6]*(x5-xfast[5])/(h) d2S6 = ddy[6]*(xfast[7]-x6)/(h) \ + ddy[7]*(x6-xfast[6])/(h) d2Sanal = np.concatenate([d2S0,d2S1,d2S2,d2S3,d2S4,d2S5,d2S6]) #Plotter de ulike funksjonene plt.figure('S(x)',figsize=(12,3)) # plt.plot(x0,S0,x1,S1,x2,S2,x3,S3,x4,S4,x5,S5,x6,S6) plt.plot(x,S) plt.title('S(x)') plt.xlabel('$x$ (mm)',fontsize=20) plt.ylabel('$S(x)$',fontsize=20) plt.ylim(0,350) plt.show() # dS = 1. derivert av S dS = np.gradient(S,x) plt.figure('dS/dx') plt.plot(x,dS,x,dSanal) plt.title('dS/dx') plt.xlabel('$x$ (mm)',fontsize=20) plt.ylabel('$dS/dx$',fontsize=20) # plt.ylim(0,350) plt.show() # Calculate slope angle beta (unit: degrees) #beta = np.arctan(dSanal) beta = np.arctan(dSanal)*180/np.pi #print(len(beta)) #print(beta[0]) # plot angle beta plt.figure('anglebeta') plt.plot(x,beta) plt.title('beta') plt.xlabel('$x$ (mm)',fontsize=20) plt.ylabel(r'$\beta$',fontsize=20) # plt.ylim(0,350) plt.show() # d2S = 2. derivert av S = 1. derivert av dS d2S = np.gradient(dS,x) plt.figure('d2S/dx2') plt.plot(x,d2S,x,d2Sanal) plt.title('d2S/dx2') plt.xlabel('$x$ (mm)',fontsize=20) plt.ylabel('$d^2 S/dx^2$ (1/mm)',fontsize=20) # plt.ylim(0,350) plt.show() # K = krumningen K = d2S/(1+dS**2)**(1.5) Kanal = d2Sanal/(1+dSanal**2)**(1.5) plt.figure('K(x)') plt.plot(x,K,x,Kanal) plt.title('K(x)') plt.xlabel('$x$ (mm)',fontsize=20) plt.ylabel('$K(x)$ (1/mm)',fontsize=20) # plt.ylim(0,350) plt.show() # R = 1/K = krumningsradius R = 1/Kanal plt.figure('R(x)') plt.plot(x,R) plt.title('R(x)') plt.xlabel('$x$ (mm)',fontsize=20) plt.ylabel('$R(x)$ (mm)',fontsize=20) plt.ylim(-5000,5000) plt.show() # Check difference between numerical and analytical # dS and d2S dSreldiff = np.absolute((dS-dSanal)/dS) d2Sreldiff = np.absolute((d2S-d2Sanal)/(d2S+1E-12)) plt.figure('dScheck') plt.plot(x,dSreldiff) plt.title('Checking dS') plt.xlabel('$x$ (mm)',fontsize=20) plt.ylabel('$(dS-dSanal)/dS$',fontsize=20) plt.show() plt.figure('d2Scheck') plt.plot(x,d2Sreldiff) plt.title('Checking d2S') plt.xlabel('$x$ (mm)',fontsize=20) plt.ylabel('$(d2S-d2Sanal)/d2S$',fontsize=20) plt.show()