Week 3 - Solving second order ODEs

AIM : To write a program to solve the ODE which represents the motion of simple pendulum with damping

INTRODUCTION :

A second order diffrential equation is an ordinary differential equation in which any derivatives with respect to the independent variable have order no greater than 2. Ordinary diffrential equations helps to solve many enggineering systems. like SHM and damping.

A simple pendulum consists of a mass m hanging from a string of length l and fixed at a pivot point. When displaced to an initial angle theta and released the pendulum will swing back and forth with periodic motion

the equation is given by

`if θ =θ1

the dθ/dt = dθ1/dt = θ2-----(i)

d^2θ/dt^2 =d^2θ1/dt^2 = d/dt(dθ/dt) = d/dt(θ2)-------(ii)

d^2θ/dt^2 + b/m*dθ/dt + g/l*sinθ =0

let θ =θ1

d^2θ1/dt^2 + b/m*θ2 +g/l*sinθ1

dθ1/dt = θ2

dθ2/dt = -b/m*θ2-g/l*sinθ.

CODE EXPLAINATION :

Inputs given in the problem are

L = 1 metre

g = 1kg

b = 0.05

g = 9.81m/s

Intial condition taken as

angular displacement = 0, angular velocity = 3 rad/sec, at time t = 0 the motion is simulated between 0-20 sec

solving ODE : ode solver ode45 is used in MATLAB program. The results are stored in an array of row of two columns for the vlaues of the displacement and velocity at diffrent time periods.

Figure 1: plotting

here the graph is plotted to get the wave form of the displacement and the velocity from the initial conditions

figure 2 : animation

here the concept of for loop as the bob moves with a velocity and displaces for the diffrent time period and recoreded it and animated the path of oscillations

%solving second order ODEs of a pendulum

clear all 
close all
clc

b = 0.05;
g = 9.81;
l = 1;
m = 0.1;

%intial condition 
%angular displacement is 0 and angular velocity is 5

theta_0 = [0;5];

%time span
t_span = linspace(0,20,600);

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

%plotting
figure(1)
hold on
plot(t,results(:,1))
plot(t,results(:,2))
title('Damped harmonic motion of a pendulum')
legend('Angular displacement', 'Angular Velocity')
ylabel('plot')
xlabel('Time(seconds)')
ct =1;
for i =1:length(results(:,1))
    x_start = 0;
    y_start = 0;
    x_end = l*sin(results(:,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 pendulum

movie(M)
videofile=VideoWriter('pendulum animation.avi' , 'uncompressed avi');
open(vidoefile)
writeVideo(videofile,M)
close(videofile)