Week 3 - Solving second order ODEs

ANIMATION - https://drive.google.com/file/d/1LIjnCtDJdg13OaqYyQ3DeQsZ22eYMyxC/view?usp=sharing

THE SIMPLE PENDULUM

A simple pendulum consists of a mass m hanging from a 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 secont law for rotational systems, the equation of motion for the pendulum may be obtained

OBJECTIVE:

Is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

PARAMETERS:

g = gravity in m/s2 = 9.81

L = length of the pendulum in m = 1

m = mass of the ball in kg = 1

b=damping coefficient = 0.05

For time period = 0-20sec

Angular displacement = 0

Angular Velocity = 3 rad/sec

at time t=0

FUNCTION CODE:

% ODE 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

EXPLANATION

* In the above function, the second order ODE is simplified to first order ODE in terms of theta1 & theta2.

* theta1 & theta2 is differentiated and dtheta_dt gives a column matrix.

MAIN CODE:

clear all
close all
clc

% Damping coefficient
b=0.05;
% Acceleration due to gravity
g=9.81;
% Length
l=1;
% Mass
m=1;

% 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
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('time')
ylabel('Velocity(red) & Displacement(blue)')

% Animation
ct= 1;
for i = 1: length(results(:,1))
    x0 = 0;
    y0 = 0;
    x1 = l*sin(results(i,1)+pi);
    y1 = l*cos(results(i,1)+pi);
    figure(2)
    
    % Plotting coordinates
    plot([x0 x1],[y0 y1],'linewidth',3,'color','r');
    hold on;
    plot(x1,y1,'o','markers',15,'MarkerFaceColor','k');
    hold on;
    plot([-3 3],[0 0],'linewidth',3,'color','b');
    axis([-3 3 -3 3])
    M(ct) = getframe(gcf);
   hold off
   ct = ct+1;
    
end
% Movie
movie(M)
videofile = VideoWriter('simple_pendulum.avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)

EXPLANATION

* Fixed parameters are length, mass, damping coefficient and acceleration due to gravity.

* Time and angle are given as initial condition.

* T_span is 500 values between 0-20 time period

* ODE function is applied using command ode45

* Velocity & Displacement is plotted against time.

* Using loop command, the mass at various points and various times are plotted.

* Using movie command, animation is created.

PLOTS: