Week 4.1 - Solving second order ODEs
AIM : To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping
In the above equation,
g = gravity in m/s2,
L = length of the pendulum in m,
m = mass of the ball in kg,
b=damping coefficient.
Input Values
L=1 metre,
m=0.1 kg,
b=0.05.
g=9.81 m/s2.
STEP1:- Convert 2nd degree ODE order into first degree ODE
Ɵ = Ɵ1
dƟ/dt = dƟ1/dt = Ɵ2
dƟ/dt+b*Ɵ/m+g*sinƟ1/L
STEP2:- Create a function of ODE
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
STEP3:- Create Simple pendumal program by usinf ODE function
-Initial condition are provided:- pi/2 is displacement and 0 is angular velocity
-Time interval 0 to 50 second in 200steps
-Solving ODE equation by using ode45 function
clear all
close all
clc
b = 0.05;
g = 9.81;
l = 1;
m = 0.1;
%initial condition
theta_0 = [0;5];
%time points
t_span = linspace(0,10,500);
%solve ODE
[t, results]= ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
ct =1;
angle = results(:,1)
for i = 1:length(angle)
ANGLE = angle(i);
x1= 5;
y1=5;
x2 = x1 +l*sin(ANGLE);
y2 = y1 -l*cos(ANGLE);
figure(1)
plot([4 6],[5 5])
hold on
plot([x1 x2],[y1 y2])
hold on
plot(x2,y2,'Marker','o','MarkerSize',20,'MarkerFaceColor','b')
hold off
axis([3 7 3 7])
M(ct) = getframe(gcf);
ct = ct+1;
end
videofile = VideoWriter('Simple_pendulum.avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)
figure(2)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('Plot');
xlabel('time');Plot :
Animation :
https://www.youtube.com/watch?v=OCaJACFuk3g
Errors :
Error faced to create video file and looping, solved using TSE.
References:
https://en.wikipedia.org/wiki/Pendulum