Solving second order ODEs

Aim:To write a program that solves the ODE which represents the equation of motion of a simple pendulum with damping and also to simulate the pendulum motion.

Software Used:Octave

Theory

In Engineering, ODE is used to describe the transient behaviour of a system. An example for this is a simple pendulum. A simple pendulum is defined to have an object that has a small mass, also known as the pendulum bob, which is suspended from a light wire or string. When displaced to an initial angle and released, the pendulum will swing back and forth with periodic motion. The way the pendulum moves depends on the Newtons second law of motion.

A simple pendulum problem can be solved using a second order differential equation and the position of the bob with respect to time can be determined.

The second order differential equation for simple pendulum is given by;

Where, b=damping coefficient

m=mass of the bob

g=acceleration due to gravity

l=length of the string

Now, to solve the second order differential equation, we have to convert this 2^nd order differential equation into two first order ODE, as it cannot be solved directly. This is because we are making use of an inbuilt solver in matlab/octave called ode45 which solves only first order differential equations.

Conversion of 2^nd order ODE to 1^st order ODE

Function Code

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

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

Main Program

clear all
close all
clc
% simulation of simple pendulum
%inputs
b=0.05;  %damping coefficient
g=9.81;  %gravity in m/s^2
l=1;     %length of pendulum in m
m=1;     %mass of the ball in kg

%initial condition
theta_0=[0;3];

%time point
time=linspace(0,20,150);

%ODE solver
[t,results]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m),time,theta_0);
figure(1)
plot(t,results(:,1));
hold on
plot(t,results(:,2));
xlabel('Time(sec)')
ylabel('Plot')
legend('Displacement','Velocity')
grid on


ct=1;
for i=1:length(results(:,1))
  x0=0;
  y0=0;
  x1=l*sin(results(i,1));
  y1=-l*cos(results(i,1));
  
  figure(2)
  plot([-0.5,0.5],[0,0],'linewidth',5,'color','r');
  line([x0,x1],[y0,y1],'linewidth',2,'color','g');
  hold on
  plot(x1,y1,'o','markersize',20,'markerfacecolor','b');
  grid on;
  axis([-1.5 1.5 -1.5 0.5]);
  hold off;
  pause(0.03);
  ct=ct+1;
  text=sprintf("pendulum%d.png",ct);
  saveas(gca,text);
   
endfor

Explanation of code

1.Start the program with clear all, close all and clc

2.Specifying the inputs and Initial conditions of the simple pendulum

3.Using ode45 solver to solve the differential equations and calling the function ‘ode_func’ which contains the code to convert the 2^nd order differential equation to first order differential equation.

4.To plot the graph for ‘displacement vs time’ and also for ‘velocity vs time’ from the results of the differential equation. Hold on is for plotting all the cases in one graph.

5.The counter is initialized using ‘ct ’ command as it is used to capture the frame of motion of the pendulum at each value of time.

6.In order to plot the moment of the pendulum we need to know the positions of its ends at each time interval.

7.One end is fixed and the coordinates of the moving end are calculated i.e. a for loop is used to calculate the varying positions of the moving end (bob) for varying values of theta.

8.For the bob at the end , marker command is used to create a marker and ‘markersize’ and ‘markerfacecolor’ commands are used to obtain the desired size of the bob and its colour.

9.“saveas” and “sprintf” commands are used to save the plots for every iteration and ‘’imageMagik” is used to stitch the plots together and create animation.

Plot

Errors Obtained:Simple errors were obtained and were rectified easily.

YouTube link for animation :

https://youtu.be/H7-3lQPSgO4

Conclusion: Program that solves the ODE representing the equation of motion of a simple pendulum with damping was written and also the pendulum motion was simulated.