Week 3 - Solving second order ODEs

AIM: To solve the ordinary differential equation (ODE )which represents the equation of motion of a simple pendulum with damping.

OBJECTIVE: 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 of output obtains by solving the ODE.

THEORY: The motion of a simple pendulum is a basic classicsal example of simple harmonic motion, consisting of small bob and a massless string. Movement of pendulum depends upon Newtons second law of motion. when this law is written down , a second order differential equation is formed that describes the position of the end of the pendulum with respect to time .

In the above equation,

g = gravity in m/s^2

L = length of pendulum in meter

m = mass of ball in kg

b = damping coefficient

since the given function is in second order differential equation we will convert into first order differential equation to solve the ODE.

Let `theta`= `theta1`
d`theta1`/dt = d`theta2`/dt

INPUT :

L = 1 metre

m =1 kg

b=0.05

g=9.81m/s^2

Time span = 0-20 sec

angular displacement =0

angular velocity =3rad/sec^2 at time t = 0

CODE:

clear all
close all
clc

%inputs 
b = 0.05;   %damping coeffiecient (unitless)
g = 9.8;    %gravity (m/s^2)
l = 1;      %lenth of pendulum string (meter)
m = 1;      %mass (kg)

%initial condition
theta_0 = [0;3];

% 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
hold on
plot(t,results(:,1))
plot(t,results(:,2))
title('damped harmonic motion of a pendulum')
legend('Angular displacement, Angular velocity')
xlabel(' Time in seconds')
ylabel('plot')

%Iteration using for loop to create animation of simple pendulum 
ct = 1;

for i= 1:length(results(:,1))
    
    %initial co-ordinates
    x_start = 0;
    y_start = 0;
  
    %final co-ordinates
    x_end = 1*sin(results(i,1));
    y_end = -1*cos(results(i,1));
    
    %plotting for animation
    figure(2)
    %base of pendulum
    plot([-1,1],[0,0],'linewidth',4,'color','k');
    %string of pendulum
    hold on
    plot([x_start x_end],[y_start y_end],'linewidth',6,'color','b')
    %bob of pendulum
    hold on
    plot(x_end,y_end,'o','markersize',25,'markerFacecolor','g')
    title('Animation of transient behaviour of a pendulum')
    axis([-1 1 -1.5 1.5])
    pause(0.05)
    hold off
    hold off
    grid on
    M(ct) = getframe(gcf);
    ct = ct+1
end
%ANIMATION
movie(M)
videofile = VideoWriter('second_order_pendulum')
open(videofile)
writeVideo(videofile,M)
close(videofile)
    

CODE FOR 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

DETAIL EXPLANATION OF CODE:

1) In the first step input and initial condition of a simple pendulum are given.

2)Then the program is written to solve the ODE in which inbuilt function (ODE45) is used with constant theta ,t , and another function is used whuch is defined separety in function program.

3)In the next line result of ODE is plotted with respect to time using plot command.

4) Now iteration of co-ordinate is given using for loop to create animation.

5) The plot command is used to plot pendulum string and bob.

6) In the next step, the movie command is used to create an animated movie and the oscilatting pendulum.

OUTPUT :

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

(.)Angular dispalcement and angular velocity plots were generated for a time range of 20 seconds.

Various errors faced during programming are:

1) missed applying comma after 'color'

Step taken to overcome this error: by applying comma at the right place.

2) spelling mistake in writing :

steps taken to overcome this error:

by correcting the spelling of length.

CONCLUSION: The simulation of transient behaviour of a simple pendulum is successfully created by using a MATLAB program.