Solving second order ODE using MATLAB

Aim:

To solves the ODE which 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 a animated video of output obtain 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.Movement of the pendulum depends on the Newtons 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.

fig: Simple Pendulum (with FBD)

Equation of motion is:

In the above equation,

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

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

Let
θ=θ1

dθ1/dt=θ2

dθ^2/dt^2=dθ2/dt

Input :

• 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

Program:

program to simulate the transient behaviour of a simple pendulum using Matlab.

%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 program:

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
        

Detail explanation each step of program:

1.In first step Input and initial condition of simple pendulum is feeded.

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

3. In next step result of solved ODE is ploted w.r.t time using plot command.

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

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

6. In next step movie command is used to create animated video of oscilating pendulum.

Output:

1. The output obtained is graphical representation velocity curve and displacement curve w.r.t Time.

• Angular Displacement and Angular Velocity plots were generated for a time range of 20 secs.

2. The second output obtained is animated video of damped motion of Simple Pendulum w.r.t Time.

click on following link to watch the video:

Various Error faced during programing are follow:

1. Spelling mistake in writing function program.

• Step taken to over come this error:

By correcting spelling mistake from function program by reviewing each line of program.

2. Typing mistake in ploting Bob of pendulum.

• Step taken to over come this error:

By replacing "o" to "0" (Zero) in, plot(x_end,y_end,'marker','o','markersize', 28,'markerfacecolor','r') command.

3. The pendulum starts oscillating in inverted direction in animated video.

• Step taken to over come this error:

The final co-ordinate, y_end = 1*cos(theta_2 ) changed to y_end = -1*cos(theta_2 )

Conclusion:

The simulation of the transient behaviour of a simple pendulum is succesfuly created by using program in MATLAB.