Week 3 - Solving second order ODEs

1.AIM: A program in Octave that will simulate the pendulum motion.

2.GOVERNING EQUATIONS IF ANY: `(d^2 θ)/(dt^2 )+b/m (dθ/dt)+g/L sinθ=0`

This equation represents the motion of a damped pendulum.

3.OBJECTIVES: A programe that solves a second order differential equation. The ODE represents the equation of motion of a simple pendulum with damping. The programe plots graph of angular velocity and angular displacement with time. The programe also simulates the motion of the pendulum with damping.

4.BODY OF THE CONTENT: This code is written in Octave.

Main code:

clear all
close all
clc

%Inputs
%damping coefficient
b = 0.05

%Acceleration due to gravity (m sec^-2)
g = 9.81

%Length of the pendulum (m)
l = 1

%Mass of pendulum ball (kg)
m = 1

%Initial condition (angular displacement,angular velocity)
theta_0 = [0;3]

%time period (divided into 100 times)
t_span = linspace(0,20,100)

%Solving ode
[results] = lsode('ode_function',theta_0,t_span);

figure(1)
hold on
plot(t_span,results(:,1),'color','b','linewidth',2)
plot(t_span,results(:,2),'color','g','linewidth',2)

%Naming X axis
xlabel('Time')
%Naming Y axis
ylabel('Angular velocity & Displacement')

%Naming both the plots
legend('Displacement','Angular velocity')

%Result of the ODE is given a variable
theta = results(:,1)

ct = 1
for i = 1:length(theta)
    T = theta(i)
    x0 = 0
    y0 = 0
    x1 = l*sin(T)
    y1 = -(l*cos(T))
    
    figure(2)

%For plotting of surface line
    plot([-1 1],[0 0],'color','g','linewidth',2)
    hold on

%For point of contact of thread and surface
    plot([0 0],[0 0],'o','markerfacecolor','r','markersize',20)

%For thread representation    
    hold on
    plot([x0 x1],[y0 y1],'color','b','linewidth',4)

%For ball representation    
    hold on
    plot(x1,y1,'o','markerfacecolor','r','markersize',50)
    hold off

%Limit of axes is defined    
    axis([-2 2 -2 2])
    grid on
    
%Pause is given after every plotting   
    pause(0.00001)
%Images are saved    
    file_text = sprintf("Damped_pendulum%001d.png",ct)  
    saveas(gca,file_text)  
    
    ct = ct+1
end

Function code:

function [dtheta_dt] = ode_function(theta,t)
%Inputs
%damping coefficient
  b = 0.05

%Acceleration due to gravity (m sec^-2)
  g = 9.81

%Length of the pendulum (m)
  l = 1

%Mass of pendulum ball (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]
end

5.PRESENTATION OF THE CODE:

Line 1-3: Clear the screen

Line 7: Damping coefficient of the pendulum which is an input.

Line 10: Acceleration due to gravity which is an input.

Line 13: Length of the pendulum thread which is an input.

Line 16: Mass of the pendulum ball which is an input.

Line 19: theta_0 is the initial condition which represents the angular displacement(=0) angular velocity (= 3 rad/sec)

Line 22: t_span is the range of time period (0 to 20 seconds) and divided into 100 times

Line 27: Figure command is used to plot the graphs separately

Line 28: hold on command is used to plot the plots in same graph with existing properties of axes.

Line 29: plot command is use to plot between time period and angular displacement

Line 30: plot command is use to plot between time period and angular velocity

Line 33: xlabel command is use to naming of X axis

Line 35: ylabel command is use to naming of Y axis

Line 38: legend command is used for giving name of both the plots with showing of plot colour.

Line 41: theta variable is used to store the result of the ODE

Line 44: for loop starts from here.It starts from i = 1 to length of theta.

Line 45: T is another variable used to store the different values of theta at different values of i.

Line 46-49: x and y coordinates are defined.

Line 51: figure command is used to plot the simulation.

Line 54: Surface line is plot with coordinates (-1,0) and (1,0). Green colour is given for the line with 2 unit of linewidth.

Line 58: The point of contact of the thread with the surface is plotted with coordinates (0,0) with a marker of red colour having 20 units.

Line 62: Thread is plotted with coordinates (x0,y) and (x1,y1). Blue color is given for the line with 4 units of linewidth.

Line 66: Ball is plotted with coordinate (x1,y1) with a red color marker having 50 units size.

Line 70: Axes limits are defined as suitable.

Line 73: Pause of 0.0001 sec is given between consecutive plots.

Line 75: Image name is defined.

Line 76: Images are saved in the directory.

Line 79: Ends the for loop.

Ode_function code

Line 2: Name of the function is defined.

Line 5-14: Inputs of the function are given.

Line 16-20: Function is defined according to the differential equation.

Line 21: Ends the function

5.ERRORS AND DEBUGGING:

1. I did not put the inputs in ode_function programe then I got the following error.

2. I did not put the hold off command at the end of the plots then I got the following type of plot.

Youtube video link of the damped pendulum simulation: https://youtu.be/13dMCjNgJYY

6.OUTPUT:

7.CONCLUSION: We have got a simulation of damped pendulum. Also I have obtained the graphs for angular displacement and angular velocity. After the correction of the errors I got the simulation correctly.