Week 3 - Solving second order ODEs

SECOND ORDER ODEs

AIM: To solve an ODE that represents the equation of motion of a simple pendulum with damping.

OBJECTIVE: To write a program that solves the following ODE which represents the equation of motion of a simple pendulum with damping and create an animated video of the output obtained by solving this ODE.

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. The 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 equation of Motion is:

Since the given function is in the second-order differential equation, we have to convert it into first order to solve the ODE. The simplification is done as shown:

INPUTS:

• L = 1 metre
• m = 1 kg
• b = 0.05
• g = 9.81 m/s^2
• Time span = 0-20 sec
• Angular displacement = 0
• Angular velocity = 3 rad/sec2 at time t = 0

CODE:

%program to simulate the transient behaviour of a simple pendulum
clear all
close all
clc

%inputs
b = 0.05;
g = 9.81;
l = 1;
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);

% ploting
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('Plot')
xlabel('Time in sec')

% 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
% creating movie from captured frames
 movie(M)
 videofile = VideoWriter('simple_pendulum.avi','uncompressed AVI');
 open(videofile);
 writeVideo(videofile,M);
 close(videofile);

FUNCTION CODE:

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:

1. In the first step Input and initial condition of a simple pendulum are fed.
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 the function program.
3. In the next step result of the solved ODE is plotted w.r.t time using the plot command.
4. Now iteration of co-ordinate is fed using for loop to create animation.
5. Then plot command is used to plot the pendulum string and bob.
6. In the next step, the move command is used to create an animated video of the oscillating pendulum.

OUTPUT:

• The output obtained is a graphical representation velocity curve and displacement curve w.r.t time.
  
• Animation (YouTube link) – https://youtu.be/0Juk2UM0bzM

ERRORS:

Function name was mistakenly written as ‘ode_function’ instead of ‘ode_func’ which is then corrected and resolved this error.

The first letter of ‘theta_1’ was small which is mistakenly written as capital and showed an error.

CONCLUSION: The simulation of the transient behavior of a simple pendulum is successfully created by using a program in MATLAB.