Solving second order ODEs Simple pendulum

Solving second order ODEs

Aim

Write a program to simulate the transient behaviour of a simple pendulum and to create an animation of it\'s motion.

Introduction

A pendulum is a weight suspended from a pivot so that it can swing freely. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the \"ball\" w.r.t time.

Governing equation

`(d^2theta)/(d t^2) + c/m*(d theta) / dt + g/l*sin(theta) = 0`

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

c=damping coefficient.

Code

%
clear all
close all
clc

%Inputs
%c = Damping coefficient, g = force of gravity, m = mass of bob, l = length
%of string
c = .05; 
g = 9.81;
m = 1;
l = 1;

%Initial condition
int = [0;3]; %initialisation condition
t_span = linspace(0,20,500); %time duration

%Solve ODE
[t,result] = ode45(@(t,result) ode(t,result,c,g,l,m),t_span,int);

%Plotting
plot(t,result(:,1))
hold on;
plot(t,result(:,2))
ylabel(\'Plot\')
xlabel(\'Time\')
legend(\'Displacement\',\'Velocity\')
hold off;

%Animation

%Defining coordinates
x0 = 0;
y0 = 0;
axis([0 2 0 2]);

ct = 1;
for i = 1:length(t_span)
    
    x1 = l*sin(result(i,1));
    y1 = -l*cos(result(i,1));
    
       %Plotting
        plot([x0,x1],[y0,y1],[-1,1],[0,0],\'linewidth\',2,\'color\',\'k\')
        hold on;
        plot(x1,y1,\'marker\',\'o\',\'markerfacecolor\',\'r\',\'markersize\',20)
        axis([-1 1 -1.5 0])
        pause(.005);
        hold off;
        
     M(ct) = getframe(gcf);
     ct = ct+1;
end

%Animation movie
movie(M)
videofile = VideoWriter(\'simple1_pendulum.avi\',\'Uncompressed AVI\');
open(videofile)
writeVideo(videofile,M)
close(videofile)
  

Function code

%Function to calculate displacement and velocity

function [dtheta_dt] = ode(t,theta,c,g,l,m)
theta1 = theta(1);
theta2 = theta(2);

%Breaking down the 2nd order ODE to two 1st order ODE

%dispalcement
dtheta1_dt = theta2;

%angular velocity
dtheta2_dt = -(c/m)*theta2 - (g/l)*sin(theta1);

%solution vector
dtheta_dt = [dtheta1_dt;dtheta2_dt];
end

Result

Link : https://youtu.be/bqOyUs_zwFo

Step by step explanation

The program start with commands clear all, close all and clc which clears commmand window and close all other previous results and windows.

Let `theta = theta1`

`(d theta)/dt = (d theta1)/dt = theta2`

`(d^2 theta)/dt^2 = (d^2 theta1)/dt^2 = (d theta2)/dt`

`(d^2theta)/(d t^2) = - (c/m*(d theta) / dt) - (g/l*sin(theta))`

`theta` = angle of pendulum at time t

`(d theta)/dt = w`(angular velocity of pendulum at time t)

Inputs are defined. Simulatig the motion between 0-20 sec, for angular displacement=0, angular velocity=3 rad/sec at time t=0. So initial conditions are inputed using the variable int. ode45 solver is used to solve this higher order differential equation. The marker command is used to create a sphere. One end of the pendulum is fixed and the other end moves according to the solution of the ODE you had solved.

Errors occured

Error 1: all images are occuring at the same time. So hold off function is used to animate properly.

Error 2: ode function failed to run. Here the ode syntax was written wrong. There are different ode solvers and ode45 can be used here.

The mistakes were easily found and corrected by reading the comments in the command window.

Conclusion

The transient behaviour of a simple pendulum was simulated and an animation of it\'s motion was created.