Week 3 - Solving second order ODEs

PROBLEM STATEMENT:

Write a program in Matlab\Octave that will simulate the pendulum motion?

use,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

SOLUTION:

SIMPLE PENDULUM DEFINITION:

A simple pendulum is a mechanical arrangement that demonstrates periodic motion. The simple pendulum comprises a small bob of mass ‘m’ suspended by a thin string secured to a platform at its upper end of length L.

The simple pendulum is a mechanical system that sways or moves in an oscillatory motion. This motion occurs in a vertical plane and is mainly driven by the gravitational force. Interestingly, the bob that is suspended at the end of a thread very light somewhat we can say it is even massless. The period of a simple pendulum can be made extended by increasing the length string while taking the measurements from the point of suspension to the middle of the bob. However, it should be noted that if the mass of the bob is changed it will the period remains unchanged. Period is influenced mainly by the position of the pendulum in relation to Earth as the strength of gravitational field is not uniform everywhere.

Matlab code explanation:

1.

The first step in building the program is to clear all the memory, which might have been present earlier.

CLEAR removes all variables from the current workspace, releasing them from system memory.

2. The next step is to input the data as provided in the problem statement.

In the above statements,

b = damping coefficient

l = length of string

g = gravitational force experienced by the bob.

m = mass of bob.

3.

As per the problem statement the motion of the pendulum should be simulated between 0-20 seconds, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

In the above statement theta_0 will contain the values of angular displacement and angular velocity and t_span will contain the values of time interval between o to 20 seconds.

4.

In order to solve second order equation in matlab, we use ODE45 in bult function which is specifically used to such type of equations.

Syntax explanation:

[t,y] = ode45(odefun,tspan,y0), where tspan = [t0 tf], integrates the system of differential equations y′=f(t,y) from t0 to tf with initial conditions y0. Each row in the solution array y corresponds to a value returned in column vector t.

By using the above ODE45 function t will store all the time interval values and results will save all the angular displacement and angular velocity values.

5.

In the above function ode_fun is the function which needs to be defined by us to solve the desired system of equation. The below statements shows how the function is being defined.

 6.

In the above statements the results are plotted on the graph with the help of figure 1 command which will store results in figure 1.

7.

In the above statements, the one end of the string will be placed at (0,0) which will be fixed on the supporting bar, and the other end of the string will be determined by (x1,y1) positions.

In the next lines the plot and the line statements are used to show them on the graph at required places.

NOTE:

• In the next statements we have made use of marker command which is used to show a particular shape at required points. 
• Title command is used to give heading to a graph which would be shown at top center position of the figure.
• Legend command is used to specify the various plots or lines which are present in a given figure.

8.

In order to store the results in the movie format, the above statements were used.

COMPLETE PROGRAM:

clear all
close all
clc

b = 0.05;
g = 9.81;
l = 1;
m = 0.1;

%initial condition
theta_0 = [0 , 3];

%time span
t_span = linspace(0,20,200);

%solve ODE

[t,results] = ode45(@(t,theta) ode_fun(t, theta,b,g,l,m), t_span, theta_0);
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))

legend('Angular Velocity', 'Angular Displacement')
ylabel('plot')
xlabel('time')
grid on

ct = 1;

figure(2)
for i =  1:length(t_span)
x0 = 0;
y0 = 0;

x1 = l*sind(results(i,1));
y1 = -l*cosd(results(i,1));

plot( [-1 1],[0 0],'linewidth',7,'Color','r')
line([x0 x1], [y0 y1], 'Linewidth',3,'Color','g')
hold on

plot(x0, y0, 'marker', 'o','MarkerSize',5,'MarkerFaceColor','w')
plot(x1, y1, 'marker', 'o','MarkerSize',10,'MarkerFaceColor','b')

title('PENDULUM');
legend('support bar','supporting wire','starting point','bob');

axis([-2 2 -1.5 1]);
pause(0.003)
hold off

M(ct) = getframe(gcf);
ct = ct+1;

end


movie(M)
videofile = VideoWriter('pendulum.avi', 'Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

 RESULTS

ERROR INCURRED WHILE PROGRAMMING:

1.

The above error occured because the ode45 function inputs are not entered in correct manner, in the ode45 function,first the time span value needs to stored in each rows and then the theta vaues will be stored in each column.

Since I have entered theta first and then t in the function I received such error. 

The above error can be resolved by correctly mentioning the inputs at right places, which is shown below:

2.

In the above images we can see that we are not getting the desired result in which we can see both the angular velocity and angular displacement and also we are receiving a warning which tells us that the lengend command used is having inputs which are not present in the graph.

The above issue occured because matlab always stores value of the command which is last used, and since we have not written a command to show both the results in same fugure window we are not having the desired output.

The above issue is resolved by using an hold on command which is used to store both the plot inputs and show them in single figure window.

YOUTUBE LINK:

The below link shows the video where the simulation of the motion of the pendulum is shown.

https://youtu.be/9CEst5FiRe0

CONCLUSION:

1. Complete error free program written which shows the results of both angular velocity and angular displacement in a single figure window.

2. Simulation of the motion of the pendulum is also performed in the next figure window which shows the motion of the simple pendulum.