Week 3 - Solving second order ODEs

AIM: To write a program in octave to simulate the transient behaviour of a simple pendulum with damping and to create an animation of it's motion.The equation of motion is:

INTRODUCTION:

An ordinary differential equation is an equation containing an unknown function of one real or complex variable x. The unknown function is generally represented by a variable which therefore depends on x For example: dy/dx= f(x). ODE can also be used for showing transient behaviour in engineering applications just like the simple pendulum.

dθ/dt=dθ1/dt=θ2-eq1

d^2θ/dt^2=d^2θ1/dt^2=d/dt(dθ1/dt)=d(θ2)-eq2

d^2θ/dt^2+b/m.dθ/dt+g/lsin⁡θ=0

d^2θ/dt^2+b/m(dθ1/dt)+g/lsin⁡θ1=0

dθ2/dt+b/mθ2+g/lsin⁡θ1=0

dθ1/dt=θ2

dθ2/dt=-b/mθ2-g/lsin⁡θ1

d/dt[θ1θ2]=[θ2(-b/mθ2-g/lsin⁡θ1)].

The simple pendulum problem can be solved by using the second order ode functions and by solving these we will get the theta1 an theta2 values of the simple pendulum.

Simple Pendulum:

The simple pendulum consists of a point mass m hanging on a string of length L from some support. One pulls the mass back to some starting angle, θ and releases it, the motion observed can be studied through Newton's law of motion.

The acceleration of the simple pendulum will be given by,

a=-gsin⁡θ

The time period of the simple pendulum is,

T= 2π/ω

This ODE represents the equation of motion of a simple pendulum with damping.

Given Parameters:

g (m/s2)	gravity	9.81 m/s2
m (kg)	mass of the ball	1 kg
b	damping coefficient	0.05
L (meter)	The length of the pendulum	1 meter
Time Interval	Simulation motion	0-20 sec

angular displacement=0,

angular velocity=3 rad/sec at time t=0.

• SECOND ORDER DIFFERENTIAL EQUATION FOR SIMPLE PENDULUM :-

The position(theta) of the bob with respect to time of a simple pendulum can be determined by using the second order differential equations using the matlab. The value theta can be solved by using the second order ode and this can be achieved by breaking the higher order ode i.e., second order ode into no of first order ode.

The second order ode for simple pendulum will be as follows,

d^2θ/dt^2+b/m.dθ/dt+g/lsin⁡θ=0.

Where, b = damping coefficient

m = mass of the bob kg,

g = acceleration due to gravitym/s^2,

l = length of the string m.

MATLAB CODES FOR SOLVING SECOND ORDER ODE:

%To write a program to simulate the transient behaviour of a simple pendulum.
clear all
close all
clc
%Inputs
%Angular velocity
theta=[0;3];
%Time interval
t_span=linspace(0,20,100);
%Calculating ODE
[result]=lsode('pendulum',theta,t_span);
rad=result(:,1);
omega=result(:,2);
%Plotting the graph
figure(1)
plot(t_span,result(:,1),'color','r');
hold on;
grid on;
plot(t_span,result(:,2),'color','g');

hold off;
xlabel('Time interval');
ylabel('Angular velocity and Angular displacement');
legend('Angular velocity','Angular displacement');
title('Time Interval vs Angular Velocity and Angular Displacement');

ct=1;
for i=1:length(rad);
theta1=rad(i);
%Animation of pendulum
x0=0;
y0=0;
l=1;
x1=(1*sin(theta1));
y1=-(1*cos(theta1));

figure(2)
plot([x0 x1],[y0 y1],x1,y1,'g.','markersize',50,'linewidth',15);

axis([-1.5 1.5,-1.5 1.5]);
hold off;
grid on;
pause(1);
text=sprintf("pendulum%d.png",ct);
saveas(1,text);
ct=ct+1;
end

FUNCTION CODES:

%Function coding
function[dtheta_dt]=pendulum_motion(theta,t);
b=0.05;
m=1;
l=1;
g=9.81;
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

EXPLANATION OF CODES:

1) Use the clear all, close all, clc command .

(2) Then, we input the value of time span using 'linspace' command to get the linearly spaced vectors and input the value of angular velocity viz. linspace (0,20,100).

(3) Then we calculate the value of both dtheta1_dt and dtheta2_dt by calling a function code and use 'lsode' command to solve the ODE.

(4) Then we plot the graph between time interval and angular velocity,angular displacement using plot command and give label, legend and title.

(5) We use for loop for iteration and create different pictures of pendulum using [x0,y0],[x1,y1] and plot the figures for different iteration and use 'pause' command to pause the frames for every iteration.

(6) Pendulum bob is created using markersize and thread by linewidth command.

(7) Axes is setup using axis command.

(8) Using 'sprintf' and 'saveas' command we save the different figures of iteration of the plot and use imagemagick finally to stitch them all together to create an animation of the pendulum.

OUTPUT PLOT:

ERRORS ENCOUNTERED:

I have made an error in writing the formula for y1 by not putting the minus(-) sign.The pendulum bob was oscillating in the reverse direction which is not possible.This led me to get wrong output animation as shown below:

ERRORS RESOLVED:

The above mentioned error was resolved by entering the correct formula as:

y1= -(1*cos(theta1)).

The correct output animated screenshot is provided below:

RESULT:

The final output animated file can be viewed by clicking on below link:

https://youtu.be/7B_tW_yvhXg