Fortran Simple Harmonic Oscillator Problem

[Dylicious]

Hello fellow computer physics nerds,

I'm trying to write a program to plot the positions of the three particles connected by two springs (one dimensional) in Fortran 90. I have a main program block and a module that calls a PGPLOT.

My problem is that the positions of the second and third particle overlap at some points as my code stands and I don't know why!

Here is my code please help meeeeeeeeeeeee

program springparticles
use ploty
implicit none

! Variable Declaration
real :: position_p1, position_p2, position_p3, velocity_p1, velocity_p2, velocity_p3, mass_p1, mass_p2, mass_p3, force_p1, force_p2, force_p3, dx, r
real, parameter :: spring_constant=0.1 , damping_constant=0.01 , viscous_drag=0.1
integer :: n
real, dimension(1000)::ap, bp, cp, tp
real :: t,dt
write (*,*)'Please enter the mass of the first particle: '
read *, mass_p1
write (*,*)'Second Particle: '
read *, mass_p2
write (*,*)'Third Particle: '
read *, mass_p3

position_p1=1.0
position_p2=2.0
position_p3=3.0

t= 0
dt = 0.1

dx = 0

force_p1=0
force_p2=0
force_p3=0

velocity_p1=0
velocity_p2=0
velocity_p3=0

print *,'Particle One |||||| Particle Two |||||| Particle Three'

do n = 1,1000

t = t +dt
dx = dx + 0.01


! force_p1 = force_p1-viscous_drag*velocity_p1
! force_p2 = force_p2-viscous_drag*velocity_p2

r= 0.5

force_p1 = spring_constant*(position_p2-position_p1-r) - damping_constant*(velocity_p1-velocity_p2)

force_p2 = - spring_constant*(position_p2-position_p1-r) - damping_constant*(velocity_p1-velocity_p2) + spring_constant*(position_p3-position_p2-r) - damping_constant*(velocity_p2-velocity_p3)

force_p3 = -spring_constant*(position_p3-position_p2-r) - damping_constant*(velocity_p2-velocity_p3)



velocity_p1 = velocity_p1+force_p1/mass_p1*dt
velocity_p2 = velocity_p2+force_p2/mass_p1*dt
velocity_p3 = velocity_p3+force_p3/mass_p3*dt

position_p1 = position_p1+velocity_p1*dt
position_p2 = position_p2+velocity_p2*dt
position_p3 = position_p3+velocity_p3*dt


! velocity_p1 = velocity_p1+dx/dt
! velocity_p2 = velocity_p2+force_p2/mass_p1p2*dt



print *,position_p1,'|| ',position_p2,'|| ',position_p3

tp(n) = t
ap(n)= position_p1
bp(n)= position_p2
cp(n)= position_p3

end do



call plotfunction(tp, ap, bp, cp)
end program springparticles


module ploty
implicit none
contains
!-----------------------------------------------------------------------
! This module uses pgplot to plot the positions of the particles
!-----------------------------------------------------------------------
subroutine plotfunction(t, a, b, c)
implicit none
integer, parameter :: d = 100
real :: x(100), y(100)
real, intent(in) :: t(:), a(:), b(:), c(:)
integer :: pgopen

! Open a plot window
IF (PGOPEN('/XWINDOW') .LE. 0) STOP


! Set-up plot axes
call PGENV(minval(t),maxval(t),min(minval(a),minval(b)),max(maxval(c),maxval(b)),0,0)
call PGLAB('Time', 'Distance from Equilibrium', 'Positions of Particles')


! Change plot colour to colour 6 ()
! Compute the function at the points
! x(d) = t(:)
! y(d) = a(:)

! write(6,*) t(::10)
! write(6,*) a(::10)
! Plot the curve
call PGSCI(6)
call PGLINE(size(a),t,a)
call PGSCI(4)
call PGLINE(size(b),t,b)
call PGSCI(7)
call PGLINE(size(c),t,c)

! Pause and then close plot window
call PGCLOS
end subroutine plotfunction
!-----------------------------------------------------------------------
end module ploty

 

[Mark44]

My problem is that the positions of the second and third particle overlap at some points as my code stands and I don't know why!

Is there some reason that you think they shouldn't overlap? I don't know what values you're putting in for the masses, and with some values for these, plus the values of the spring and damping constants and viscous drag, it might be that two or more of the particles collide from time to time.

 

[Dylicious]

I just figured that conceptually the spring would compress to a point and never allow the particles to move past one another, mathematically I see how its possible but doesn't that just mean I'm modelling the physical situation inaccurately?

 

[Mark44]

I just figured that conceptually the spring would compress to a point and never allow the particles to move past one another, mathematically I see how its possible but doesn't that just mean I'm modelling the physical situation inaccurately?

I don't recall seeing anything that sets a minimum length for the springs. Your program doesn't understand the physical characteristics of springs such as minimum length when it is compressed, or maximum length at which the spring deforms or breaks when it is stretched. If you want your program to model that behavior, you need to have code that does this.