Solving second order ODEs for simulation of damped oscillation of a simple pendulum

---

---

Code for simple pendulum simulation:

% program for solving pedulum ODE
clear all
close all
clc

% inputs

b = 0.05;  % damping coefficient
g = 9.81;  % acceleration due to gravity
l = 1.0;   % Length of the string
m = 1.0;   % mass of the bob

% initial condition given that velocity at time 0 is 3 m/s 

initialcond = [0;3];  

% timespan for oscillation of pendulum to be calculated

time = linspace(0,20,500); 

% calling the ode 45 solver function for pendulum motion and angular velocity

[t,results] = ode45(@(t,theta)pendulum_ode(t,theta,b,m,g,l),time,initialcond);

f = 1; % to count the frames generated

%  simulating the pendulum motion from the results generated 

for i = 1:length(t)
    
    % to plot the free end of pendulum
    x = l*sin(results(i,1));
    y = l*cos(results(i,1));
     
    % plotting of string for pendulum 
    plot([0 x],[0 -y],'color','b','linewidth',2)
    hold on
    
    %plotting of the bob 
    plot(x,-y,'marker','o','markersize',20,'markerfacecolor','r')
    hold off
    grid on;
    
    % fixing axis
    axis([-1 1 -1.5 0]); 
    
    % display the time for simulation 
    title([num2str(t(i)),'sec']);
    xlabel('time');
    
    
    % to generate frames
    pause(0.05);
    
    
    % To store each frame generated into an array
    M(f) = getframe(gcf);
    
    % increasing array storage to store next frame
    f = f+1;
    
end

% plotting the graph of time vs angular displacement and angular
% velocity

figure(1);
plot(t,results(:,1));
hold on;
plot(t,results(:,2));
Xlabel('time');


% generating movie using the frames captured in M
movie(M)
videofile = VideoWriter('damped_simple_pendulum.avi','uncompressed.AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);

Explanation:

• initially we define the specifications like length, damping coefficient,mass and acc. due to gravity.
• we then provide the inital condition at time zero.
• then we determine the time-span for the simulation and divide it.
• then we call the ode45 solver to solve the second order differential equation which is broken down to 2 first order differential equations.
• we feed in the specifications along with time and initial condition which then gives us the values for position and angular velocity.
• now we create a simulation using the position of the free end obtained from the above result and taking its projections along x and y axis to mark the bob.
• we can also plot a graph using the results obtained

Pendulum_ode function

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

theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta2;
dtheta2_dt = (-b/m)*theta2 - (g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt;dtheta2_dt];
 

Function;

• the function takes in the specifications and also the initial conditions.
• the initial condition is store in theta array.
• we substitute the values of theta array into the first order differential equations and solve for timespan given.
• the results are then returned to the main function.

Plots for damped simple pendulum

• Simulation:

• Plot of time Vs Angular displacement (blue) and Angular velocity (red)

From the above plot we can find that the velocity and displacement of the pendulum is damped depending on the damping coefficient given.

---

Simulation along with plot

To generate simulation along with the plots we insert the following code in the for loop of the above main function.

for i = 1:length(t)
    
    % to plot the free end of pendulum
    x = l*sin(results(i,1));
    y = l*cos(results(i,1));
    
    
    subplot(2,4,[1:2,5:6]);
     
    % plotting of string for pendulum 
    plot([0 x],[0 -y],'color','b','linewidth',2)
    hold on
    
    %plotting of the bob 
    plot(x,-y,'marker','o','markersize',20,'markerfacecolor','r')
    hold off
    grid on;
    
    % fixing axis
    axis([-1 1 -1.5 0]); 
    
    %display the time for simulation 
    title([num2str(t(i)),'sec']);
    xlabel('time');
    
    pause(0.05);
    
    
    figure(1);
    clf;
    subplot(2,4,[3:4]);
    plot(t,results(:,1),'color','g');
    hold on;
    plot(t(i),results(i,1),'ko');
    hold off;
    xlabel('time');
    ylabel('Angular displacment (rad)');
    title('Simple pendulum damped');
    
    subplot(2,4,[7:8]);
    plot(t,results(:,2),'color','m');
    hold on;
    plot(t(i),results(i,2),'ko');
    hold off;
    xlabel('time')
    ylabel('Angular velocity (rad/s)');

• the code divides the outout window into a 2x4 grid.
• plots the simulation in the blocke1,2,5,6.
• while the plot for radians and angular velocity are plotted in blocks 3,4 and 7,8 respectivley.

---

ERRORS

The following error occured due to an invalid syntax