Solving Second Order ODEs

AIM

Simulation of the ,motion of a simple pendulum by solving the second order ODE and to create an animation of its motion.

OBJECTIVES

• To simulate the motion of the pendulum by solving the second order ODE for motion between 0-20 seconds for initial angular displacement = 0, and initial angular velocity =3 rad/s.
• To create an animation of the motion

THEORY

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

where

θ = Angular displacement (in radians)

b = Damping coefficient (in Nm/s)

m =Mass of the ball (in kg)

g =Acceleration due to gravity (in m/s²)

l =Length of the string (in m)

Steps to solve the second order ODE

1.Consider θ=θ₁

2.Hence

3.Now we get two first order ODEs

This vector on solving gives the value of θ₁ and θ₂

INPUT DATA

Initial angular displacement(θ at t=0) = 0,

Initial angular velocity (dθ/dt at t=0) = 3 rad/sec,

Damping coefficient, b =0.05 Nm/s,

Mass of the ball, m =1 kg,

Length of the string, l =1 m,

Acceleration due to gravity, g =9.81 m/s²

Time period t =0 to t =20 sec

CODE

Programming Language used- MATLAB

For solving the ODE and simulation

%Solve a second order ODE
clear all
close all
clc

%Input variables
b=0.05;%Damping coefficient(in Ns/m)
g=9.81;%acceleration due to gravity(in m/s)
l=1;%length of the string(in m)
m=1;%mass of the ball(in kg)

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

%Time span
time_span=linspace(0,20,100);

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

%Plotting the Angular displacement and Angular velocity
figure(1)
plot(time_span,results(:,1));
hold on
plot(time_span,results(:,2));
legend('Angular displacement','Angular velocity');
title('Plot for Angular displacement and Angular velocity');

%creating animation
ang_dis=results(:,1);
%counter
ct=1;
for i=1:length(ang_dis)
  theta_x=ang_dis(i);
  
  %starting coordinates
  x0=0;
  y0=0;
  
  %ending coordinates
  x1=l*sin(theta_x);
  y1=-l*cos(theta_x);
  
  %plotting displacement and velocity
  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,'marker','o','markersize',15,'markerfacecolor','r')
  hold off
  axis([-1 1 -1 1])
  grid on
  pause(0.05)
  file_text=sprintf('pendulum%d.png',ct)
  saveas(gca,file_text);
  ct=ct+1;
 endfor

For function ‘ode_func’

function [dtheta_dt]=ode_func(theta,t)
  %damping coefficient(in Ns/m)
  b=0.05;%Damping coefficient(in Ns/m)
  g=9.81;%acceleration due to gravity(in m/s)
  l=1;%length of the string(in m)
  m=1;%mass of the ball(in kg)
  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

Plot-

Animation of motion of Pendulum-

The animation can be viewed from the following link-

https://drive.google.com/file/d/100akvzIJLJjjYRlfWVka_6Y7ba8D2EWC/view?usp=sharing

EXPLANATION

Code-

The ‘ode_func’ is used to pass the two set of single order ODEs and the ‘lsode ‘fuction is used to solve the ODE using the initial conditions.

The plot is made by taking the position of pendulum at each value in time.The position is taken from the solution of the ODE i.e. the angular displacements.

Output-

The Angular Displacement and Angular Velocity plot is obtained. The animation shows the position of pendulum at 100 instances in between 0 to 20 seconds.

ERRORS ENCOUNTERED

This warning was displayed each time the position was plotted because of the error in syntax of the marker command.

plot(x1,y1,’o’,’markers’,15,’markerfacecolor’,’r’); was used[without ‘markersize’]

plot(x1,y1,'marker', 'o','markersize',15,'markerfacecolor','r'); was the correct command.

CONCLUSION

So, the second order ODE was solved and the animation of the motion of the simple pendulum was successfully created.