Week 3 - Solving second order ODEs
AIM - To solve 2nd Order differential equations using matlab.
OBJECTIVE - To write a program in Matlab\Octave that will simulate the pendulum motion.
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.
THEORY - An ordinary differential equation (also abbreviated as ODE), in Mathematics, is an equation which consists of one or more functions of one independent variable along with their derivatives. A differential equation is an equation that contains a function with one or more derivatives. But in the case ODE, the word ordinary is used for derivative of the functions for the single independent variable.
Simple Pendulum - Simple pendulam is another mechanical system that moves in an osciallatory motion. It consist of a point mass m suspended by mean of light inextensible string of length L from a fixed support.
EQUATION - d^2(theta)/(dt)+(b/m)*d(theta)/(dt)+(g/l)*sin(theta) = 0
b = damping coefficient, m = mass(kg), g = gravity(m/s^2), l = length of string(m)
BODY OF CONTENT - Below mentioned is the program that we used in solving 2nd order differential equation for simple harmonic motion of a pendulum.
% Solving 2nd Order ODE Using Matlab.
clear all
close all
clc
% Inputs
% Damping Coefficient
b = 0.05;
% Gravity
g = 9.81;
% Length of rope
l = 1;
% Mass of ball
m = 1;
% Initial Condition
% Angular displacement = 0;
% Angular velocity = 3;
theta_0 = [0;3];
% Time Points
% Generating 200 values between 0-20(sec) using linspace command
t_span = linspace(0,20,200);
% Solve Odes
% using ode45 which is inbuilt ode solver in matlab
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
% Plotting
% Angular displacement
plot(t,results(:,1))
hold on
% Angular velocities
plot(t,results(:,2))
xlabel('time')
ylabel('plot')
legend({'displace','velocity'})
% To defined pivot point of pendulam & plot it.
x0 = 0;
y0 = 0;
plot(x0,y0,'Marker','o')
% To extract the first coloumn of a result and assigned it to 'y'
y = results(:,1); % Theta/Angular displacement
% To create images for each time value and stored it in array 'm'.
ct = 1;
for i = 1:length(t_span)
x1 = l*sin(y(i));
y1 = -l*cos(y(i));
plot([x0 x1],[y0,y1],'LineWidth',3,'Marker','o')
hold on
plot(x1,y1,'marker','o','Color','g','markersize',30,'MarkerFaceColor','g')
hold off
title({'simple harmonic motion'; 'of a pendulam'},'color','r')
axis([-1 1 -2 0])
pause(0.03)
xlabel('time')
ylabel('angular displacement & angular velocity')
% All figures are converting into frames and are getting stored in array called 'M'
M(ct)= getframe(gcf);
ct=ct+1;
end
% Movie Creation
movie(M);
videofile = VideoWriter("second_order_odes.avi","Uncompressed AVI");
open(videofile);
writeVideo(videofile,M);
close(videofile);
EXPLANATION -
1) Program start with clear all, close all, clc commands.
2) Explained b,g,l,m with their units.
3) Inputs the initial condition for angular displacement & angular velocity.
4) Defined Time points t_span using 'linspace command' taken 200 values between 0-20 sec.
5) Used ODE45 inbuilt solver in matlab to calculate results of angular displacement & angular velocity for as above mentioned time period by giving initial conditions, also we have defined function which was saved in anotrher .m file and we called function inside ODE45 inbuilt solver.
6) Plots the result that came out from the solver in two coloums i.e first for angular displacement & second for angular velocity after that we had given labels legend & title to the plot.
NOW FOR ANIMATION PART -
1) we defined pivot points of pendulam & after that we plotted it using 'marker command'.
2) we extract the 1 coloum of the result (angular displacement) & assigned it to variable called y , as shown in code/program.
3) Now for each & every value of theta we had calculated x1 & y1 using for loop.
4) For loop used here for repetative action & it continues till the length of t_span as mentioned in code.
MOVIE PART -
1) Created a movie by using M(ct) = getframe(gcf)
where, M = array which store all the frame in sequence
ct = stands for countor - indexing all frames
getframe = This command is used to convert all the figures into frames.
2 ) After the end of for loop we called Movie(M) to create movie for all the frame, then we used variable called video file assigned to videowriter which is used to give empty file & profile name.
3) open(videofile) & close(videofile) is used to open & close the file in between writevideo command is used which fills the movie into empty vediofile -- writevideo(videofile,M);
RESULT - Below is the link for the animation of "Damped harmonic simple pendulum" that i had uploaded in my youtube channel.