Week 3 - Solving second order ODEs

Second Order ODEs Challenge (Simple Pendulum)

Introduction:

A pendulum is a weight suspended from a pivot so that it can swing freely. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth. The time for one complete cycle, a left swing and a right swing, is called the period. The period depends on the length of the pendulum and also to a slight degree on the amplitude, the width of the pendulum's swing.

Objective:

1. Simulating Second Order ODEs of Simple Pendulum

2. Plotting Results of Angular Velocity and Displacement

3. Creating Animation of the Simple Pendulum using Matlab Plot.

Solution:

Notations:

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/s2.

Equation of Simple Pendulum

Programe Code

1. ODE Function

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

Explaination:

The above function simply transform the second order differential equation to first order differential equation by introducing the variable `theta2` and the equation is transformed to first order ODE in terms of 'theta1 and theta2 such that they are differentiated to yield a colum matrix.

2. Main Program

clear all
close all
clc

b=0.05; % Damping Cofficient
g=9.81; % Acceleration Due to Gravity
l=1;  % Length of the Pendulum
m=0.1; % Mass of the Bob of the Pendulum

% initial condition
theta_0=[0;5];

% time points
t_span =linspace(0,10,500);

% solve ODE

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

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

Explaination:

In the above code the results are been plot by solving the Transformed differential Equation by using ODE45 and ODE_func function as shown. This results indicated the Velocity and the Displacement of the Pendulum with respect to Time.

3. Animation

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',20,'markerFacecolor','y')
     title('Animation of transient behaviour of a Pendulum')
     axis([-1.5 1.5 -1.5 1])
     pause(0.03)
     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)

Explaination

The above code represents the process of Animation of the Simple Pendulum in which multiple plots of simple Pendulum with different orientation are generated keeping one end of the Pendulum fix as shown in the Plot. This animation is further saved as an AVI file named Pendulum Animation.avi

Result Plot

The above plot shows the variation in the Velocity and the Displacement of the Pendulum.

Youtube Animation Link:

https://www.youtube.com/watch?v=00-WBqQwgDU