Week 3 - Solving second order ODEs

AIM :- To write a code and solve the ODE which represents the equation of motion of the simple pendulum on MATLAB.

OBJECTIVE : The objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

`(d^2theta)/(dT^2) + (b)/(m)(d theta)/dt + (g)/(l) sintheta = 0`

EXPLANATON : A simple pendulum consists of a ball (point-mass) m hanging from a (massless) string of length L and fixed at a pivot point P. When displaced to an initial angle and released, the pendulum will swing back and forth with periodic motion. By applying Newton's second law for rotational systems, the equation of motion for the pendulum may be obtained

`τ=Iα`⇒ `−mgsinθL=(mL^2)(d^2θ)/(dt^2)`

and rearranged as

`(d^2θ)/(dt^2)+((g)/(L))sinθ=0` (without damping force)

`(d^2theta)/(dt^2) + (b)/(m)(d theta)/dt + (g)/(l) sintheta` (with damping force)

• The fundamental idea solving ODE is analytically and computational is when we have higher ODE ,we break down into set of first order ODE

`(d^2theta)/(dt^2) + (b)/(m)(d theta)/dt + (g)/(l) sintheta = 0`

Where θ = θ₁

`(d theta)/(dt^2) = (d theta1)/(dt) = θ2` ------------------------(1)

`(d^2 theta)/(d t^2) = (d^2 theta 1)/(d t^2) = (d)/(dt) ((d theta 1)/(d t)) = (d)/(dt) (theta 2)`------------------(2)

`(d^2theta)/(dt^2) + (b)/(m)(d theta)/dt + (g)/(l) sintheta = 0`

`(d^2theta 1)/(dt^2) + (b)/(m)(d theta 1)/dt + (g)/(l) sintheta 1 = 0`

`(d theta 2)/(dt) + (b)/(m)( theta 2) + (g)/(l) sintheta 1 = 0`

• We have created two ODE’S

1. `(d theta 1)/(d t) = theta 2`
2. `(d theta 2)/(dt) = - (b)/(m)( theta 2) - (g)/(l) sintheta 1`

This is breaking down into 2^nd order ODE into series of 1^st order ODE

`(d)/(dt) ((theta1),(theta 2)) = ((theta 2),(- (b)/(m)( theta 2) - (g)/(l) sintheta 1))`

INPUTS : This ODE represents the equation of motion of simple pendulum with damping

`(d^2theta)/(dT^2) + (b)/(m)(d theta)/dt + (g)/(l)sintheta = 0`

In above equation,

g = gravity in m/s2

L = length of the pendulum in m,

m = mass of the ball in kg,

b = damping coefficient.

• Inputs of the matlab program
• L = 1 metre,
• m = 1 kg,
• b = 0.05.
• g = 9.81 m/s2
• motion between 0-20sec,
• angular displacement = 0
• angular velocity = 3 rad/sec at time t = 0

PROGRAM :This code will shows the transient nature of pendulum

clear all
close all
clc

%inputs
b = 0.05;
g = 9.81;%gravity in m/s^2
l = 1;%length of the string in m
m = 1;%mass in kg

%intial condition, angular disp =0,angular velocity = 3rad/s
theta_0 = [0;3];

%time points
time = linspace(0,20,600);

%[t,Results]=ode45('challange_2nd_ode',t_span,theta_0);
% solve ODE
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),time,theta_0);

figure(1)
plot(t,results(:, 1))
hold on 
plot(t,results(:,2))
xlabel('Position-Radians');
ylabel('Angular Velocity');
legend('velocity vs time','displacement vs time')
grid on;

ERRORS :

UPDATED PROGRAM :

clear all
close all
clc

%inputs
b = 0.05;
g = 9.81;%gravity in m/s^2
l = 1;%length of the string in m
m = 1;%mass in kg

%intial condition, angular disp =0,angular velocity = 3rad/s
theta_0 = [0;5];

%time points
time = linspace(0,20,600);

%[t,Results]=ode45('challange_2nd_ode',t_span,theta_0);
% solve ODE
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),time,theta_0);

figure(1)
plot(t,results(:, 1))
hold on 
plot(t,results(:,2))
xlabel('Position-Radians');
ylabel('Angular Velocity');
legend('velocity vs time','displacement vs time')
grid on;

OUTPUT :The plot describes that as time increases the velocity and displacement of the SHM decreases since we have used the damping factor.

PROGRAM 2: This code which shows animation of motion of the pendulum.

clear all
close all
clc

%inputs
b = 0.05; %damping coefficient
g = 9.81;%gravity in m/s^2
l = 1;%length of the string in m
m = 1;%mass in kg

%intial condition angular disp=0,angular velocity=3rad/s
theta_0 = [0;3];

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

%[t,Results]=ode45('challange_2nd_ode',t_span,theta_0);
% solve ODE
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),time,theta_0);

figure(1)
plot(t,results(:, 1));
hold on 
plot(t,results(:,2));
xlabel('Position-Radians');
ylabel('Angular Velocity');
legend('velocity vs time','displacement vs time')
grid on;

%Simulation of pendulum
figure(2)
ct=1;
for i = 1:length(results(:,1))
    x0 = 0;
    y0 = 0; 
    x1 = l*sin(results(i,1)); 
    y1 = -l*cos(results(i,1));
    
    plot([-0.5,0.5],[0,0],'linewidth',5,'color','g');
    line([x0 x1],[y0 y1],'linewidth',2,'color','r');
    hold on;
   %plotting the bob using marker and facecolor   
    plot(x1,y1,'o','markers',15,'markerfacecolor','b');
    grid on ;
    axis([-2 2 -2 1])
    pause(0.00001)
    hold off;
    M(ct) = getframe(gcf); 
    %file name is pendulum and specify the location
    Text = sprintf("pendulum%d.png",ct);
    saveas(1,Text);%saveas(fig,filename)
    ct = ct + 1;
end

% creating  a movie
movie(M)
videofile = VideoWriter('pendulum_animation.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Function used in this program: Funtion is nothing but peace of code that u can use just by calling its name and passing in right parameters

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/1)*sin(theta1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];

end

OUTPUT :

https://youtu.be/K1b_J66RJvY