Week 4.1 - Solving second order ODEs
Aim
Solving the 2nd order ode and simulating the motion of the simple pendulum
Simple pendulum
Simple pendulum is mechanical system that exhibits periodic motion. It consists of bob(weight) attached to the end of string or rod which is suspended from a fixed pivot joint. When displaced from its equilibrium position and released, the pendulum will swing back and forth under the influence of gravity.
The differential equation given is the differential equation for damped pendulum
The above given 2nd order ode is broken down into two
Substituting this in 2nd order ode
Then
Explanation
First, we provide the variables that are required
% Define parameters
b = 0.05; % Damping coeffeicent
g = 9.81; % Acceleration due to gravity
l = 1; % Length of rope
m = 0.15; % Mass of bob
theta_0 = [0;5]; % Initial condition
t_span = linspace (0,30,500); % Time range
The function named ode_func is called and ode45 command is used to compute the differential equation
% Solving ODE
[t, results] = ode45(@(t,theta) odefunc(t,theta,b,g,l,m),t_span,theta_0);
Then the resulting values are plotted to obtain angular displacement and angular velocity vs time graph
% Plotting results
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
hold off
title('2nd orde ODE angular displacement / angular velocity vs time')
xlabel('Time')
legend('Angular Displacement','Angular Velocity')
Position of simple pendulum at each instance are found using for loop and each frame of it are captured
ct = 1;
% Simulating the motion of simple pendulum
for i = 1:length(t)
x = 1*sin(results(i,1)); % x-coordinate of bob
y = -1*cos(results(i,1)); % y-coordinate of bob
figure(2)
plot([-0.5 0.5],[0 0],'LineWidth',1,'Color','k')
hold on
plot(x,y,'o','markersize',10,'Color','r')
plot([0 x],[0 y],'LineWidth',1,'Color','b')
plot(0,0,'.','Color','k')
xlabel('x axis')
ylabel('y axis')
title('Simple Pendulam Simulation')
legend('Fixed Line','Bob','Rope')
axis([-2 2 -2 2])
pause(0.03)
% Capture frames for the movie
M(ct) = getframe(gcf);
ct = ct + 1;
hold off
end
Then the frames are collected to form an animation then it is saved
% Making a movie from collected frames
movie(M)
% Creating and saving the file
videofile = VideoWriter('Simple_Pendulum_Simulation.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
Main program
clear all
close all
clc
% Define parameters
b = 0.05; % Damping coeffeicent
g = 9.81; % Acceleration due to gravity
l = 1; % Length of rope
m = 0.15; % Mass of bob
theta_0 = [0;5]; % Initial condition
t_span = linspace (0,30,500); % Time range
% Solving ODE
[t, results] = ode45(@(t,theta) odefunc(t,theta,b,g,l,m),t_span,theta_0);
% Plotting results
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
hold off
title('2nd orde ODE angular displacement / angular velocity vs time')
xlabel('Time')
legend('Angular Displacement','Angular Velocity')
ct = 1;
% Simulating the motion of simple pendulum
for i = 1:length(t)
x = 1*sin(results(i,1)); % x-coordinate of bob
y = -1*cos(results(i,1)); % y-coordinate of bob
figure(2)
plot([-0.5 0.5],[0 0],'LineWidth',1,'Color','k')
hold on
plot(x,y,'o','markersize',10,'Color','r')
plot([0 x],[0 y],'LineWidth',1,'Color','b')
plot(0,0,'.','Color','k')
xlabel('x axis')
ylabel('y axis')
title('Simple Pendulam Simulation')
legend('Fixed Line','Bob','Rope')
axis([-2 2 -2 2])
pause(0.03)
% Capture frames for the movie
M(ct) = getframe(gcf);
ct = ct + 1;
hold off
end
% Making a movie from collected frames
movie(M)
% Creating and saving the file
videofile = VideoWriter('Simple_Pendulum_Simulation.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
Program for breaking down differential equation
function[dtheta_dt] = odefunc(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
Conclusion
The 2nd order ode of simple pendulum is solved and its motion is also animated. The animated video is available in the link shown below.
https://youtu.be/mOq13dP0gvE?si=7cV6XYDnSFxpH10M
