Search This Blog

Tuesday, June 21, 2011

Damped Harmonic oscillator


#Damped harmonic oscillator
#dx/dt=v
#dv/dt=-w*w*x-b*v
#Feynman Newton algorithm
from pylab import*
from math import*
x0=10
v0=20
w=2
b=0.05
t=0 #initial time
h=0.01 #time step size

xdat=[]
vdat=[]#velocity data store here
time=[]#time stored here

v0=v0+h*(-w*w*x0-b*v0)/2

while(t<=250):
   
    x1=x0+h*v0
    v1=v0+h*(-w*w*x1-b*v0)
   
    xdat.append(x1)
    vdat.append(v1)
    time.append(t)
    x0=x1
    v0=v1
    t=t+h
   
figure(1)
title("Damped harmonic oscillator motion")
xlabel(" Time")
ylabel("x")
plot(time,xdat)
grid(True)
figure(2)
title("Damped harmonic oscillator motion")
xlabel("Time")
ylabel("Velocity")
plot(time,vdat)
grid(True)
figure(3)
title("Damped harmonic oscillator motion")
xlabel("x")
ylabel("velocity")
plot(xdat,vdat)
grid(True)

show()

Projectile motion


#projectile motion
from pylab import*
from math import*
ax=0
ay=-9.8
theta=30
v=20
vx0=v*cos(theta*pi/180)
vy0=v*sin(theta*pi/180)
x0=0
y0=0
t=0 #initial time
h=0.001 #time step size

xdat=[x0]
ydat=[y0]#position data stored here
vdat=[v]#velocity data store here
time=[t]#time stored here

while(y0>=0):
    x1=x0+h*vx0
    y1=y0+h*vy0
    vx1=vx0+h*ax
    vy1=vy0+h*ay
    xdat.append(x1)
    ydat.append(y1)
    vdat.append(sqrt(vx1*vx1+vy1*vy1))
    time.append(t)

    x0=x1
    y0=y1
    vx0=vx1
    vy0=vy1
    t=t+h

figure(1)
title("Motion of Projectile")
xlabel(" x")
ylabel("y")
plot(xdat,ydat,"r-")


figure(2)
title("Motion of projectile")
xlabel(" Time")
ylabel("Velocity")
plot(time,vdat,"k-")
show()


   

Harmonic oscillator


#Harmonic oscillator
#dx/dt=v
#dv/dt=-w*w*x
#Feynman-Newton method
#x(t+h)=x(t)+h*v(t+h/2)
#v(t+h/2)=v(t-h/2)+h*a(t)
#v(1/2)=v0+(h/2)a0
from pylab import*
from math import*
x0=5
v0=5
w=1
t=0 #initial time
h=0.001 #time step size

xdat=[]
vdat=[]#velocity data store here
time=[]#time stored here
v0=v0+h*(-w*w*x0)/2
while(t<=100):
    x1=x0+h*v0
    v1=v0+h*(-w*w*x1)
   
   
    xdat.append(x1)
    vdat.append(v1)
    time.append(t)
    x0=x1
    v0=v1
    t=t+h
   
figure(1)
title("Harmonic oscillator motion")
xlabel(" Time")
ylabel("x")
plot(time,xdat)
grid(True)

figure(2)
title("Harmonic oscillator motion")
xlabel("Time")
ylabel("Velocity")
plot(time,vdat)
grid(True)

figure(3)
title("Harmonic oscillator motion")
xlabel("x")
ylabel("Velocity")
plot(xdat,vdat)
grid(True)

show()

Advanced Python programming notes by Purushothaman sir

Simulation falling through viscous medium


#Falling of a body through viscous medium
#Loss of weight   V*d1*g-V*d2*g=V*d1*g(1-d2/d1)=m*a where a=g(1-d2/d1)
#Drag in the case of ball 6*pi*eta*a*v=k*v
#dy/dt=v
#dv/dt=a-kv
from pylab import*
from math import*
d1=2000.0
d2=750.0
g=9.8

height=2000
y0=0
v0=0
a=(1.0- d2/d1)*g #including buoyant force
k=0.2 #drag coefficient
h=0.001
t=0
y=[y0]
v=[v0]
time=[t]
 
while height-y0>0:
    v1=v0+h*(a-k*v0)
    y1=y0+h*v0
 
    t=t+h
    y.append(height-y1)
    v.append(fabs(v1))
    time.append(t)
    y0=y1
    v0=v1

 
figure(1)
title("Fall through viscous medium")
xlabel(" Time")
ylabel("Position")
plot(time,y)


figure(2)
title("Fall through viscous medium")
xlabel(" Time")
ylabel("Velocity")
plot(time,v)
show()


 

Friday, June 10, 2011

Simulation free fall of a body

#Free fall of a body
#We are solving two first order differential equations by Euler method
#1. dy/dt=v initial condition y0
#2. dv/dt=a initial condition v0


from pylab import*
from math import*

a=-9.8 #acceleration 
v0=0.0 #initial velocity
y0=100 #initial height
t=0 #initial time
h=0.01 #time step size
ydat=[y0]#position data stored here
vdat=[v0]#velocity data stored here
time=[t]#time stored here

#Implementation Euler method
while(y0>0):
    y1=y0+h*v0
    v1=v0+h*a
    t=t+h
    ydat.append(y1)
    vdat.append(fabs(v1))
    time.append(t)
    y0=y1
    v0=v1
   
figure(1)
title("Free fall of a body")
xlabel(" Time")
ylabel("Position")
plot(time,ydat)


figure(2)
title("Free fall of a body")
xlabel(" Time")
ylabel("Velocity")
plot(time,vdat)
show()