Week 3 - Solving second order ODEs

To simulate

Simple pendulm motion .

Given

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient

Program

clear all
close all
clc

b = 0.5;
g = 9.81; - inputs
m = 1;
l = 1;

%initial condition

theta_0 = [0;3];

%time span

t_span = linspace (0,20,250)

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

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('time')
ylabel('theta') - plotting & specifications
legend('angular displacement','angular velocity')
hold off

ct = 1;
for i = 1:length(t)
x0=0;
y0=0;

x1=1*sin(results(i,1));
y1=-1*cos(results(i,1));

figure(2)

plot([x0 x1],[y0 y1],'color','r')
hold on
plot([1 -1],'color','b')

plot(x1,y1,'.',MarkerSize',50) - plotting & specifications
axis([-1 1 -1 1])
hold off

M(ct) = getframe(gcf);

ct = ct+1;

end

%making video file

movie(M)
Videofile = VideoWriter('simple_pendulum.avi')
open(videofile) - creation of video files
writeVideo(videofile,m)
close(videofile)

Program for ode_func

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

theta1=theta(1);

theta2=theta(2);

dtheta1_dt = theta2; - program for ode_func

dtheta2_dt = (b/m)*theta2-(g/l)*sin(theta1);

dtheta_dt = [dtheta1_dt ; dtheta2_dt];

end