Week 3 - 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.

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

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

Input:- 

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.

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Code:-

clear all
close all
clc

% Solving 2nd order ODE
b = 0.05;
g = 9.81;
l = 1;
m = 0.1;

%initial condition 
theta_0 = [0;3];

% time points
t_span = linspace(0,20,100);

% solve ODE
[t, result] = ode45(@(t,theta) ode_fun(t,theta,b,g,l,m), t_span, theta_0);

figure(1)
plot(t,result(:,1))
hold on
plot(t,result(:,2))
xlabel('plot')
ylabel('time')

ct = 1;
for i = 1:length(result(:,1))
x0=0;
y0=0;
x1=l*sind(result(i,1));
y1=-l*cosd(result(i,1));
%plotting
figure(2)
plot([-0.2,0.2],[0,0],'linewidth',8)
hold on
plot([x0,x1],[y0,y1],'linewidth',2,'color','r')
hold on
plot(x1,y1,'linewidth',8,'marker','o')
hold off 
axis([-1.5 1.5 -1.5 1.5])
pause(0.03)
grid on
M(ct) = getframe(gcf)
ct = ct+1;
end

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

Function:-

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

theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta2;
dtheta2_dt = -(b/m)*theta2 - (g/1)*sin(theta1);
dtheta_dt = [dtheta1_dt;dtheta2_dt];
end

Output:-

Links:-

https://www.youtube.com/watch?v=dXVbF8_qcJs