Week 3 - Solving second order ODEs

AIM:- To solve Second order ODE and animation of the pendulum using MATLAB.

OBJECTIVE:- To write code for the simulate the pendulum motion using the ODE function and creating an animation.

Simple Pendulum:-

A pendulum is a body suspended from a fixed support so that it swings freely back and forth under the influence of gravity. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging it back and forth.

ODE represents the equation of motion of a simple pendulum with damping

In the above equation,

g = 9.81 = gravity in m/s2,

L = 1 = length of the pendulum in m,

m = 1 = mass of the ball in kg,

b= 0.05 = damping coefficient.

Simulate the motion between

time = 0-20 sec,

angular displacement=0,

angular velocity=3 rad/sec at time t=0.

Error in line 20 & 21 are of not proper bracket place.

Line 18 explains the ODE function used in the function is ode45(@(t, theta) it's main inputs are t, theta & they are dependent and independent variable of the ode. It is defining the ode_func which has input parameter l, m, b & g, and after that, we provide time array (t_span) and initial condition to the starting of the time. In this, we are integrating the time into the initial condition.

function [dtheta_dt] = ode_func(t,theta,l,m,b,g)

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

clear all
close all
clc

% Inputs
b = 0.05; % damping co-efficient
g = 9.81; % gravitional force
l = 1; % length
m = 1; % mass

% Initial condition
theta_0 = [0;3]; %inital value of displacement and velocity.

% Time
t_span = linspace(0,20,500); % 0 sec to 20 sec with 500 value of time.

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

% PLotting
figure(1)
hold on
plot(results(:,1)); % Plotting the results of diplacement w.r.t time.
plot(results(:,2)); % Plotting the results of velocity w.r.t time.
hold off;
xlabel('plot'); % Nameing the x-axis.
ylabel('time'); % Naming the y-axis.
legend('displacement','velocity') % Naming the plotted line in the graph.

theta2 = results(:,1); % displacement of the pendulum
ct = 1;
for i = 1:100 
    
    theta3 = theta2(i);
    
     x0 = 0; %position of the fixed point of the pendulum.
     y0 = 0;
     
     x1 =  sin(theta3); % position of the another end which is oscillating at angel theta3.
     y1 = -cos(theta3);
     
     % plotting
     figure(ct)
     hold on;
     plot([-1 1],[0 0],'linewidth',3,'color','g'); % Fixed horizontal line where fixed end is  joint.
     plot([x0 x1],[y0 y1],'linewidth',3,'color','b'); % Moving Vertical line at theta3
     plot(x1,y1,'o','MarkerFaceColor','r','MarkerSize',10); % placing the bob at moving end of the vertical line.
     hold off;
     axis([-2 2 -1.5 1]); % Command for arrange of the axis.
     pause(0.03); % command for time duration between two frame.
     M(ct) = getframe(gcf); % To get all frame in one frame.
     ct = ct + 1;
end

% to make a movie or video for animation of the moving part.
movie(M)
videofile = VideoWriter('Simple Pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Youtube link for animation:- https://youtu.be/SMM1n04K6H8