PROGRAM Damped_harmonic_oscillator
integer::i,n,k
real::t,dt,a,x,v,k1,k2,k3,k4
f(t,v)=-a*v-x
f1(t,x)=v
open(15,FILE="p2out3.dat",STATUS="unknown")
a=1.0
t=0;x=1;v=1;n=300
dt=0.1
write(15,*)t,x
do j=1,n
k1=dt*f(t,v)
k2=dt*f(t+0.5*dt,v+0.5*k1)
k3=dt*f(t+0.5*dt,v+0.5*k2)
k4=dt*f(t+dt,v+k3)
v=v+k1/6+k2/3+k3/3+k4/6
k1=dt*f1(t,x)
k2=dt*f1(t+0.5*dt,x+0.5*k1)
k3=dt*f1(t+0.5*dt,x+0.5*k2)
k4=dt*f1(t+dt,x+k3)
x=x+k1/6+k2/3+k3/3+k4/6
t=t+dt
write(15,*)t,x
enddo
END PROGRAM
integer::i,n,k
real::t,dt,a,x,v,k1,k2,k3,k4
f(t,v)=-a*v-x
f1(t,x)=v
open(15,FILE="p2out3.dat",STATUS="unknown")
a=1.0
t=0;x=1;v=1;n=300
dt=0.1
write(15,*)t,x
do j=1,n
k1=dt*f(t,v)
k2=dt*f(t+0.5*dt,v+0.5*k1)
k3=dt*f(t+0.5*dt,v+0.5*k2)
k4=dt*f(t+dt,v+k3)
v=v+k1/6+k2/3+k3/3+k4/6
k1=dt*f1(t,x)
k2=dt*f1(t+0.5*dt,x+0.5*k1)
k3=dt*f1(t+0.5*dt,x+0.5*k2)
k4=dt*f1(t+dt,x+k3)
x=x+k1/6+k2/3+k3/3+k4/6
t=t+dt
write(15,*)t,x
enddo
END PROGRAM
Comments
Post a Comment