SHO-Feynman Newton Method

#Ideal simple 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()


output

#DAMPED  HARMONIC OSCILLATOR
#dx/dt=v(t)
#dv/dt=-w*w*x(t)-c*v(t)
#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*
x=10
v=20
w=2
c=input('enter the damping coefficient[enter the zero for undamped]: ')
t=0
h=0.01
xdat=[]
vdat=[]
time=[]
v=v+h*(-w*w*x-c*v)/2
while(t<=150):
    xdat.append(x)
    vdat.append(v)
    time.append(t)
    x=x+h*v
    v=v+h*(-w*w*x-c*v)
    t=t+h
figure(1)
title('damped harmonic oscillator')
xlabel('time')
ylabel('displacement')
plot(time,xdat)
grid(True)
figure(2)
title('damped harmonic oscillator')
xlabel('time')
ylabel('velocity')
plot(time,vdat)
grid(True)
figure(3)
title('damped harmonic oscillator phase space plot')
xlabel('displacement')
ylabel('velocity')
plot(xdat,vdat)
grid(True)
show()


OUTPUT: