Simple Pendulum_simulation_solving Second ODE

Aim: To Write a code that solves second order ODE Which reflects with equation of motion of a simple pendulum with damping.

To simulate the motion of the pendulum for 0-20 sec for angular displacement=0,angular velocity=3rad/sec^2 at time t=0.

To create an animation video of simulation for simple pendulum, subplots for the Angular displacement of bob with respect to angular velocity.

Given Data:

A simple Pendulum consists of a particle of mass attached to one end of a light inextensible string of length where the other end of the string is attached to afixed point.The particle when at rest hangs vertically below the fixed point.The particle and string when displaced from the equilibrium position oscillate in a circular arc in the same vertical plane.This physical description suggests the particle moves in 2D and therefore the motion will not behave like simple harmonic motion.

We consider a simple pendulum that consists of a mass m suspended by a string or rod of length L(and negligible mass).The angle it makes with the vertical varies with time as a sine or cosine.

Given Equation is:   

Looking at the forces on the pendulum bob,we see that the restoring force is proportional to sin theta, where as the restoring force for a spring is proportional to the displacement(which is theta in this case).

Where 

g=9.81;     % (m/s^2) Gravity.

m=1kg;       % (Mass of the bob).

l=1kg;         % (length of the String).

b= 0.05;      % (damping co-efficient).

CONVERSION OF SECOND ORDER ODE IN TO FIRST ORDER ODE:

A Second order ODE cannot be solved directly, so we have to convert this ODE into two first order ODE's.if it is first order simple ODE we can solve it using ananymous function.

Splitting Second Order ODE:

theta1=theta(1);

theta2=theta(2);

dtheta1/dt=theta2;

dtheta2/dt=-(b/m)*theta2 - (g/l)*sin(theta1)

Explanation to the code:

* Initially a function(ode_pendulum)is created to the call the ordinary differential equation of the simple pendulum.The code for function:

Function:

-for [dtheta_dt] = ode_dynamic_fun(t,theta,g,l,b,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.

* Then initial conditions of angular displacement and angular velocity are also given.

*Time period is declared using the 'linspace' function to get multiple time points.

*The ode_pendulum function is then called using ode45 a matlab solver for ode functions.

*Now under the forloop we define the co-ordinates(X0,Y0) and (X1,Y1) using the equations.

*The plot command is used to plot the links on the graph.And to animate the plot,getframe command is used(getframe returns amovie frame.The frame is a snapshots(pixmap)of the current axes or figure).And the counter variable is incremented.

*Both for loops are then ended and the syntax for movie is written.

Syntax for the Movie:

 movie(M)    %Plays the movie in matrix M once.

videofile=VideoWriter ('videoname.avi',file type)  %Assigns file name and file format.

open(videofile)   %Opens the videofile.

writevideo(videofile,M)  %Writes the video with plot frames in the array M.

close (videofile)   %closes the videofile.

Function Code:

function [dtheta_dt] = ode_dynamic_fun1(t,theta,g,l,b,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

Program Code:

% Simple Program for the Pendulum using the ode45.

clear all
close all
clc

% Basic Inputs

g = 9.81;  %(m/s^2) Gravity accelartion.
l = 1.0;   %(m)     Length of the Pendulum.
b = 0.05;  %(1/s)   Viscous Damping Co-efficient.
m = 1.0;   %(kg)    Mass of the Bob.
ct=1;

%Initial condition.
w_d = 0;
w_v = 5;
theta_0 = [w_d;w_v];
t_span = linspace(0,20,500);
dtmax=0.1;
ct=1;

[t,outputs] = ode45(@(t,theta)ode_dynamic_fun1(t,theta,g,l,b,m),t_span,theta_0);

wd_output1 = outputs(:,1);

for i=1:length(wd_output1);
    result = wd_output1(i); 
    
   % initial position of pendulum.
    x0=0;
    y0=0;
    x1= l*sin(result);
    y1= -l*cos(result);
    figure(1)
    subplot(2,4,[1,2,5,6]);
  
%Generate the Plot1 for the simulation of pendulum.
    plot(x0,y0,'k');
    hold on
    plot([x0 x1],[y0 y1],'color','b','linewidth',2)
    hold on
    plot(x1,y1,'k.','markersize',32,'markerfacecolor','g')
    hold off
    grid on
    axis equal
    axis([-1.5 1.5 -2 0]);  %axis[xmin xmax ymin ymax]
    title([num2str(t(i)),'s'])
    xlabel('x[m]')
    ylabel('y[m]')
    
 
  
    %Generate a subplot1 of the result.
    subplot(2,4,3);
    plot(t,outputs(:,1));
    hold on
    plot(t(i),outputs(i,1),'Ko');
    hold off
    xlabel ('Time(s)'); ylabel('Angle (rad)');
    title('Simple Damped Pendulum');
    
    %Generate a subplot2 of the result.
    subplot(2,4,4);
    plot(t,outputs(:,2));
    hold on
    plot(t(i),outputs(i,2),'Ko');
    hold off
    xlabel('Time(s)');
    ylabel('Rate (rad/s)');
    pause(dtmax);
    M(ct) = getframe(gcf);
    ct=ct+1;
    
end

movie(M);
videofile = VideoWriter('ODE_Simple_Pendulum.avi','Uncompressed avi');
open(videofile);
writeVideo(videofile,M);
close(videofile);

You tube link for an animation in your report:

https://youtu.be/870cGySDT9o