Week 3 - Solving second order ODEs
Question- To solve an ordinary differential equation (ODE) of second order using Matlab and with the help of calculated information to
animate the motion of the simple pendulum.
Theory-
Pendulum have many applications and were utilized often before the digital age. They are used in clocks and metronomesw due to the regularity of their period,in wrecking balls and playground swings,sue to their simple way of building up and keeping energy. A simple pendulum consists of a point mass suspended on a string or wire that has negligible mass. If the pendulum weight or bob is pulled to a relatively small angle from the vertical and let go, it will swing back and forth at a regular period and frequency.
Below figure shows a simple pendulum model
The figure above shows a idealized pendulum,with a massless string/rod of length L and a bob of mass m. The open circle shows the rest position of the bob. When the bob is moved from its rest position and let go, it swings back and forth.The time it takes the pendulum to swing from its farthest right position to its farthest left position and back to its farthest right position is the period of the pendulum.
The primary forces acting on the bob are the gravitational force that make it move in the first place and the force exerted by the string to keep it be a damping force from friction at the pivot or air resistance or both. We will construct a model to describe how the angle theta of the pendulum varies as the function of the time t. let s(t) be the distance along the arc from the lowest point to the position of the bob at time t, with displacment to the right considered positive. Let theta(t) be the corresponding angle with respect to the vertical.The figure shows tangential and radial components of the gravitational force on the pendulum bob. The radial component is exactly balanced by the force exerted by the string,so the only relevant force producing the motion is the tangential component of the gravitational force. For the moment,we ignore the damping
force,if any. An ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variables
and the derivatives of those functions.The ODE equation and the calcualtions are shown below in figure .
Function Code Explanation-
Below is the function code
function [dtheta_dt]= ode_fun(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
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
A function is a group of statements that together perform a task. In matlab, functions are defined in seprate files. The name of the file and the functions should be the same .Functions operate on variables within their own workspace,which is also called as local workspace,seprate from the workspace you acess at the matlab command prompt which is called base workspace. Function handles are simple to function codes that can be created without using the regular function definition using the 'function' and 'end'. They do not require any seprate '.m' files or not needed to be declared seprately below/above the main code. Function handles are directly declared in the main code itself and also accessible from any point in that particular code,i.e theya are localized to a particular conditional. The declaration of a function handle is very similar to that of a function decleration. It can be distinguished by the means of '@'symbol before your value/expression allocation.the syntax for creating function handles is.a function was defined in order to establish a relation between all the given input values,the function code is based on the equation in above figure.
Main Code of the Project
%This program is to solve the differntial equation
% Solving differential equation of a simple pendulum for theta and angular
% velocity with time frame
clear all
close all
clc
% Inputs
b = 0.02;% damping coefficient
g = 9.81;% gravity(m/sec^2)
l = 1;% length of the pendulum(m)
m = 0.05;% mass of the bob(kg)
%intial condition
theta_0=[0;5]; % Intial displacement=0; intial velocity=5;
%time span
t_span=linspace(0,20,500);% time array in which motion is observed
%solve ode
[t,results]=ode45(@(t,theta)ode_fun(t,theta,b,g,l,m),t_span,theta_0);
%Plotting
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
title('Damped Harmonic motion of the pendulum')
legend('Angular displacement','Angular velocity')
ylabel('Plot')
xlabel('Time(second)')
%for loop
% Before starting the loop, a counter value was defined
ct=1;
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 the animation
figure (2)
plot([-1 1],[0,0],'linewidth',5)
hold on
plot ([ x_start x_end],[y_start y_end],'linewidth',3)
hold on
plot(x_end,y_end,'o','markersize',25,'markerFacecolor','m')
axis([-1.5 1.5 -1.5 1 ])
hold off
M(ct)= getframe(gcf);
ct=ct+1;
end
% Animation of the pendulum
movie(M)
Videofile=VideoWriter('Pendulum_animation.avi', 'UNCOMPRESSED AVI');
open(Videofile)
writeVideo(Videofile,M)
close(Videofile)
% Solving differential equation of a simple pendulum for theta and angular
% velocity with time frame
clear all
close all
clc
% Inputs
b = 0.02;% damping coefficient
g = 9.81;% gravity(m/sec^2)
l = 1;% length of the pendulum(m)
m = 0.05;% mass of the bob(kg)
%intial condition
theta_0=[0;5]; % Intial displacement=0; intial velocity=5;
%time span
t_span=linspace(0,20,500);% time array in which motion is observed
%solve ode
[t,results]=ode45(@(t,theta)ode_fun(t,theta,b,g,l,m),t_span,theta_0);
%Plotting
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
title('Damped Harmonic motion of the pendulum')
legend('Angular displacement','Angular velocity')
ylabel('Plot')
xlabel('Time(second)')
%for loop
% Before starting the loop, a counter value was defined
ct=1;
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 the animation
figure (2)
plot([-1 1],[0,0],'linewidth',5)
hold on
plot ([ x_start x_end],[y_start y_end],'linewidth',3)
hold on
plot(x_end,y_end,'o','markersize',25,'markerFacecolor','m')
axis([-1.5 1.5 -1.5 1 ])
hold off
M(ct)= getframe(gcf);
ct=ct+1;
end
% Animation of the pendulum
movie(M)
Videofile=VideoWriter('Pendulum_animation.avi', 'UNCOMPRESSED AVI');
open(Videofile)
writeVideo(Videofile,M)
close(Videofile)
Code Explanation-
First of all, it's important to write clear all,close all,clc commands since
clear all- removes data from the workspace ,freeing up the system memory
close all- deletes the current figure or specified figure
clc- clears all the input and output from the Command Window Display,giving you a clear screen.
In Next step we define the inputs and variables.i.e damping coefficient(b),gravity(g),length of the spring(l) and mass of the bob (m)
Next step for solving any differential equation is to define the intial conditions,the angular dispalcement and the angular velocities has
been defined. Timespan in which the motion is being visualized is also defiend.
Next step is to solve the ODE by using the ode45 solver, it is function is matlab that solve ordinary or partial differential equation using a numeric method.the ode45 id taking three inputs here , the function to be integrated,the time span,the intial condition.[ t,results] means that the result that have been calculated w.r.t to time.
After calculating the results, a graph has to be plotted, the graph contained two plots. One for angular displacment w.r.t time and one for
angular velocity w.r.t time. The legend and the Graph axis labels ar also defined.
Below figure shows the graph obtained
Below is the screenshot of the command window
The second part of the code is to establish a for loop, A for loop for the repetition control structures that allow you to efficiently write a loop
that needs to execute a specific number of times,a counter loop was added before for loop , the start and the end points for x and y axis has been defiend.
The loop is used where all the values of timespan are assigned to a variable'i' one at the time. The co-ordinates of [x0 x1 y0 y1] and their corresponding
relations are entered to form a 'string' of the pendulum.Marker size commands was used to form the bob of the pendulum.
A pause command was used for pausing the command for a certain time,for this project a pause of 0.03 was used. The M command was used to capture the current frame
every time the for loop run. Each command was used for the plot and for the loop to stop both from running.
The final step of the code is to create a movie of the animation of the simple pendulum,so the movie command was used.The Videowriter command was used to open the
video.writevideo command was used to write the video with plot frames.Finally the close commands was used to close the video file.
Conclusion - the ode is sucessfully solved and the graph has been plooted correctly the link of the animation of the pendulum is attaached below. The ODE for the pendulum is solved correctly and the graph and the animation video are also showing that vibration are damped in the nature and are diminshing with respect to the time. No variables and results are visible in the command window because semi colon is used almost evey step.
link - https://youtu.be/WQX5UBNUqNo