MATLAB PROGRAM FOR 2nd ODE of Pendulum motion.

1. AIM: -

To simulate the 2^nd ODE of Pendulum motion.

2. OBJECTIVE

• Write a program in Matlab that will simulate the pendulum motion.
• Create an animation for plot obtain in MATLAB program.

3. GOVERNING EQUATION

2^nd ODE of simple pendulum with damping.

`(d^2 θ)/(dt^2 ) + b/m*(dθ)/dt + g/l sinθ`,

were,

b = damping coefficient

m = mass of bob

g = Acceleration due gravity (m/s^2)

l = length of string in meter

4. DESCRIPTION

A simple pendulum can be described as a device where its point mass is attached to a light inextensible string and suspended from a fixed support. The vertical line passing through the fixed support is the mean position of a simple pendulum. The vertical distance between the point of suspension and the centre of mass of the suspended body (when it is in mean position) is called the length of the simple pendulum denoted by L. This form of the pendulum is based on the resonant system having a single resonant frequency.

4.1 Term related to simple pendulum

• Length (L): Distance between the point of suspension to the center of the bob

• Time period (T): Time taken by the pendulum to finish one full oscillation

• Displacement (x): Distance travelled by the pendulum bob from the equilibrium position to one side. The angle described by the pendulum with an imaginary axis at the equilibrium position is called the angular displacement (θ).

• Amplitude (x _max): Maximum distance travelled by the pendulum from the equilibrium position to one side before changing its direction. For angle, it is denoted by θ _max.

4.2 2^nd order differential equation

`(d^2 θ)/(dt^2 )+b/m*(dθ /dt)+g/l sinθ`

We had to solve this because we know that Matlab can only solve 1^st order differential equation.

Assume

`θ=θ1`

1. `(dθ)/dt=(dθ1)/dt=θ2`
2. `(d^2 θ)/(dt^2 )=(d^2 θ1)/(dt^2 )=d/dt [(dθ1)/dt]=(dθ2)/dt`

From equation (1) & (2)

`(dθ2)/dt+b/m θ2+g/l sinθ1=0`

`(dθ2)/dt=-b/m θ2-g/l sinθ1`

In matrix form

`d/dt [(θ1.θ2)]=[θ2;(-b/m-g/l sinθ1)]`

From above equation we can get i.e. angular displacement.

5. PROGRAM FOR SOLVING DIFFERENTIAL EQUATION

We are using function command here for solve differential equation as follow,

function[dtheta_dt] = ode_func(t,theta,b,g,l,m)

theta1= theta(1);
theta2 = theta(2);
dtheta1_dt = theta2;
dtheta2_dt = -(b/m)*theta2-(g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt;dtheta2_dt];
end

Now, we are going to use this Program in our first program for animation correctly, but using total program to our first program is not correct way therefore, we are creating a function command here. To use this command in maim program, we have to use same name for function. here we only put our equation in function command to solve equation so we can get angular displacement.

6. EXPLANATION OF PROGRAM

• Define Input

All values are given in problem statement as follow,

b = 0.05; % Damping coefficient

g = 9.81; % Gravity

l = 1; % Length of rod

m = 1; % Mass of ball

• Initial Condition

theta_0 = [0;3];

Here we are use initial angular displacement is 0 and initial angular velocity 3 (m/s^2).

• Time Interval

t_span = linspace(0,20,300);

Linspace(X1,X2,N): - N is equally space value between X1 and X2. here we were used time array for pendulum oscillation 20 second with 300 intervals.

Here we are taking multiple time interval. Initially, we are providing the inputs for the time with the help of linspace command as above.

• ODE45 Solver

Most of the time. ode45 should be the first solver you try. ode45 is based on an explicit Runge-Kutta (4,5) formula. in ode solver we were used function command for appropriate results.

[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

here ode45 is internal solver of Matlab. There are other constant also but (t) and (theta) is dependant and independent variable of our ode system. So, we are said ode is function of (t, theta). Here (t) is time and results would be and displacement and velocity.

• Plotting Figure 1

% Plotting

figure(1)

plot(t, results(:,1));

hold on

plot(t,results(:,2));

legend('Displacement', 'Velocity')

title('Damping of Simple Pendulum')

xlabel('Time')

ylabel('Plot of Displacement and Velocity')

we used following command,

plot command to plot figure graph.

hold on command to creat graph in same graph.

grid on command to display grid in graph.

title command for give title to graph.

legend command is used for give legend to graph.

Xlabel and ylabel command for label x and y axis.

6.1 OUTPUT

• Explanation

Here we can see the pendulum oscillation are damping gradually. We are damp only for 20 second if we increase time, we can saw complete damping of pendulum.

7. PROGRAM FOR PLOT PENDULUM ANIMATION

For animation we were created following command,

To change the result i.e. displacement and velocity, we are using For-loop command as shown below. Now, here is the value of (i) changes we are taking the values of results from array (i).

We use the value of (i) to take particular value for results and then we are going to perform animation.

% For loop

for i = 1:length(results(:,1))

x_start = 0;

y_start = 0;

x_end = l*sin(results(i,1));

y_end = -l*cos(results(i,1));

• plotting for animation

here we are plot support line, string line and bob.

figure(2)

% Base line

plot([-1.5 1.5],[0 0],'linewidth',10);

hold on

% Pendulum line

plot([x_start x_end],[y_start y_end],'LineWidth',5);

hold on

% Pendulum ball

plot(x_end,y_end,'o','markersize',15,'markerfacecolor','r');

axis ([-2 2 -2 1]);

pause(0.03)

grid on

7.1 Output

Here we get our 1^st error. We put hold on command for got all plot on same graph so we have to use hold off command at the end. Also, we used pause command to get require output for limited time.

7.2 Error correction

% Pendulum ball

plot(x_end,y_end,'o','markersize',15,'markerfacecolor','r');

axis ([-2 2 -2 1]);

grid on

pause(0.03)

hold off

7.3 Output

We get animation of pendulum as below,

8. MOVIE FOR ANIMATION: -

To animate we are adding figure(ct) before plotting and outside for-loop ct=1 And after the plot command type ct=ct+1 (i.e. we incrementing the value each time program gets executed in the for-loop).

To create movie type M(ct)=getframe(gcf) just below the pause(0.03) and at below end command type Movie(M).

Where,

ct : - every frame

getfram : - getting the current figure in the gcf window.

Further, we are going to save all the figures in the array called (M).

To create movie, we are using following command,

movie(M);

videofile = VideoWriter('ODE_Pendulum.avi');

open(videofile);

writeVideo(videofile,M);

close(videofile);

we can use ‘uncompressed AVI’ command but it gives us the large avi file.

9. FINAL PROGRAM

clearvars
close all
clc

b = 0.05; % Damping coefficient
g = 9.81; % Gravity
l = 1; % Length of rod
m = 1; % Mass of ball

%initial condition
theta_0 = [0;3];

%time point
t_span = linspace(0,20,300);

%solve ODE
[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

% Plotting
figure(1)
plot(t, results(:,1));
hold on
plot(t,results(:,2));
legend('Displacement', 'Velocity')
title('Damping of Simple Pendulum')
xlabel('Time')
ylabel('Plot of Displacement and Velocity')
ct = 1;

% For loop
for i = 1:length(results(:,1))
    x_start = 0;
    y_start = 0;
    x_end = l*sin(results(i,1));
    y_end = -l*cos(results(i,1));
    
    % plotting for animation
    figure(2)
    % Base line
    plot([-1.5 1.5],[0 0],'linewidth',10);
    hold on

    % Pendulum line
    plot([x_start x_end],[y_start y_end],'LineWidth',5);
    hold on

    % Pendulum ball
    plot(x_end,y_end,'o','markersize',15,'markerfacecolor','r');
    axis ([-2 2 -2 1]);
    grid on
    hold off
    pause(0.03)

    M(ct) = getframe(gcf);
    ct = ct+1;

end

% Creating video
movie(M);
videofile = VideoWriter('ODE_Pendulum.avi');
open(videofile);
writeVideo(videofile,M);
close(videofile);

10. OUTPUT OF PROGRAM

Here is the link for movie of simple pendulum

https://drive.google.com/file/d/1pkdTsPpNHsK3Wrfaj9vWgqDB3A7hnOl3/view?usp=sharing

11. CONCLUSION

Hence, we can conclude that we are able to Program 2^nd ODE of simple pendulum. with the help of this program. By small adjustments in the program, we can use it in multiple applications. Such as damping coefficient, time, angular velocity, angular displacement.