%% pendul - Program to compute the motion of a simple pendulum % using the Euler or Verlet method clear all; help pendul % Clear the memory and print header %% * Select the numerical method to use: Euler or Verlet NumericalMethod = menu('Choose a numerical method:', ... 'Euler','Verlet'); %% * Set initial position and velocity of pendulum theta0 = input('Enter initial angle (in degrees): '); theta = theta0*pi/180; % Convert angle to radians omega = 0; % Set the initial velocity %% * Set the physical constants and other variables g_over_L = 1; % The constant g/L time = 0; % Initial time irev = 0; % Used to count number of reversals period = []; % Used to record period estimates tau = input('Enter time step: '); %% * Take one backward step to start Verlet accel = -g_over_L*sin(theta); % Gravitational acceleration theta_old = theta - omega*tau + 0.5*tau^2*accel; %% * Loop over desired number of steps with given time step % and numerical method nstep = input('Enter number of time steps: '); for istep=1:nstep %* Record angle and time for plotting t_plot(istep) = time; th_plot(istep) = theta*180/pi; % Convert angle to degrees time = time + tau; %* Compute new position and velocity using % Euler or Verlet method accel = -g_over_L*sin(theta); % Gravitational acceleration if( NumericalMethod == 1 ) theta_old = theta; % Save previous angle theta = theta + tau*omega; % Euler method omega = omega + tau*accel; else theta_new = 2*theta - theta_old + tau^2*accel; theta_old = theta; % Verlet method theta = theta_new; end %* Test if the pendulum has passed through theta = 0; % if yes, use time to estimate period if( theta*theta_old < 0 ) % Test position for sign change fprintf('Turning point at time t= %f \n',time); if( irev == 0 ) % If this is the first change, time_old = time; % just record the time else period(irev) = 2*(time - time_old); time_old = time; end irev = irev + 1; % Increment the number of reversals end end %% * Estimate period of oscillation, including error bar AvePeriod = mean(period); ErrorBar = std(period)/sqrt(irev); fprintf('Average period = %g +/- %g\n', AvePeriod,ErrorBar); %% * Graph the oscillations as theta versus time figure(1); clf; % Clear figure window #1 and bring it forward plot(t_plot,th_plot,'+'); xlabel('Time'); ylabel('\theta (degrees)');