Week 3 - Solving second order ODEs

Function code:

function [dtheta_dt] = challenge_simple_pendulum_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

Main code:

clear all
close all
clc

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

% Initial conditions
theta_0 = [0;3];

%Time intervals
t_span = linspace(0,20,500);

% Calling ode function
[t,result] = ode45(@(t,theta)challenge_simple_pendulum_func(t,theta,b,g,l,m), t_span, theta_0);

%Plotting the results
figure(1);
plot(t,result(:,1), 'linewidth', 0.75,'color', 'r');
hold on
grid on
plot(t,result(:,2), 'linewidth', 0.75, 'color', 'b');
xlabel('Plot');
ylabel('Time');
axis([0 20 -3 3]);

% Creating animation of the pendulum

x0 = 1; y0 = 1; ct = 1;
theta = result(:,1);
for i = 1:length(theta)
figure(2);
x1 = 1-l*sin(theta(i));
y1 = 1 -l*cos(theta(i));
title('Motion of Simple Pendulum');
plot([0 2], [1 1], 'linewidth', 5, 'color', 'k');
hold on
plot([1 x1],[1 y1],[x1 x1],[y1 y1], 'o', 'markers', 20, 'markerfacecolor', 'r');
hold off
axis([0 2 -0.5 1.5]);
pause(0.02);
grid on;
M(ct) = getframe(gcf);
ct = ct+1;
end

% Making the animation into a video
movie(M);
videofile = VideoWriter('Simple pendulum.avi', 'Uncompressed AVI')
open(videofile);
writeVideo(videofile,M);
close(videofile);