Week 3 - Solving second order ODEs

AIM :

Your objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

%inputs
l = 1 %meter
m = 1 %kg
b = 0.05 %damping coefficient
g = 9.81 %m/s2

%initial conditions
theta_0 = [0;3];

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

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

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

ct=1;
for i = 1:length(result(:,1))    
    x0 = 0;
    y0 = 0;
    
    x1 = l*sin(result(i,1));
    y1 = -l*cos(result(i,1));
        
    figure(2)           
    plot([x0,x1],[y0,y1],'r')
    hold on
    plot([1 -1],[0 0],'g','linewidth',3)
    plot(x1,y1,'marker','o','markersize',10,'markerfacecolor','b');
    hold off    
    axis([-1.5 1.5 -1.5 1.5])    
    M(ct) = getframe(gcf);   
    ct=ct+1;
end

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

To solve the ode by using ode45 we have to convert a second order differential equation to two first order differential equation.

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

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

results

Link - https://drive.google.com/file/d/14k3j3z0kWgP0pC0DU0INEN4NgspwDNvZ/view?usp=sharing