Solving a second order ODE using Octave

Aim:

To model a dampened simple pendulum using octave.

Description:

Any continuous physical system can be modeled approximately using a single or a system of differential equations for predicting its behavior. A dampened simple pendulum is one such example system. In this project this system is modeled using a 2^nd order differential equation :

$d^2\left(\Theta \right)\mathrm{/}d{t}^2+\left(b\mathrm{/}m\right)d\left(\Theta \right)\mathrm{/}dt+\left(g\mathrm{/}l\right)sin\Theta =0$d^2(\varTheta)/dt^2+(b/m)d(\varTheta)/dt+(g/l)sin\varTheta=0

Here,

$\Theta$\varThetaΘ - Angular Displacement,

b- Damping coefficient,

m- Mass of bob,

g- Acceleration due to gravity

l- Length of pendulum

Procedure for Solving the ODE :

For solving the ODE, it is first converted into two first order ODE's which are simultaneously solved to obtain a numerical solution.The 2^nd order ODE is converted by following substitution:

$\Theta =\Theta$\varTheta=\varTheta

Θ=Θ₁

d(Θ)/dt=
Θ₂

So, the 2^nd order ODE gets converted into two 1^st order ODE's:

d(Θ)/dt=
Θ₂

d(Θ₂)/dt=-(b/m)Θ₂-(g/l)sinΘ₁

Solving these equations numerically gives us angular position and angular velocity respectively.

Methodology :

1. The differential equations are defined in a separate function.
2. Now a new file is created in same folder.
3. In this new file initial the constants are defined.
4. Next a model of simple pendulum is made by using 'plot' and 'marker' command.
5. One end of simple pendulum is defined at origin and other end is defined as function of angular displacement.
6. The bob is made by marker command.
7. Next system of differential equations are solve using inbuilt ODE solver.
8. The results are used to obtain different plots which are stitched together to form the animation.

Code :

Main code :

clear all;

close all;

clc;

%defining the inputs or constants

l=2;

g=9.81;

damp_coeff=0.05;

mass=1;

%defining boundaries

theta_0=[0;3]; %initial condition of [position ; velocity]

t_0=linspace(1,20,500);%defining the Time domain

%calling the inbuilt ode solver

[t,output]=ode45(@(t,theta)diff_1(t,theta,l,g,damp_coeff,mass),t_0,theta_0);

%configuring the output

ang=output(:,1);

x_1=-0.5;%

x_2=0.5;

x1=0;

y1=0;

ct=0;

for i=1:length(ang)

y2=-l*cos(ang(i));

x2=l*sin(ang(i));

plot([x_1,x_2],[0 0],'linewidth',10,%creating base

[x1,x2],[y1,y2],'linewidth',3,'color','b',%creating string of pendulum

x1,y1,'marker','o','markersize',10,'markeredgecolor','g',%creating the holding ring

x2,y2,'marker','o','markersize',15,'markerfacecolor','r','markeredgecolor','r');%creating the bob

axis ([-5 5 -3.5 3.5])%defining the axis

pause(0.003);

%creating and saving an animation

text=sprintf('Output%d.png',ct);

saveas (gca,text);

ct=ct+1;

grid on;

end