MATLAB Answers

How to stop a simple harmonic oscillator after one period

2 views (last 30 days)
Ash matlab
Ash matlab on 23 Jan 2020
Commented: Walter Roberson on 24 Jan 2020
Hi Matlabers
I have this code and I was wondering how should I stop the code once the pendulum makes one period. I thought an if statement saying when t_plot = 2*pi break out of the code but nothing worked....
I appreciate all help
Sorry for the long code....
theta0 = input('Enter initial angle (in degrees): ');
theta = theta0*pi/180; % Convert angle to radians
omega = 0; % Set the initial velocity
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)');

  3 Comments

Walter Roberson
Walter Roberson on 24 Jan 2020
You have finished a half period when the pendulum changes direction.
Walter Roberson
Walter Roberson on 24 Jan 2020
If you start from negative theta then the pendulum sweeps towards more positive values of theta until it reaches the peak and then starts reversing direction. When it starts reversing direction toward negative then theta_old > theta becomes true.
But watch out for the case where it starts out positive...

Sign in to comment.

Answers (1)

KSSV
KSSV on 24 Jan 2020
Knowing mass and spring constant of the pendulum, you can calculate the frequency, time period. Once time period is known, you can fix the time and stop wherever you want.
clc; clear all ;
m = 10 ;
k = 5. ;
% frequency
f = 1/(2*pi)*sqrt(k/m) ;
% Time period
T = 2*pi*sqrt(m/k) ;
% Amplitude
A = 1. ;
t = linspace(0,T) ;
x = A*sin(2*pi*f*t) ;
plot(t,x) ;

  0 Comments

Sign in to comment.

Sign in to answer this question.

Tags