Solving second order ODE & generating Pendulum simulation using MATLAB

I) Code and Function Code :

Function code file name : 'ode_fun.m'

% ODE Function
function dthheta_dt = ode_fun(t,thheta,b,g,l,m);
thheta1 = thheta(1);
thheta2 = thheta(2);
dthheta1_dt = thheta2;
dthheta2_dt = -(b/m)*thheta2-(g/l)*sin(thheta1);
dthheta_dt = [ dthheta1_dt ; dthheta2_dt];
end

% Solving second order derivative 
% *********  All numerical values in SI Units  ************
 
clear all
close all
clc

% Fixed Parameters =========================>>>

% 1.Damping Co-eff.:
b = 0.05;
% 2.Acceleration due to gravity(m/s^2):
g = 9.81;
% 3.Mass(kg) :
m = 1;
% 4.Length(m) :
l = 1;

% Initial condition
thheta_0 = [0;3];

% Time points
tspan = linspace(0,20,500);

% Solving ODE =======================>>>
[t,results] = ode45(@(t,thheta) ode_fun(t,thheta,b,g,l,m), tspan, thheta_0);

% Plotting ================>>>
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Time')
ylabel('Velocity (Orange) & Displacement (Blue)')

% Animation =====================>>>
ct = 1;
for i = 1:length(results(:,1))
x0 = 0;
y0 = 0;
x1 = l*sin(results(i,1)+pi);
y1 = l*cos(results(i,1)+pi);
figure(2)

plot([x0 x1],[y0 y1],'-o','MarkerSize',10,'MarkerFaceColor','r','MarkerEdgeColor','r','linewidth',2)
axis([-1.5 1.5 -1.5 1.5])
hold on
plot([-1 1],[0 0],'linewidth',4,'color','[0 0.6 0]')
M(ct) = getframe(gcf);
hold off
ct = ct + 1;
end

% Movie & Videofile ===============>>>
movie(M)
v_f = VideoWriter('Pendulum_Challenge.avi','Uncompressed AVI')
open(v_f)
writeVideo(v_f,M)
close(v_f)

II) Outputs :

1. Velocity and Displacement vs Time graph :

The orange line represent velocity and the blue line shows the displacement.

2. Pendulum motion (Screenshot of animation) :

3.Video of the pendulum motion: https://youtu.be/tjfeM5TS6Tg

III) Stepwise explanation of code :

1. Initially used commands - 'clear all' used to clear all the assigned values to the variables.

'close all' - force deletes all figures

'clc' - clears the command window

2. Fixed Parameters like acceleration due to gravity, length of the pendulum, mass of the ball, damping coefficient are the inputs provided at the begining.

3. Initial conditions which are displacement and velocity are given in the variable 'thheta_0'. These conditions are further used to solve the ODE.

4. The next step includes giving the time points for the ODE. These points are used for integrating while solving the ODE. The 'linspace' command is used to divide the time interval into required parts.

5. The function ' ode_fun ' is callled for solving the ODE using the 'ode45' command. Using this function we get Time, Displacement and Velocity values. A second order ODE is converted in first order ODE. The general syntax of the 'ode45' is as follows :

[t,results] = ode45(@(t,thheta) function_name(t,thheta,b,g,l,m), integration_span, initial_conditions)

6. Plotting the 'Velocity & displacement vs Time' graph. The 'plot' and the 'hold' commands are used to plot the graph.

7. To create an animation for the pendulum a conter variable 'ct' which is equaled to 1 is created. Then a for loop has to be created from 1 to the length of time span and results array generated earlier.To do this the 'length' command is used. The following main code was used while plotting:

x0 = 0;
y0 = 0;
x1 = l*sin(results(i,1)+pi);
y1 = l*cos(results(i,1)+pi);

The loop also consists of the equation ' ct = ct +1' to increase the value of the counter variable.

8. The 'getframe' command is used to create a movie. For this the counter variable 'ct' is used to generate the movie array.

9. The created movie is used to create a video file using the 'VideoWriter' command. The 'open' and 'close' commands are used to create the open and close the file.Then 'writeVideo' is used to write the video into the file.The file is saved in the '.avi' format.

IV) Errors faced while coding :

1. Pendulum generated was facing upward:

Correction : Added 'pi' (+pi) to the angle which was responsible for generating the final point position of the pendulum (as discussed in the 7th step mentioned above) .

2. Didn't use the 'hold off' command:

Correction : Used the hold off command just before the end of the for loop.

V)Video Link: https://youtu.be/tjfeM5TS6Tg