Abhishek S. Jadhav
A) Aim:- Solving second order ODEs of simple pendulum with damping & write a program in Matlab that will simulate the pendulum motion.
B) Objective:-
- Solving second order ODEs of simple pendulum with damping using MATLAB
-
Write a program in Matlab that will simulate the pendulum motion as shown below .
C) Theory:-
- In Engineering, ODE is used to describe the transient behaviour of a system. A simple example is a pendulum
- The way the pendulum moves depends on the Newtons second law.
- When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time. Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of theses systems.
1) Governing Equation:-
`(d^2 θ)/(dt^2 )+b/m×dθ/dt+g/L sinθ=0 …..(1)`
Where,
b=damping coefficient.
m=mass of the ball `(kg)`
g=gravity `(m/s^2)`
L=length of the pendulum `(m)`
2) Converting 2nd order ODE equation to 1st order ODE equation
`θ=θ_1 …..(2)`
`dθ/dt=(dθ_1)/dt=θ_2 …..(3)`
`(d^2 θ)/(dt^2 )=(d^2 θ_1)/(dt^2 )=d/dt ((dθ_1)/dt)=d/dt (θ_2 ) …..(4)`
Now, put `(d^2 θ)/(dt^2 )=(d^2 θ_1)/(dt^2 ), dθ/dt=θ_2 & θ=θ_1`in eq. (1)
`(d^2 θ)/(dt^2 )+b/m×(dθ_1)/dt+g/L sinθ=0 …..(5)`
Again put, `(d^2 θ)/(dt^2 )=(dθ_2)/dt & (dθ_1)/dt=θ_2` in eq. (5)
`(dθ_2)/dt+b/m θ_2+g/L sinθ=0 …..(6)`
Above equation break into two different ODE
`(dθ_2)/dt=-b/m θ_2-g/L sinθ` And `(dθ_1)/dt=θ_2`
Converting above ODE eq. in matrix form
C) Given data for programing:-
L=1 m,
m=1 kg,
b=0.05
g=9.81 `m/s^2`
D) Coding:- Coding will be done in two part as following
1) Main programing
2) Function programing
1) Main programing:-
- First we take input the initial condition of pendulum I.e. initial angular displacement=0,initial angular velocity=3 rad/sec at initial time t=0 & total time is t=20sec
- We use MATLAB inbuilt function ‘ode45’ to solve ODE
- ODE solver ‘ode45’ uses input data t, theta, b, g, l, m which is defined in function program
- Then we solve 2nd order ODE , plot graph time (sec) vs velocity (rad/sec)
- For creating animation, we uses iteration of coordinate is feed using ‘for’ loop
- To create animation string & bob, we uses ‘plot’ & ‘hold on’ command
- Creating movie by using movie, VideoWriter, open, WriteVideo & close command
- Program
clear all
close all
clc
b=0.05;
g=9.81;
l=1;
m=1;
%initial condition
theta_0=[0;4];
%time point
t_span= linspace(0,20,500);
%solve ODE
[t,result] = ode45(@(t,theta) ode_func(t, theta, b,g,l,m), t_span, theta_0);
plot(t,result(:,1))
hold on
plot(t,result(:,2))
ylabel('Velocity (rad/sec)')
xlabel('Time (sec)')
%Create animation of simple pendulum
ct=1;
theta_1 = result(:,1);
for i= 1:length(result(:,1))
%initial coordinate
x_start= 0;
y_start= 0;
theta_2 = theta_1(i);
%final coordinate
x_end = 1*sin(theta_2);
y_end = -1*cos(theta_2);
%ploting for animation
figure(2)
%base of pendulum
plot([-1,1],[0,0], 'linewidth',4,'color','black')
hold on
%string of pendulum
line([x_start x_end],[y_start y_end],'linewidth',6,'color','b')
hold on
%bob of pendulum
plot(x_end,y_end,'marker','o','markersize', 25,'markerfacecolor','r')
axis([-2 2 -2 2])
pause(0.03)
hold off
M(ct)=getframe(gcf);
ct=ct+1;
end
%creating movie from captured frames
movie(M)
videofile = VideoWriter('Second_ODE.avi','uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);
- Screenshot of main program & output graph
- Output graph
- Output animation
2) Function programing:-
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
- Screenshot of main program
D) Error encountered:-
- Solution for error:- As show above screenshot error occured due to wrong spelling of command of linespace instead of linspace.
E) Conclusion:- By using we MATLAB we show transient behavior of simple pendulum.