Program to solve second order ODE to simulate the motion of a simple pendulum

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

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.

But first we need to explain you what is a pendulum

Simple Pendulum - A simple pendulum consists of a mass hanging from a string of length and fixed at a pivot point .When displaced to an initial angle and released, the pendulum will swing back and forth with periodic motion.

we have written a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

`(d^2theta)/dt^2+b/m (d theta)/dt+g/Lsintheta=0`

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

To solve second order ODE we have to break it into two equations of first order and create a fuction file named as pendulum_damping_ode

Function code:

function thetadot = pendulum_damping_ode(theta,t,g,l,m,b)
  b= 0.05; %damping coefficient
  g= 9.81;  %gravity in m/s^2
  l=1;  %length of pendulum in m
  m=1;  %mass of ball in kg

  thetadot(1)= theta(2);
  thetadot(2)= (-b/m * theta(2))- (g/l *sin(theta(1)));
  endfunction

we have used thye values,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

And Simulated the motion between 0-20 sec,

angular displacement=0,

angular velocity=3 rad/sec at time t=0

Command lsode is used intead of ode in this program because this program is been written in octave instead of MATLAB and for the same reason the syntax for calling the function has also been changed which is clearly stated in comments in the code itself.

Saveas command is used to create animation by saving the images in the file named pendulum simulation and stitching the images afterwards with the help of imagemagick software.

Main code:

clear all 
close all
clc

%input conditions
theta_0=[0;3];

%time
t= linspace(0,20,500)

%solving ODE
%[t,result]= ode45(@(t,theta)ode_function(t.theta,b,g,l,m),t_span,t_theta_0) is the syntax used for function calling in MATLAB
[results]= lsode( 'pendulum_damping_ode',theta_0,t)

%plotting displacement and velocity
hold on
plot(t,results(:,1))
plot(t,results(:,2))
xlabel('time')
ylabel('plot')
title('Angular displacement and angular velocity wrt time')
legend('Angular displacement','angular velocity')
theta2= results(:,1);
ct=1;
for i=1:length(theta2)
  theta_x = theta2(i)
  x0=0;
  y0=0;
  x1= 1*sin(theta_x)
  y1= -1*cos(theta_x)
  
  %plotting pendulum
  figure(2)
  plot([-0.5 0.5],[0 0],'color','k','linewidth',3)
  hold on
  plot([x0 x1],[y0 y1],'color','b','linewidth',3)
  hold on
  plot(x1,y1,'o','markers',10,'markerfacecolor','r')
  hold off
  grid on
  axis([-2 2 -2 2]);
  
  %creating animation
  file_text = sprintf('pendulum%d.png',ct);
  saveas(gca, file_text)
  pause(0.001)
  ct=ct+1
  end

Graph for displacement velocity against time:

For animation you can check the following youtube link:

https://www.youtube.com/watch?v=76VOvB2FUzY

ERRORS:

Encountered an error 'vector lengths must match' because of not calling the function properly

Overcomed this error by calling the function correctly as function name was not stated properly and lsode command was not used correctly.