Linear oscillator solution

We can use Equations (47) to solve the linear oscillator problem for no external forces applied. These equations give for p=2:

The initial conditions for the position are:

The FORTRAN code to solve the linear oscillator differential equation is shown below. The main program contains the routine initial that initializes the variables from the input of the user. These variables are: (the unit system is not specified). Typical values for the constants are given, but of course you can select different values.

k: spring constant, is a positive value and k=1.0 is a typical value. The mass of the oscillator was taken as 1. With k=1.0the period is

: initial amplitude, a value like 1.0 is a good value. Coded as y(1).

: initial velocity, for example 0.0. Coded as y(2).

h: time step. If the period is

, a selection can be 0.1 (see the problems at the end of this section).

tmax: is an integer, and is the number of time steps. Since for k=1 the period is 6.28.., and h=0.1, tmax=250 will give some few periods.

The routine integrate integrates the linear differential equations. The input for the routine are the time step h, the spring constant k, the initial values and contained in the two dimensional array y.

The routine func returns the values of and as ff(1) and ff(2), respectively, from equations (75).

The routine rk4 has as input h; k; values and at the previous time step as y(1), y(2) in array y; and returns $(y_{n+1})^(1)$ and $(y_{n+1})^(2)$ in the array ynp1. This routine is called by integrate many times, 250 times in our example.

Two disk files are produced: linosx.dat and linosv.dat. The first contains x(t) versus time, and the second v(t) versus time. They can be plotted on the same graph.

	PROGRAM linos
C-----------------------------------------------------------|	
C Program for the solution of the linear harmonic oscillator|
C d^2x/dt^2=-k x(t)/m.      Let us call :                   |
C         y1(t)=x(t),        y2(t)=dy1/dt                   |
C         f1(t,y)=y2(t)      f2(t,y)=-k y1(t)/m             |
C using 4th order Runge-Kutta                               |
C h : time step size      k: spring constant    mass m=1    |
C-----------------------------------------------------------|	
        IMPLICIT NONE
        INTEGER tmax
        DOUBLE PRECISION h,y(2),k,ynp1(2)
C       file linosx.dat will contain x vs t data        
        OPEN(9,FILE='linosx.dat',STATUS='NEW')
c       file linosv.dat will contain v vs t data        
        OPEN(8,FILE='linosv.dat',STATUS='NEW')
c       values are initialized        
        CALL initial(h,k,tmax,y)               
        CALL integrate(h,k,tmax,y,ynp1)
        CLOSE(8)
        CLOSE(9)
        STOP
        END
        
        SUBROUTINE initial(h,k,tmax,y)  
c       y(1) initial position     y(2) initial velocity           
c       h: time step     k: spring constant                     
c       computes the period T=2*pi sqrt(m/k)                     
        IMPLICIT NONE
        INTEGER tmax
        DOUBLE PRECISION y(2),h,k,period,pi
        pi=3.14159265358979323846
        WRITE(*,*)'Enter spring constant k (>0.0, like 1.0)'  
        READ(*,*)k
        WRITE(*,*)'you entered the spring constant k=',k
C       for m=1 the period is:
        period=2*pi/sqrt(k)   
        WRITE(*,*)'The (theoretical) period is=',period
        WRITE(*,*)'Enter maximum number of time steps'
        WRITE(*,*)'must be an integer (like 50)'
        READ(*,*)tmax
        WRITE(*,*)'you entered tmax=',tmax
        WRITE(*,*)'Enter the time step (like 0.1)'
        READ(*,*)h
        WRITE(*,*)'you entered time step= ',h    
        WRITE(*,*)'Enter the initial position (like 1.0)'
        READ(*,*)y(1)
        WRITE(*,*)'you entered initial position',y(1)
        WRITE(*,*)'Enter the initial velocity (like 0.0)'
        READ(*,*)y(2)
        WRITE(*,*)'you entered initial velocity',y(2)           
        RETURN
        END
        
	SUBROUTINE func(k,y,ff)
C       computes f1(t,y)  and f2(t,y)   returns ff(1) for f1 and  
C         ff(2) for f2(t,y)                                     
        IMPLICIT NONE
	DOUBLE PRECISION y(2),ff(2),k
	ff(1)=y(2)
	ff(2)=-k*y(1)
	RETURN
	END
	
	SUBROUTINE rk4(h,k,y,ynp1)
C       Fourth order Runge-Kutta routine to solve two ODE         
C         input:  h (time step), y: y(1),y(2)  returns yn+1 at the  
C         new time step.
C       k1,k2,k3,k4 : each dimension 2, in formula:               
C       yn+1 = yn +(k1+2k2+2k3+k4)/6                        
C       ytem: temporary array stores  =y+k1/2 or y+k2/2 or y+k3/2 
        IMPLICIT NONE
	INTEGER i
	DOUBLE PRECISION t,h,k1(2),k2(2),k3(2),k4(2),k
	DOUBLE PRECISION y(2),ynp1(2),ff(2),ytem(2)
	CALL func(k,y,ff)
c       k1's values and ytem	
	DO 10 i=1,2
	  k1(i)=h*ff(i)
	  ytem(i)=y(i)+k1(i)/2.D0
 10	CONTINUE
c       k2's and ytem 
        CALL func(k,ytem,ff)
	DO 20 i=1,2
	  k2(i)=h*ff(i)
	  ytem(i)=y(i)+k2(i)/2.D0
 20     CONTINUE
	CALL func(k,ytem,ff)
	DO 30 i=1,2
	  k3(i)=h*ff(i)
	  ytem(i)=y(i)+k3(i)
 30     CONTINUE	  
	CALL func(k,ytem,ff)
	DO 40 i=1,2
	  k4(i)=h*ff(i)
	  ynp1(i)=y(i)+(k1(i)+2.D0*k2(i)+2.D0*k3(i)+k4(i))/6.D0		
 40     CONTINUE
	RETURN
	END

        SUBROUTINE integrate(h,k,tmax,y,ynp1)
C   integrates the differential equations for the harmonic 
c   oscillator using the 4 order Runge-Kutta routine rk4   

        IMPLICIT NONE
        INTEGER i,tmax        
        DOUBLE PRECISION k,h,y(2),ynp1(2),time
        time=0.D0
        DO 10 i=1,tmax
          CALL rk4(h,k,y,ynp1)  
c         new y is now ynp1
c         send values to diskfile for plotting x vs t
          WRITE(9,*)time,ynp1(1)
          WRITE(8,*)time,ynp1(2)
          y(1)=ynp1(1)
          y(2)=ynp1(2)
          time=time+h          
 10     CONTINUE         
        RETURN
        END

The graphs of the displacement x(t) versus time and the velocity v(t) versus time are shown in next figure. 0 is the equilibrium position for both positive and negative displacements.

Figure 3: The displacement x(t) x(t) and velocity v(t) v(t) for the linear oscillator.