Pendulum problem (main program source for PC and Mac workstation)

%
%   This program is for 4.170 in "An Introduction to Dynamics, 3rd ed."
%   By David J. McGill and Wilton W. King.
%

%
%   necessary input constant 
%

m = 1;                      %  mass of the bar
g = 9.81;                   %  gravity
l = 1;                      %  length of the bar

t_final = 5;                %  time duration for calculation
                            %  initial time is always 0 

%
%   Initial conditions for theta and theta_dot
%

theta = 1e-4;               %  initial displacement of the bar
                            %  initial disturbance in other words  
theta_dot = 0;              %  initial angular velocity of the bar

%
%   Setup integration solver (ODE45)
%

tol = 1e-8;                 %  tolerance for ODE45 solver
x0 = [theta,theta_dot];        

%
%  Perform ODE45 solver
%

[t,x] = ode45('p4_170_solver',0,t_final,x0,tol);
[n,m] = size(x);

Theta = x(:,1);
Theta_dot = x(:,2);

%
%  Calculation of reaction force at pin
%

Rx = -(m*g)*3/2*(1-cos(Theta));  %  reaction force at pin in x-dir
Ry = m*g*(1-3/4*sin(Theta)); %  reaction force at pin in y-dir

R = sqrt(Rx.^2+Ry.^2);       %  magnitude of reaction force at pin      

%
%  Plot the results
%

figure(1); clf;
plot(Theta*180/pi,Theta_dot);  
      
title('Problem 4.170');
xlabel('Displacement (degrees)');
ylabel('Angular velocity (radians/sec.)');
mx_theta_dot = max(Theta_dot);
axis([0 360 0 1.1*mx_theta_dot]);  grid;

figure(2); clf;
plot(Theta*180/pi,R);  
      
title('Problem 4.170');
xlabel('Displacement (degrees)');
ylabel('Reaction force at pin (Newton)');
mx_R = max(R);
axis([0 360 0 1.1*mx_R]);  grid;