Week 4.1 - Solving second order ODEs

1] AIM -

To solve equation of motion of a simple pendulum with damping

2] OBJECTIVE-

1. To solve the second order ODE with graph of Velocity v/s Displacement.

2. To create the animation of the simple pendulum.

3] INTRODUCTION -

An ordinary differential equation is an equation containing one or more functions of one independent variable and the derivatives of those functions.The behavior of a system can be written in the form of a differential equation.One such example is the motion of a simple pendulum.
The ODE representing the motion of a simple pendulum is:

`(d^2 θ)/(dt^2 )+b/m (dθ)/(dt)+g/L sin⁡θ=0`

`theta`= angular displacement

b= damping coefficient

m= mass of the pendulum

g= acceleration due to gravity

L= length of the pendulum

4] PROCEDURE -

EQUATION BEING SOLVED -

First we have to convert the given second order ODE into two first order ODE'S as MATLAB can only sove first
order ODE's.

`theta_1 = theta` `theta_2 = (d theta)/dt`

`(d theta_1)/dt =( d theta)/dt` `(d theta_2)/dt = (d^2 theta)/dt`

`(d theta_1)/dt = theta_2` `(d theta_2)/dt = -b/m (d theta)/dt - g/l sin(theta)`

`(d theta_2)/dt = -b/m (theta_2) - g/l sin(theta_1)`

`theta_1` = Angular Displacement

`theta_2` = Angular Velocity

Therefore,we have our two first order ODE'S.

5] EXPLAINATION OF THE CODE

5.1) Main Code -

clear all
close all
clc
%INPUTS
l=1;       %LENGTH OF STRING
m=1;       %MASS OF BALL
b=0.05;    %DAMPING COEFFICIENT
g=9.81;    %GRAVITATIONAL ACCELERATION

t_span= linspace (0,20,658);%TIME STEP

theta_0=[0 ;3];%INITIAL CONDITION
%SOLVING ODE AND CALLING THE FUNCTION
[t,result]=ode45(@(t,theta)pendulum_function(t,theta,b,g,l,m),t_span,theta_0);
figure(1)
plot(t,result(:,1))% ANGULAR DISPLACEMENT
grid
hold on
plot(t,result(:,2)) %ANGULAR VELOCITY
xlabel('time')
ylabel('displacement,velocity')
legend('displacement','velocity')

% PENDULUM ANIMATION

x2=-0.7;
x3=0.7;
for i=1:length(result(:,1))
    x0=0;
    y0=0;
    x1=l*sin(result(i,1));
    y1=-(l*cos(result(i,1)));
    
    figure(2)
    plot([x0 x1],[y0 y1],'linewidth',2,'color','b') % CHORD
    hold on
    plot(x1,y1,'r.','markersize',70,'linewidth',2,'color','r') % BOB
    plot([x2 x3],[y0 y0],'linewidth',2,'color','b') % HORIZONTAL LINE
    hold off
    axis([-2 2 -2 2])
    pause(0.01)
   movie(i)=getframe(gcf);
  
end

%VIDEO FILE IN MPEG-4

V=VideoWriter('pendulum','MPEG-4');
open(V);
writeVideo(V, movie);
close(V);

5.2) Function -

%Function to solve pendulum ODE

%Function Definition
function [dtheta] = pendulum_function(t,theta,b,g,l,m)
dtheta=[theta(2);((-b/m)*theta(2))-((g/l)*sin(theta(1)))];
end

First all the given inputs were writen i.e

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b = damping coefficient.

Initial conditions(`theta_0`),

Angulardisplacement = 0 m

Angular velocity = 3 `(rad)/sec^2`

After that a function called as pendulun_function was created in separate script file in order to keep the length of final code small.In this function an array was created ,named as 'dtheta',which contains the both of our ODE's i.e 'theta1' = angular displacement and 'theta2' = angular velocity.And this function was saved by the same name i.e. pendulun_function

To solve the ODE'S MATLAB's inbuilt function called as 'ode45'is used.'ode45' can only solve first order ODE's using fourth order Runga-Kutta method.'ode45' requires function containing ode's ,time steps(t_span) to determine for how long and in how many steps we want to solve the ODE's,and initial condition(theta_0). After solving the ODE's output gets saved in variable called 'result',which contains values for angular displacement and angular velocity.

Then I plotted single graph plot comparing velocity and displacement with each other w.r.t. time .

From the graph we can see that both velocity and displacement gets damped as the time span increases.

5.3) Animation

To create the animation for the pendulum,first i used for loop and inside it created the 2D pendulum using the plot command and giving values for spacial co-ordinates. As we want to map change of displacement for all the time steps,I used simple trigonometry to get the value for x1 & x2,which represents the displacement of the ball.

Then first I plotted vertical line of the pendulum geometry with linewidth 2 and color blue,then using the hold
on command to stay in the same plot created ball using markersize command to give the ball radius and gave it
color red.And last created horizontal line of blue color.

Then I used hold off command so that, I get position of pendulum only at a peticular time step instead of getting
pendulum positions at all time steps at the same time and getting a messy animation.I then gave axis values ,which defines the teritory under which my bob should move,use pause of 0.01 sec to
smooth the animation between current and next time step.

5.4) Video file

I created a variable with index 'i' inside the for loop and used 'getframe' to create frame for each snapshot.Then used VideoWriter to write video data to a MPEG-4 file and store that in variable 'V'.Then recorded the video using writeVideo and closed the video.

Video link -https://www.youtube.com/watch?v=vtth_A0F66s

6] CONCLUSION

* The graph of Displacement and Angular velocity is plotted in Velocity v/s Time diagram.

* Animation of Simple pendulum is created.