Week 3 - Solving second order ODEs

Aim:

To solve the ODE which represents the equation of motion of a simple pendulum with damping.

Objecctive:

To write a Program that solves the following ODE which represents the equation of Motion of a simple Pendulum with Damping and create an animated video output obtained by solving this ODE.

Theory:

The Motion of a simple Pendulum is a Basic classical example of simple harmonic motion, consisting of a small bob and a mass less stirng. Movement of the pendulum depends on Newton's second law. When this law is written down, a second - order Ordinary Differential Equation is formed that describes the position of the end of pendulum with respect to time.

Fig : Simple Pendulum (with FBD)

Equation of Motion:

In the above Equation,

• g = gravity in m/s^2.
• L = length of the pendulum in meter.
• m = mass of the ball in kg.
• b = damping coefficient.

Since the Given function is in second order differential equation,we will convert it into first order to solve the ODE.

Let

`theta = theta1`

`d theta1/dt = theta2`

`d theta^2/dt^2 = d theta2/dt`

Input:

• L = 1 meter
  
• m = 1 kg
  
• b = 0.05
  
• g = 9.81m/s^2
  
• Time span = 0-20seconds
  
• Angular Displacement = 0
  
• Angular velocity = 3rad/sec^2 at time t = 0

Code:

clear all
close all
clc

% Inputs 
b = 0.05;
g = 9.81;
l = 1;
m = 1;

% Initial Condition
theta_0 = [0;4];

% Time Points
t_span = linspace(0,20,500);

% Solve ODE
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

% Plotting 
plot(t,results(:,1))
hold on 
plot(t,results(:,2))
xlabel('Time in second')
ylabel('Plot')

% Iteration using for loop to create animation of a simple pendulum 
ct = 1;
theta_1 = results(:,1);

for i = 1:length(results(:,1))
    
    % Initial Co ordinates
    x_start = 0;
    y_start = 0;
    theta_2 = theta_1(i);
    
    % Final co ordinates 
    x_end = 1*sin(theta_2);
    y_end = -1*cos(theta_2);
    
    % Plotting for Animation
    figure(2)
    % Base of Pendulum 
    plot([-1,1],[0,0],'linewidth',4,'color','black')
    hold on 
    grid 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',28,'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_order_ode.avi','uncompressed AVI');
open(videofile);
writevideo(videofile,M);
close(videofile);

Code for Function Program:

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 

Detail Explantion of the code:

1. In the first step,Input and initial condition of a simple penduluum are fed.
2. Then the program is written to solve the ODE in which matlab inbuilt function (ode45) is used with constant theta, t, and another function is used, which is defined seperately in function program.
3. In the next step result of solved ODE is plotted with respect to time using plot command.
4. Now Iteration of co-ordinate is fed using for loop to create animation.
5. Then plot command is used to plot the pendulum string and bob.
6. In the Next step, movie command is used to create an animated video of the oscillating pendulum.

Output:

* The output obtained is graphical representation velocity curve and displacement curve with respect to time.

* Angular Displacement and Angular Velocity plots were generated for a time range of 20 seconds.

Conclusion:

The Simulation of the transient behaviour of a simple pendulum is successfully craeted by using a program in Matlab.