Week 3 - Solving second order ODEs

INTRODUCTION:

The 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.

AIM:

To write a program in Matlab based on the second-order differential equation and to simulate the pendulum motion,

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

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

FIXED PARAMETERS:

Lenth=1 metre,

mass=1 kg,

damping co-efficient=0.05.

acceleration due to gravity=9.81 m/s2.

motion has to be between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

FUNCTION CODE:

% ODE Function
function dtheta_dt = ode_fun(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

CODE EXPLANATION:

The function ' ode_fun ' is called for solving the ODE using the 'ode45' command. Using this function we get Time, Displacement, and Velocity values. A second-order ODE is converted in first order ODE. The general syntax of the 'ode45' is as follows :
[t,results] = ode45(@(t,theta) function_name(t,theta,b,g,l,m), integration_span, initial_conditions)

the conversion of second to first order was done by introducing the theta1 and theta2 variables.

MAIN CODE :

% Solving second order derivative
%  All numerical values in SI Units
clear all
close all
clc

% Fixed Parameters

% 1.Damping Co-efficein:
b = 0.05;
% 2.Acceleration due to gravity(m/s^2):
g = 9.81;
% 3.Mass(kg) :
m = 1;
% 4.Length(m) :
l = 1;

% Initial condition
theta_0 = [0;3];

% Time points
tspan = linspace(0,20,500);
% Solving ODE
[t,results] = ode45(@(t,theta) ode_fun(t,theta,b,g,l,m), tspan, theta_0);


% Plotting
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Time')
ylabel('Velocity (Orange) & 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)

plot([x0 x1],[y0 y1],'-o','MarkerSize',10,'MarkerFaceColor','r','MarkerEdgeColor','r','linewidth',2)
axis([-1.5 1.5 -1.5 1.5])
hold on
plot([-1 1],[0 0],'linewidth',4,'color','[0 0.6 0]')
M(ct) = getframe(gcf);
hold off
ct = ct + 1;
end

% Movie & Videofile
movie(M)
v_f = VideoWriter('Pendulum_Challenge.avi','Uncompressed AVI')
open(v_f)
writeVideo(v_f,M)
close(v_f)

 

Initially used commands - 'clear all' used to clear all the assigned values to the variables.
'close all'
'clc' - clears the command window
Fixed Parameters like acceleration due to gravity, length of the pendulum, the mass of the ball, damping coefficient are the inputs provided at the beginning.

Initial conditions which are displacement and velocity are given in the variable 'theta_0'. These conditions are further used to solve the ODE.

then the function program code was introduced in order to obtain the pendulum simulation.

Plotting the 'Velocity & displacement vs Time' graph. The 'plot' and the 'hold' commands are used to plot the graph.

PLOT:

DISPLACEMENT & VELOCITY VS TIME GRAPH:

PENDULUM MOTION:

ANIMATION: