7 views (last 30 days)

I have a code written out to simulate a pendulum, but there seems to be an error after the loop, and I can't seem to resolve it on my own. I would appreciate any advice on how to solve this issue. My professor can't seem to solve this issue either, and as a result, has told me that I must solve it without further assistance. I have everything written out with comment lines to help the user. I have also included the error message I get. Thanks!

Unrecognized function or variable 'period'.

Error in pendul (line 59)

AvePeriod = mean(period);

%pendul - Program to compute the motion of a simple pendulum

%using either the Euler or Verlet methods

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 the number of reversals

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, just record the time

time_old = 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 oscillations, 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

clf; figure(gcf); %Clear and forward figure window

plot(t_plot,th_plot, '+');

xlabel('Time'); ylabel('\theta (degrees)');

Walter Roberson
on 27 Sep 2020

for istep=1:nstep

if nstep is less than 1 then the loop will not be executed at all and period will not be assigned to.

If nstep is at least 1 then we need to investigate whether we can be certain that period is assigned to.

if( theta*theta_old < 0)

If we have not executed enough steps for the angle to have crossed zero then period would not have been assigned to.

if(irev == 0)

The first time we cross 0, we do not assign to period.

Putting these together: period will not be assigned to until you have executed enough steps to have crossed zero twice.

if( theta*theta_old < 0)

No crossing will be detected if theta or theta_old are exactly 0 because then their multiple would be exactly 0 rather than less than 0.

You need to be careful about this case so that you do not accidentally detect the same crossing twice.

Walter Roberson
on 27 Sep 2020

Increase the number of steps you allow. Each step only goes a fraction of a cycle.

But also before you output check that irev is at least 2, because if not then you have not found the period yet.

Opportunities for recent engineering grads.

Apply TodayFind the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 0 Comments

Sign in to comment.