Week 3 - Solving second order ODEs

AIM

To find out the plot of the second order ODE and stimulate the motion of a simple pendulum using Octave.

FORMULA

`d^2/dt^2theta+b/md/dttheta+g/lsintheta`

`theta=theta1`

`d/dttheta=d/dttheta1=theta2` (1)`

`d^2/dt^2theta=d^2/dt^2theta1=d/dt(d/dttheta1)=d/dttheta2` (2)`

`d^2/dt^2theta+b/md/dttheta+g/lsin theta=0`

`d^2/dt^2theta1+b/md/dttheta1+g/lsintheta1=0`

`d/dttheta2+b/mtheta2+g/lsintheta1=0`

where,

b=damping coefficient

m=mass in kg

g=gravity in `m/s^2`

l=length in m

THEORY

The ODE is an equation which can explain the motion of mechanical things which involve the springs and dampers.

A simple pendulum is an arrangement consisting of an string which is considered to be inextensible, having its length constant, and having a bob of a fixed mass at its end.

OBJECTIVE

1) To solve the second order ODE and plot its graph

2) To understand the motion of the simple pendulum and plot its animation using octave

CODE

Main code

clear all
close all
clc

 
%initial conditions
theta_0=[0,4];
l=1;
%time points
t=linspace(0,30,500);
%solve ODE
[results]=lsode('ode_function',theta_0,t);
%plotting
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('time')
ylabel('plot')
hold off


figure(2)
plot(results(:,1),results(:,2))
xlabel('angular displacement')
ylabel('angular velocity')

%simulation of simple pendulum
ct=1;
for i=1:length(results(:,1));
 theta3=results(:,1);
 THETA3=theta3(i);
 
 x0=0;
 y0=0;
 
 x1=l*sin(THETA3);
 y1=-l*cos(THETA3);
 
  % horizontal support for string
 plot([-0.5 0.5],[0 0],'linewidth',2,'color','k');      
 hold on
 grid on 
 
  % string
 plot([x0,x1],[y0,y1],'linewidth',1,'color','r');       
 hold on
 
  % bob tied at the end of string
 plot(x1,y1,'-o','markersize',30,'markerfacecolor','k');
 axis([-1.5 1.5 -1.5 1]);
 
   file_text = sprintf("armanimation%05d.png",ct)
      saveas(gca, file_text)
 pause(0.03);
 hold off
 
 ct=ct+1;
 
 end

Function

function[dtheta_dt]=ode_function(theta,t)
%inputs
b = 0.08;
g = 9.81;
l = 1;
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

ERRORS

1

2

HOW I SOLVED THIS ERRORS

I just edited the code and changed the name of the file name of the function code to ode_function, to get the proper output.

OUTPUT

The simple harmonic motion animation is here in this link

https://drive.google.com/file/d/1ltnNYc4JP5ASAakZkHAmu5Pe-QpxPDAW/view?usp=drivesdk

The graph of the 2nd order ODE is

CODE EXPLANATION

The code is written in octave

Firstly, write the 3 basic lines common to every octave code, ie, clear all, close all and clc.

Then we input the code comment wherever necessary using the % sign

Then write the function code to write the 2nd order ODE code with all the required data

The initial condition and the time points are given in the main code

The output from the function code is inputted in the main code

Then the graph is plotted

Now write the code for the simulation of the simple pendulum

Code is written for the horizontal support of the string

Code is written for the string

Code is written for the bob at the end of the string

Since the program is done in octave, the process is slightly different from that of MATLAB.

CONCLUSION

After rectifying all the errors and typing the code, we get all the required outputs, like the plot of the ODE and the Simple Pendulum animation.