Week 3 - Solving second order ODEs
AIM: Write a program in MATLAB\Octave that will simulate the pendulum motion.
DESCRIPTION: In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum.
Figure1: A simple pendulum.
A simple pendulum is one that can be considered to be a point mass suspended from a string or rod of negligible mass. It is a resonant system with a single resonant frequency.
GOVERNING EQUATIONS: The way the pendulum moves depends on Newton's 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.
`(dθ/dt)^2+b/m*dθ/dt+g/L*sin(θ)=0`
Parameters of the system will be defined by
g = gravity in m/s2
L = length of the pendulum in m,
m = mass of the ball in kg,
b=damping coefficient.
GIVEN:
L=1 metre,
m=1 kg,
b=0.05.
g=9.81 m/s
OBJECTIVES OF THE PROJECT
To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.
Figure2: Given the ODE of the pendulum(damped).
Figure3: Flow of the project.
CODE
Function to def
function [dtheta_dt] = ode_function(t, theta,l,m,b,g)
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
Main code
clear all
close all
clc
l = 1; %length of string
m = 1; %mass of bob
b = 0.05;%damping coefficient
g = 9.81;
%initialization
theta0 = [0;3];
%span of motion
t_motion = linspace(0,20,500);
%solving Ode using ode45 function
[t, results] = ode45(@(t,theta) ode_function(t,theta,l,m,b,g), t_motion, theta0);
figure(1)
plot(t,results(:,1),'color','g')
hold on
plot(t,results(:,2),'color','k')
plot(results(:,1),results(:,2))
ylabel('Angular Velocity')
xlabel('time')
legend('position Radian','angular velocity')
hold off
grid on
figure(2)
ct=1;
for i = 1:length(t_motion)
x0 = 0;
y0 = 0;
x1 = 1*sin(results(i,l));
y1 = -1*cos(results(i,1));
plot([x0 x1],[y0 y1],[-0.5 0.5],[0 0],'linewidth',3)
hold on
plot(x1, y1,'o','MarkerSize',30,'MarkerFacecolor','k')
axis([-1.5 1.5 -1.5 0.5])%defining axis
grid on
pause(0.0001)
hold off
M(ct) = getframe(gcf);%collecting frames
ct = ct+1;
end
movie(M)
videofile = VideoWriter('ODE_Pendulum.avi','Uncompressed AVI');% making movie
open(videofile)
writeVideo(videofile,M)
close(videofile)
EXPLANATION OF CODE
- Give input data as mentioned in the problem, as l,m,b,g.
- For initialization time and angular velocity as '0' & '3' used respectively.
- For defining the time span as 20 seconds and dividing it into 500 interval linspace command is used.
- By calling the function ode_function written in beginning the differential equation is solved.
- Plots are plotted with plot function and specified using the legend function in the plot.
- "Markersize" function is used to make the bob of the pendulum circular and "MarkerFacecolor" to the colour the same.
- getframe is used to collect frames and movie function to make a movie.
CONCLUSION
Figure4: variation of angular velocity with respect of time
- The amplitude of vibration keeps decreasing as it should in the damped vibration.
- The objective is achieved.
- This is the link for output video
https://drive.google.com/file/d/12JtuUApPkpuVpUK9aiy4olrnbxG-QlCF/view?usp=sharing