! motion of a real pendulum with damping and a driving force ! with the Euler-Cromer method ! Program to accompany "Computational Physics" by N. Giordano ! Copyright Prentice Hall 1997 program nonlinear option nolet library "sgfunc*", "sglib*" dim theta(50000),omega(50000),t(50000) call initialize(theta,omega,t,length,dt,damping,force,drive_frequency) call calculate(theta,omega,t,length,dt,damping,force,drive_frequency) call display(theta,omega,t,length,dt,damping,force,drive_frequency) end ! initialize variables ! theta() = pendulum angle omega() = pendulum angular velocity ! damping = damping factor, force = amplitude of drive force ! drive_frequency = angular frequency of drive force ! all angles in radians sub initialize(theta(),omega(),t(),length,dt,damping,force,drive_frequency) input prompt "initial pendulum angle (in radians) -> ": theta(1) input prompt "initial angular velocity of pendulum (in radians/s) -> ": omega(1) t(1) = 0 input prompt "length of pendulum (in m) -> ": length input prompt "time step -> ": dt input prompt "damping constant -> ": damping input prompt "amplitude of driving force -> ": force input prompt "driving frequency -> ": drive_frequency end sub ! use the Euler-Cromer method for a damped, nonlinear, driven pendulum sub calculate(theta(),omega(),t(),length,dt,q,drive_force,drive_frequency) i = 0 g = 9.8 period = 2 * pi / sqr(g/length) do i = i + 1 t(i+1) = t(i) + dt omega(i+1) = omega(i) - (g/length) * sin(theta(i)) * dt - q * omega(i) * dt + drive_force * sin(drive_frequency * t(i)) * dt theta(i+1) = theta(i) + omega(i+1) * dt if theta(i+1) > pi then theta(i+1) = theta(i+1) - 2 * pi ! keep theta in if theta(i+1) < -pi then theta(i+1) = theta(i+1) + 2 * pi! range -pi to pi loop until t(i+1) >= 60 ! 10 * period mat redim omega(i+1),theta(i+1),t(i+1) end sub sub display(theta(),omega(),t(),length,dt,damping,force,drive_frequency) set background color "white" set color "black" call setcanvas("white") clear open #1: screen 0,.9,0,.48 call settitle("Angular velocity vs. time") call sethlabel("Time(s)") call setvlabel("Omega") call datagraph(t,omega,1,1,"black") close #1 open #1: screen 0,.9,.52,1 call settitle("Theta vs. time") call sethlabel("Time(s)") call setvlabel("Theta") call datagraph(t,theta,1,1,"black") set cursor 3,10 print "length = "; length set cursor 4,10 print "time step = "; dt close #1 end sub