Solution procedure of pendulum problem using Matlab

This is an example of Matlab program. The simple programming techniques in Matlab to solve the problem and some results will be shown.
Problem Statement

The drawing below shows the bar which is pin-jointed to the vertical wall. When the bar stands at vertical position, the bar is in equilibrium position. If there exists even a very small disturbance at that position, the bar will starts rotating. What we will do is to describe the motion of the bar and the forces which act on the bar when it rotates.

Free Body Diagram and Equation of Motion

The approach is same as we did in a piston problem . Kinematics first, then setting up the equations of motion from kinematics, force and moment balance. In this problem, we have kinematic restraint like the bar is rotating about a pin joint, although kinematics does not play an important roll as it did in a piston problem.

You can see the procedure for writing a equation of motion for this problem below

In this problem, we have a pendulum pin-jointed to a vertical wall falling from equilibrium due to a small disturbance. There is not much to the kinematics of the problem. The only kinematic constraint we have is that the bar cannot be separated from the pin joint, but it does not affect the equation of motion.

The force balance equation from the FBD above can be written as folows.
\begin{displaymath}I\alpha=I\ddot{\theta}=mg\frac{l}{2}\mbox{sin}\theta\end{displaymath} (1)
where the moment of inertia of the bar about the pin joint is
\begin{displaymath}I=I_G+m\left(\frac{l}{2}\right)^2=\frac{ml^2}{3}\end{displaymath} (2)
Therefore the equation of motion is
\begin{displaymath}\ddot{\theta}-\frac{3g}{2l}\mbox{sin}\theta=0\end{displaymath} (3)

Closed Form Solution

Notice that the equation of motion is non-linear due to the sin(theta) term. We could linearize this equation for small motions, since theta is approximately equal to sin(theta) if theta is between 0 and 10 degrees. For general motions, the closed form solution in terms of elliptic integrals. This is not a "nice" closed form solution as in the linear vibrations problems, so it will be best for our purposed to use a numerical solver. Again we go to Matlab.

Integrating EOM with Matlab

How do we solve Ordinary Differential Equations (ODE) using Matlab? If you don't don't know, please read Matlab Assistance for using ODE45. There, you will find most of basic things to use ODE45 solver in Matlab.

We will follow the same procedure that has been used in part 3, piston problem. The following is a Matlab program source codes. The first few lines are for constant input values like mass of the bar, etc. The numbers given are all arbitrary. 

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

tspan = [0,5];              %  time duration for calculation
Now, plug in the initial conditions for main variable, theta in this case, the rotation angle of bar. 
theta = -1e-4;              %  initial displacement of the bar
                            %  initial disturbance in other words  
theta_dot = 0;              %  initial angular velocity of the bar
I used small disturbance (as we see the above, 1e-4 rad.) to initiate the motion of the bar, which was in equilibrium. The reason for the negative initial value for theta is that theta is set up as positive in counter-clockwise direction and the bar rotates clockwise direction. You can use a velocity input also, instead of disturbance in displacement. 
x0 = [theta,theta_dot]; 

options = odeset('Refine',6,'RelTol',1e-4,'AbsTol',1e-7); 
  
%  options 'Refine' was used to produce good-looking plots
%  default value of Refine was 4.
%  Refine doesn't apply if length(TSPAN) > 2.
%  scalar relative error tolerance 'RelTol' (1e-3 by default) and 
%  vector of absolute error tolerances 'AbsTol' (all components 1e-6 by
%  default).
'Refine' option does not apply if time span is more than 2. After we specify all the inputs for ODE45 solver, we write 
[t,x] = ode45('pend_sol',tspan,x0,options);
Plotting Results

After finishing the calculation, you plot them. Here are some of the plotting commands. 

figure(1); clf;
plot(abs(Theta)*180/pi,abs(Theta_dot),'b');  
title('Pendulum problem');
xlabel('Displacement (degrees)');
ylabel('Angular velocity (radians/sec.)');
mx_theta_dot = max(abs(Theta_dot));
axis([0 360 0 1.1*mx_theta_dot]);  grid;
Absolute values for Theta and Theta_dot are taken to have nice looking plot. Now, you've got them all.
Please see Matlab tutorial piston problem , if you have no idea what the commands above mean.

Here, the complete source codes of main program and subprogram for Matlab v5.0.

In the plot, the angular velocity of the bar is shown as the rotation angle varys.

If you need some more information about commands, try "help". For example, 

help plot
 

Return to Example Table