Numeriske pendelsvingninger

Her kan du stille spørsmål vedrørende problemer og oppgaver i matematikk på høyskolenivå. Alle som har kunnskapen er velkommen med et svar. Men, ikke forvent at admin i matematikk.net er spesielt aktive her.

Moderatorer: Vektormannen, espen180, Aleks855, Solar Plexsus, Gustav, Nebuchadnezzar, Janhaa

Svar
Gjest

Hei, ganske langt spørsmål så dette blir for de spesielt interesserte.

Hadde en oppgave å numerisk beregne vinkelutslaget til en vertikal pendel, men fikk det ikke helt til også lurte jeg på om det var noen kloke hoder som kunne hjelpe meg å lokalisere feilen i programmet mitt.
Jeg har selv en liten mistanke om at feilen er knyttet til fysikken og ikke kodingen, men det er bare en antagelse.

Her er et bilde av pendelen (legg merke til at det indre festet har utstrekning r):
Bilde

Oppgaveopplysninger:
Pendelen består av ei tynn metallstang med 4 lodd, hver med masse 50 g, festet nær stangas ende.
Akslingen er ei metallstang med radius r = 8 mm. Avstanden fra akslingens senterakse til de 4
loddenes massesenter er målt til L = 636 mm. Dette er strengt tatt en fysisk pendel, noe vi tar hensyn
til når friksjonskreftene tas i betraktning. (fµ har arm r, mens de andre kreftene har arm L; derfor
må Newtons 2. lov for rotasjon brukes.) Vi har imidlertid ikke tatt hensyn til at pendelstanga har en
viss masse, og at de fire loddene har en viss utstrekning. Det kan derfor tenkes at du i beregningene
dine må variere pendellengden L litt omkring den målte verdien for å få best mulig samsvar mellom
eksperiment og beregninger.

Det ble også oppgitt eksperimentelle verdier for 3 forsøk i tre ulike tekst filer:
track1.txt
t x y theta
0 26.214 -611.505 -87.5
0.01 83.001 -607.841 -82.2
0.02 139.787 -596.85 -76.8
0.03 191.079 -582.196 -71.8
0.04 249.697 -563.878 -66.1
0.05 300.988 -538.232 -60.8
0.06 344.952 -510.755 -56
...
...
...

Kode: Velg alt

from scitools.std import *		#import exp, sin, pi ...

'-----------------------------definies intial functions------------------------------------'

def N(theta, s):
	s = int(s) #for some reason s changes to float on the second interation
	return m*(g*cos(theta)+sign(v[s])**2/float(L))
	
def F(theta, t, s):
	s = int(s)
	return (-m*g*sin(theta)-D*sign(v[s])**2)*L-(mu*N(theta, t))*r #-G_x-F_D-f_mu. the sum of the forces is the momentum
 
def	velocity(theta, t, s): #s parameter carries on from previous foor loop
	theta = theta*pi/180. #converts to radians
	s = int(s) 
	v[s+1] = v[s] + dt*F(theta,t, s)/float(m) #creates next v
	return v[s] #applies previous v

'--------------------------definies inital parameters---------------------------------'

m = 0.2 #4*50g in kg
g = 9.81 #m/s^2
L = 0.636 #m
r = 0.008 #m
D = 0.47 #estimated from sphere like bodies in kg/m
mu = 0.6 #estimated from steel on steel dry kinetic friction coefficient


'-----------------------------reads from file 1----------------------------------------'
time_lenght = 7				#in seconds
n = time_lenght*100			#number of frames					
t = linspace (0, time_lenght, n) #0.01 intervals

theta = [0]*n
v = [0]*n

movie1eksp = []
time1 = []
counter = 0
with open('movie1.txt') as infile:
	for p in range(2):
		infile.readline() #skips header
	for line in infile:
		counter += 1
		lane = line.split()
		movie1eksp.append(float(lane[3]))
		time1.append(float(lane[0]))
		if counter == 1:
			theta[0] = (float(lane[3])+90)*pi/180 #starting angle
		elif counter == 2:
			dt = float(lane[0]) #reads of timestep
			theta[1] = (float(lane[3])+90)*pi/180
			v[0] = L*(abs(theta[0]-theta[1]))/dt #angular velocity times lenght of pendulum

'----------------constructs array for theta and plots for pendulum 1-------------------'
for s in range(n-1):
	dt = t[s+1]-t[s]
	theta[s+1] = theta[s] + (dt*velocity(theta[s],t[s], s)/float(L))*180/pi

movie1eksp = movie1eksp[:700] #only compares the first 7 seconds
time1 = time1[:700] 		#only compares the first 7 seconds
movie1eksp = [x + 90 for x in movie1eksp] #adds on 90 degrees to all values
plot(t, theta, time1, movie1eksp,
	xlabel='time[s]',
	ylabel='displacement[degree]',
	legend=('computed', 'experimental'),
	legend_loc='upper right',
	axis=[0, time_lenght],
	title='Pendulum1')
raw_input()


'-----------------------------reads from file 2----------------------------------------'
time_lenght = 6						#in seconds
n = int(time_lenght*100)			#number of frames					
t = linspace (0, time_lenght, n) #0.01 intervals

theta = [0]*n #
v = [0]*n #v[0] = 0

movie2eksp = []
time2 = []
counter = 0
with open('movie2.txt') as infile:
	for p in range(2):
		infile.readline() #skips header
	for line in infile:
		counter += 1
		lane = line.split()
		movie2eksp.append(float(lane[3])) #creates eksperimental array
		time2.append(float(lane[0]))
		if counter == 1:
			theta[0] = (float(lane[3])+90)*pi/180 #starting angle
		elif counter == 2:
			dt = float(lane[0]) #reads of timestep
			theta[1] = (float(lane[3])+90)*pi/180
			v[0] = L*(abs(theta[0]-theta[1]))/dt #angular velocity times lenght of pendulum
'----------------constructs array for theta and plots for pendulum 2-------------------'
for s in range(n-1):
	dt = t[s+1]-t[s]
	theta[s+1] = theta[s] + (dt*velocity(theta[s],t[s], s)/float(L))*180/pi

movie2eksp = movie2eksp[:600] #only compares the first 6 seconds
time2 = time2[:600] #only compares the first 6 seconds
movie2eksp = [x + 90 for x in movie2eksp] #adds on 90 degrees to all values
plot(t, theta,time2, movie2eksp,
	xlabel='time[s]',
	ylabel='displacement[degree]',
	legend=('computed', 'experimental'),
	legend_loc='upper right',
	axis=[0, time_lenght],
	title='Pendulum2')
raw_input()
'-----------------------------reads from file 3----------------------------------------'
time_lenght = 6						#in seconds
n = int(time_lenght*20)			#number of frames					
t = linspace (0, time_lenght, n) #0.01 intervals

theta = [0]*n #
v = [0]*n #v[0] = 0

movie3eksp = []
time3 = []
counter = 0
with open('movie3.txt') as infile:
	for p in range(2):
		infile.readline() #skips header
	for line in infile:
		counter += 1
		lane = line.split()
		movie3eksp.append(float(lane[3]))
		time3.append(float(lane[0]))
		if counter == 1:
			theta[0] = (float(lane[3])+90)*pi/180 #starting angle
		elif counter == 2:
			dt = float(lane[0]) #reads of timestep
			theta[1] = (float(lane[3])+90)*pi/180
			v[0] = L*(abs(theta[0]-theta[1]))/dt #angular velocity times lenght of pendulum
'----------------constructs array for theta and plots for pendulum 3-------------------'

for s in range(n-1):
	dt = t[s+1]-t[s]
	theta[s+1] = theta[s] + (dt*velocity(theta[s],t[s], s)/float(L))*180/pi

movie3eksp = movie3eksp[:120] #only compares the first 6 seconds
time3 = time3[:120] #only compares the first 6 seconds
movie3eksp = [x + 90 for x in movie3eksp] #adds on 90 degrees to all values
plot(t, theta, t, movie3eksp,
	xlabel='time[s]',
	ylabel='displacement[degree]',
	legend=('computed', 'experimental'),
	legend_loc='upper right',
	axis=[0, time_lenght],
	title='Pendulum3')
raw_input()

Svar