Solving second order ODEs

Aim:

Write a program to simulate the transient behaviour of a simple pendulum and to create an animation of its motion.

Objective:

• Stimulate the motion between 0-20 seconds for angular displacement and angular velocity =3 rad/s at time t=0.
• Solve second order ODE of a simple pendulum with damping using lsode function in octave.
• Create an animation of the simple pendulum

Theory:

Consider a simple pendulum having mass ‘m’ attached to a string of length ‘l’ and having a damping coefficient ‘b’ is acting under the acceleration due to gravity ‘g’.

The ordinary differential equation for the motion of simple pendulum is given by :

in the given equation g=9.81m/s^2, l=1m, m=1 kg, b=0.05 N-s/m

Steps to solve the problem:

Now ,the code will be solved by using lsode function in ocatve

Code:

%To solve a second order ODE
clear all
close all
clc

%Input variables

%Damping coefficient unit-Ns/m
b=0.05;
%acceleration due to gravity unit-m/s^2
g=9.81;
%length of the string unit=m
l=1;
%mass of the ball unit=kg
m=1;

%Initial condition | angular displacement= 0  and angular velocity= 3 rad/s
theta_0=[0;3];

%time points
time_span=linspace(0,20,100);

%solving ODE
[results]=lsode('ode_func',theta_0,time_span);

%plotting
figure(1)
plot(time_span,results(:,1))
hold on
plot(time_span,results(:,2))

%creating animation
ang_dis=results(:,1);
%counter
ct=1;
for i=1:length(ang_dis)
  thetax=ang_dis(i);
 
  %starting coordinates
  x0=0;
  y0=0;
  
  %ending coordinates
  x1=l*sin(thetax);
  y1=-l*cos(thetax);
  
  %plotting
  figure(2)
  plot([-0.5 0.5],[0 0],'color','r','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
  axis([-1 1 -1 1])
  grid on
  pause(0.05)
  file_text=sprintf('output%d.png',ct)
  saveas(gca,file_text);
  ct=ct+1;
endfor

Code for ode_func:

function [dtheta_dt]=ode_func(theta,t)
  %Damping coefficient unit-Ns/m
  b=0.05;
  %acceleration due to gravity unit-m/s^2
  g=9.8;
  %length of the string unit=m
  l=1;
  %mass of the ball unit=kg
  m=1;
  theta1=theta(1);
  theta2=theta(2);
  dtheta1_dt=theta2;
  dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);
  dtheta_dt=[dtheta1_dt;dtheta2_dt];
endfunction

Output:

Error occurred:

The animation of pendulum was upside so, I had to change the y1 to negative y1 to bring it to down position.

Conclusion:

So, using the ordinary differential equation of simple pendulum and assuming angular velocity and displacement, an animation of simple pendulum was created.

Link for animation:

https://drive.google.com/open?id=1T1_RisjvMEx-ooCoBabxHZop1q9J7af4