Week 3 - Solving second order ODEs

TITILE

Motion of a simple pendulum

OBJECTIVE

To calculate the angle and position of a simple pendulum and making an animation using.

EQUATIONS

`(d^2theta)/dt^2 + (b/m)*(d theta)/dt+(g/l)*sintheta =0`

where `theta` is the angle made by the pendulum with respect to an imaginary axis.

b = Damping coefficient

m = Mass of the pendulum

g = Gravity (9.81 m/s)

l = length of the pendulum in metres

CODE LOGIC

> We start off the code with 'clear all', 'close all', and 'clc'.

> Key in all the constant variables with the values.

> Give the initial conditions of the problem.

> This problem requires a time constrain (t_span), which uses a function (linspace()) to get the values of time from 0-20 seconds with 500 intervals for good resolution.

> We call in the ode45() function from which we have defined as below. This function recieves values of time,angle,b,g,l and m.

> We assign a name for the array of variables which have to be stored ([t,results]). The data will be stored in the array called 'results'.

> Now we plot the graph of angular velocity and angular displacement in a plot vs time graph. The values for which will be stored in the array 'results'.

> We use the 'legend()' function to label the respective curves.

> Then begins the for loop iteration for creating the animation.

> Before we start the for loop we give the desired part of the 'results' a variable, say in this code 'theta2'.

> ct=1 for the initial frame.

> We assign another variable 'THETA2' in order to help in the iteration with the iteration variable 'i'.

> x0,y0 are the points of the top of the pendulum.

> X1,Y1 are the coordinates with which we need a line to represent the length of the pendulum visially. Their values are w.r.t the x axis with angle 'THETA2'.

> Then we begine the code for the animation of the pendulum.

> We use 'figure(1)' and 'figure(2)' to differenciate the two results.

> We use the pause function with which rate we want the osillation to run.

> The syntax for making the movie is as follows.

movie(M) %plays the movie in matrix 'M' once
videofile = VideoWriter('pendulum.avi','uncompressed AVI'); %assigns a name and format for the video file
open(videofile) %opens the video file
writeVideo(videofile,M) %writes the video in said file
close(videofile) %closes the file

CODE

% solving second order ode using MATLAB

clear all
close all
clc
%inputs 

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

%initial conditions
theta_0 = [0;5];

% time points

t_span = linspace(0,20,500);

% Solving the ode

[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);


% plotting 
figure(1)

plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('time')
ylabel('plot')
legend('angular_dispacement','angular_velocity')

% for loop iteration
theta2 = results(:,1);
ct = 1;

for i = 1:length(theta2)
    THETA2 = theta2(i)


 x0 = 0;
y0 = 0;

X1 = l*sin(THETA2)
Y1 = -l*cos(THETA2)
% animation code
figure(2)
 plot([-1 1],[0 0],'linewidth',2,'color','b')
 axis([-2 2 -2.5 2.5])
 hold on
 plot(x0,y0,'s','markersize',15,'markerfacecolor','m')
 hold on
 plot([x0 X1],[y0 Y1],'linewidth',7,'color','y')
 plot(X1,Y1,'o','markersize',20,'markerfacecolor','k')
 grid on
 hold off
 
 pause(0.007)
 
 M(ct) = getframe(gcf);
 ct = ct+1
 
end

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

The code for definition the function

% ode function definition 

function [dtheta_dt] = ode_func(t,theta,b,g,l,m)

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

RESULTS

The below video shows the oscillation of a simple pendulum coded in MATLAB.

GRAPH

The above graph shows the relation between time and plot with angular velocity and angular displacement.

ERRORS

This error was caused because there was a ')' missing in the line of the code. The syntax was wrong. After realizing the same the parenthesis was put and was able to run the code.