Solving Pendulum 2nd order ODE

AIM 

 The aim of the project is to create MATLAB program to solve the ODE which represent the equation of the motion of the simple pendullum with damping.

Second order differential equation :

In the above equation 

g = gravity in m/s^2

L = length of the pendulum in m

m = mass of the ball in kg

b = damping coefficient 

Given values :

L = 1m

m = 1kg

b = 0.05

g = 9.81 m/s^2

angular velocity = 3 rad/s^2

angular displacement = 0 rad/s

Program :

The inputs for the pendulum is provided & function to calculate the theta values ode45 to calculate the value at different time step

ode45 Solve non-stiff differential equations, medium order method. [TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the system of differential equations y\' = f(t,y) from time T0 to TFINAL with initial conditions Y0. ODEFUN is a function handle. For a scalar T and a vector Y, ODEFUN(T,Y) must return a column vector corresponding to f(t,y). Each row in the solution array YOUT corresponds to a time
returned in the column vector TOUT. To obtain solutions at specific times T0,T1,...,TFINAL (all increasing or all decreasing), use TSPAN = [T0 T1 ... TFINAL].

The getframe command is used to store the plot of each step & last four lines are used to create a video file.

% Pendulum Program

% Inputs

b=0.5; % damping coefficient
g=9.81; % acceleration due to gravity
l=1; % length of pendulum
m=1; % mass of pendulum

% Initial condition
theta_0=[0;3];

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

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

% Graphs

plot(t, r(:,1))
hold on 
plot(t,r(:,2))
legend(\'displacement\', \'velocity\')

%Animation

    c= 1;
for i = 1:length(r(:,1))

    x0 = 0;
    y0 = 0;
    x1 = sin(r(i,1));
    y1 = -l*cos(r(i,1));

    figure(2)
    plot([x0,x1],[y0,y1],\'linewidth\',2)
    axis([-2 2 -1.5 0.5]);
    viscircles ([0 0], 0.02); % hook
    viscircles([x1 y1],0.20); % pendulum ball
    pause(0.001);
    Movie(c) = getframe(gcf); 
    c= c + 1;
 
end

movie(Movie)
videofile = VideoWriter(\'Pendulum_Movie\',\'Uncompressed AVI\');
open(videofile)
writeVideo(videofile,Movie)
close(videofile)

FUNCTION :

The function is used to provide the output of the value theta along the x & y axes.

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

PLOT :

The plot which depicts the value of the angular displacement & angular velocity of the pendulum.

Video File of Pendulum oscilation :