pendulv.m

% pendul - Program to compute the motion of a simple pendulum
% using the Verlet method
clear;  help pendulv     % Clear the memory and print header
theta = input('Enter initial angle (in degrees) - ');
theta = theta*pi/180;    % Convert angle to radians
omega = 0;               % Set the initial velocity
tau = input('Enter time step - ');
g_over_L = 1;            % The constant g/L
time_old = -1;           % Fake value (see below)
irev = 0;                % Used to count number of reversals
nstep = 300;             % Number of time steps
% Take one backward Euler step to start Verlet
theta_old = theta - omega*tau;       
time = 0;
%%%%% MAIN LOOP %%%%%
for istep=1:nstep      
  t_plot(istep) = time;            % Record time and angle
  th_plot(istep) = theta*180/pi;   % for plotting
  accel = -g_over_L*sin(theta);    % Gravitational acceleration
  theta_new = 2*theta - theta_old + tau^2*accel;
  theta_old = theta;				% Verlet method
  theta = theta_new;  
  time = time + tau;
  if( theta*theta_old < 0 )  % Test position for sign change
    fprintf('Turning point at time t= %f \n',time);
    if( time_old < 0 )       % If this is first change,
      time_old = time;       % just record the time
    else
      irev = irev + 1;       % Increment the number of reversals
      period(irev) = 2*(time - time_old);
      time_old = time;
    end
  end
end
fprintf('Average period = %g +/- %g\n',mean(period),...
                          std(period)/sqrt(irev));
% Graph the oscillations
plot(t_plot,th_plot,'+');
xlabel('Time');
ylabel('Theta (degrees)');
title('Simple pendulum');