Week 4.1 - Solving second order ODEs

Objective:

To solves the ODE which represents the equation of motion of a simple pendulum with damping and to animate the same.

Theory:

The motion of a simple pendulum is a basic classical example of simple harmonic motion, consisting of a small bob and a massless string. Movement of the pendulum depends on Newton's second law. When this law is written down, a second-order Ordinary Differential Equation is formed that describes the position of the end of the pendulum w.r.t time.

The differential equation y'' + p(x)y' + q(x)y = f(x) is called a second order differential equation with variable coefficients if the functions p(x) and q(x) are not constant functions and are functions of x.

EQUATION OF MOTION FOR SIMPLE PENDULUM

where,

• g = gravity in m/s^2.
• L = length of the pendulum in the meter.
• m = mass of the ball in kg.
• b = damping coefficient.

Since the given function is in the second-order differential equation we will convert it into first order to solve the ODE.

Detail explanation each step of the program:

1. In the first step Input and initial condition of a simple pendulum are fed.

b = 0.05;
g = 9.81;
l = 1;
m = 1;

2. Then the program is written to solve the ODE in which Matlab inbuilt function (ode45) is used with constant theta, t, and another function is used, which is defined separately in function program.

3. Solved ODE is plotted w.r.t time using plot command.

4. Now iteration of co-ordinate is fed using for loop to create animation.

5. Then plot command is used to plot pendulum string and bob.

6. The movie command is used to create an animated video of the oscillating pendulum.

Main code:

clear all
close all
clc


b = 0.05 % damping coefficient
g = 9.81 % gravity (m/s^2)
l = 1   % length of the spring (m)
m = 1  % mass of the bob (kg)

% initial condition
% angular displacement is 0 and angular velocity is 5

theta_0 = [0;5]

% time span
% Pendulum oscillates between 0-20 seconds

t_span = linspace(0,20,600)

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

%Plotting
plot(t,results(:,1))
hold on

plot(t,results(:,2))
title('Damped Harmonic motion of a Pendulum')
legend('Angular displacement', 'Angular velocity')
ylabel('Plot')
xlabel('Time(seconds)')



% Iteration using for loop to create animation of simple pendulum
ct=1;
theta_1 = results(:,1);

for i= 1:length(results(:,1))
    
    % initial co-ordinates
    x_start= 0;
    y_start= 0;
    theta_2 = theta_1(i);
    
    % final co-ordinates
    x_end = 1*sin(theta_2);
    y_end = -1*cos(theta_2);
    
    % ploting for animation
    figure(2)
    % base of pendulum
    plot([-1,1],[0,0], 'linewidth',4,'color','black')
    hold on
    % string of pendulum
    line([x_start x_end],[y_start y_end],'linewidth',6,'color','b')
    hold on
    % bob of pendulum
    plot(x_end,y_end,'marker','o','markersize', 28,'markerfacecolor','r')
    axis([-2 2 -2 2])
    pause(0.03)
    hold off
    M(ct)=getframe(gcf);
    ct=ct+1;
end

%saving animation file
movie(M);
videofile=videoWriter('Simple_pendulum.avi','MPEG-4');

Function code for ODE:

%Calling FUnction for solving second order differential
%ode function

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

theta1=theta(1)
theta2=theta(2)
dtheta1_dt=theta2;
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);   %converting second order diffrential to first order.
dtheta_dt=[dtheta1_dt;dtheta2_dt];
end

Plot screenshot:

Below is the graphical representation angular velocity curve and angular displacement curve w.r.t Time.

Youtube animation link:

https://youtu.be/v1jVubb33m0