Week 3 - Solving second order ODEs

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

b = damping

m = mass

g = gravity

l = length

`theta = theta1`

`d(theta)/dt = d(theta1)/dt = theta2` .......1

`d^2(theta)/dt^2 = d^2(theta1)/dt^2 = d/dt (d(theta1)/dt) = d/dt (theta2)` .......2

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

`d^2(theta1)/dt^2 + b/m d(theta1)/dt +g/l sin(theta1)`

`d(theta2)/dt +b/m (theta2) +g/l sin(theta1) = 0`

`d(theta2)/dt = - b/m(theta2) - g/l sin(theta1)`

STEP 1 : write ode function

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

STEP 2 :write a program that solve ode .

close all
clear 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);

plot(t, results(:,1))
hold on
plot(t,results(:,2))
ylabel('displacement')
xlabel('time')

STEP 3: FINAL PROGRAM

close all
clear all
clc

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

% Initial condition
theta_0 = [0;5];

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

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


% plotting
subplot(2,1,1);
plot(t,results(:,1),'Color','b')
hold on
plot(t,results(:,2),'Color','r')
hold off
subplot(2,1,2)
plot(results(:,1),results(:,2))
ylabel('position_radians')
xlabel('displacement & velocity')


ct = 1; 

for i =1: length(results(:,1))

    x0 = 0;
    y0 = 0;
    x1 = l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    figure(2)
    plot([-1 1],[0 0],'linewidth',6,'Color','r');
    axis([-2 2 -2 2]);
    line([x0 x1],[y0 y1],'linewidth',2,'color','g');
    hold on;
    hold off;
    M(ct) = getframe(gcf);
    ct = ct+1;

end

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