Week 3 - Solving second order ODEs

clear all
close all
clc

% input values
b=0.05 
g=9.81
m=1
l=1

% initial condition
% initial displacement = 0
% angular velocity = 3 at theta=0
theta_0 = [0; 3]

% time points
% time period = 0 to 20 seconds
t_span = linspace(0, 20, 500);

% solve ODE
% Matlab code for ode function created
[t, results]= ode45(@(t, theta) ode_func(t, theta, b, g, l,m), t_span, theta_0);

% plotting graph
% plotting displacement
plot(t, results(:, 1))
hold on % to retain the existing plot

% plotting angular velocity
plot(t, results(:, 2))
ylabel('plot') 
xlabel('time')
legend('Displacement', 'Angular velocity')

ct=1
displacement=results(:,1) % stores displacement values in seperate matrix for loop
for i = 1:length(displacement)
    angle = displacement(i)
    figure(2)
    title('simple pendulum')

    % fix points for Pendulum
    x0 = 0
    y0 = 0
    x1 =1*sin(angle);
    y1 =-1*cos(angle);
    plot([-0.5, 0.5], [0, 0], 'linewidth', 2);% plotting fixed support
    hold on
    plot([x0, x1], [y0, y1], 'color', 'g', 'linewidth', 2) % plotting string
    plot(x1, y1, 'o', 'markers', 20, 'markerfacecolor','r') % plotting bob position
   
    % removing previous position
    hold off
    grid on
    axis([-1 1 -2 0.2])
    pause(0.01)
    M(ct)= getframe(gcf) % get current figure
    ct = ct+1 % represents every frame
end

% creating movie
    movie(M) % it will create movie using all the frames
    % create new AVI file where we can dump all these frames.
    videofile = VideoWriter('simple_pendulum.avi', 'uncompressed AVI');
    open(videofile)
    writeVideo(videofile,M)
    % M represents the frames that we have captured.
    close(videofile)

%ode function code
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];