Week 3 - Solving second order ODEs

%This program is to solve the ODE representing the motion of a simple pendulum with damping and hence simulate the motion for 0 to 20 secs.
close all 
clear all 
clc

%%
%Problem Statement

%Governing Differential Equation is given by θ"+(B/M)θ'+(g/L)sinθ=0.
%where,
%θ is angular displacement
%θ' and θ" is 1st and 2nd derivative of angular displacement with respect to time respectively. 
%B is the damping coefficient.
%M is the mass of the ball in kg.
%g is the acceleration due to gravity in m/s^2.
%L is the length of the string attached between the wall and the ball in m.

%%
%Solving for ODE

%Inputs
M=1;                        %mass of the pendulum sphere in kgs.
B=0.05;                     %damping coefficient.
L=1;                        %length of the pendulum in m.
g=9.81;                     %acceleration due to gravity in m/s^2.
ini_cond=[0;3];             %initial condition of θ=0 and θ'=3 rad/s as column vector.
t_span=linspace(0,20,500);  %time span for which we need to solve for the ODE is 0 to 20 sec.

[t,results]=ode45(@(t,theta) pend_ode(M,B,L,g,theta,t),t_span,ini_cond);

%%
%Plot generation

hold on
plot(t,results(:,1));
plot(t,results(:,2));
xlabel('time');
ylabel('magnitude');
legend('displacement','velocity');
title('Plot for "Dispacement vs time" & "Velocity vs time"');

%%
%Animation generation

x0=0;                                                     %x coordinate for the fixed end of the pendulum.
y0=0;                                                     %y coordinate for the fixed end of the pendulum.
figure(2)                                                 %figure 2 for the animation of the pendulum motion.
cntr=1;
for i=1:length(t_span)                                    %for loop for generating the plot for displacement. of the pendulum for 20 secs interval
    x1=L*sin(results(i,1));                               %x coordinate for the moving end of the pendulum.
    y1=-L*cos(results(i,1));                              %y coordinate for the moving end of the pendulum.
    plot([x0,x1],[y0,y1],x1,y1,'r.','MarkerSize',30,'LineWidth',2);     %plotting the moving end of the pendulum.  
    axis([-2 2 -2 0]);                                    %axis range for the figure is set for x from -1  to 1 and y from -2 to 0.
    s1='Time=';                                          
    s2=num2str(t_span(i));                                %converting value in t(i) to string.
    text(x1+0.05,y1-0.05,[s1,s2 'seconds']);              %labeling time for a given position of the pendulum.
    MV(cntr)=getframe(gcf);                               %each plot is stored in an array MV.
    cntr=cntr+1;                                          %index for the array MV is increased by 1.
    z(cntr)=x1^2+y1^2;                                    %verifying length of string is constant.
end
movie(MV);                                                %movie is generated from array of figures.
Animation=VideoWriter('Pendulum_motion.avi');             %video file is created in .avi format.
open(Animation);                                          
writeVideo(Animation,MV);                                 %the movie MV is written into the .avi file.
close(Animation);


        

%This function returns the array 'detheta' as output,
%while taking mass (M), damping coefficient (B), length of the string (L), angular displacement (theta) and time (t) as input. 

function [dtheta]=pend_ode(M,B,L,g,theta,t)
dtheta1_dt=theta(2);
dtheta2_dt=(-g/L)*sin(theta(1))-(B/M)*theta(2);
dtheta=[dtheta1_dt;dtheta2_dt];
end

This program solves the 2nd order ODE representing the motion of a simple pendulum with damping.

The program is divided into 3 code sections:

1) Line(9 to 16) Defining the Problem Statement

2) Line(19 to 29) Solving the ODE for given inputs like mass, damping coefficient, length of string,

gravitational constant, time span and initial conditions.

3) Line(32 to 40) Generating the plot for "displacement vs time" and "Velocity vs time" of the

simple pendulum.

4) Line(43 to 65) Generating the animation for the motion of the pendulum with time for the interval

of 20 secs.

While solving the ode, ode45 command was used with the following syntax:

[t results]=ode45(@(t theta) pend_ode(M,B,L,g,theta,t), t_span,ini_cond);

here ode45 is the inbuilt function in matlab that solves a 1st order non stiff differential equaions.

since our equation is 2nd order differeential equation we convert it into two 1st order differential equations an solve them simultaneously.

the syntax above shows that it ode45 takes t and theta as independent and dependent parameters repectively. it takes another defined function "pend_ode", time limits and initial condition as inputs while giving 2 outputs: t and results.

'results' is a 2D array or matrix with column 1 representing displacements and column 2 representing velocity, at different time values, repsectively.

The animation is also uploaded on youtube which can be accessed by going to the following URL:

https://youtu.be/FTee3F0s8Ag