Week 3 - Solving second order ODEs

b = 0.8; damping coefficient
g = 9.81; acceleration due to gravity
l = 1;  length of the pendulum
m = 1.5; mass of the bob

now the initial conditions

%initial conditions   initial displacement =0 and velocity =5
theta_0 = [0;5];

%time_points
t_span = linspace(0,10,200); time interval between 0 to 10 sec

now calling the ode function

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

to solve this we have to use a function

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

theta1 = theta(1)
theta2 = theta(2)
dtheta1_dt = theta2;               this is our first ode
dtheta2_dt = -(b/m)*theta2-(g/l)*sin(theta1);   this is our second ode
dtheta_dt = [dtheta1_dt;dtheta2_dt];
end

now plotting

plot(t,results(:,1))
hold on
plot(t,results(:,2))
disp = results(:,1)

the curve showing displacement and velocity w.r.t time

red curve is the velocity and the blue one is the displacement

now to plot the pendulum

figure(2)
x0 = 0 ; these are the initial points
y0 = 0 ;
ct = 1;

now to get the motion of the pendulum

for i=1:length(disp)
    x1=l*sin(results(i,1));  these are the final points
    y1=-l*cos(results(i,1));
     plot([x0 x1],[y0 y1]);  joining the two points
    hold on
    plot(x0,y0,'s','Markersize',10)  to create a box at the starting point
    plot(x1,y1,'.','Markersize',50)  to create a bob at the ending point

    
     axis([-2 2 -2 2]);
   
    hold off
    M(ct)=getframe(gcf);
    ct=ct+1;
     
    
end

to create this into a video

movie(M)
videofile=VideoWriter('pendulum','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

the link of the video is given below

https://drive.google.com/file/d/1dtWTNROHJD72NU-WXJMoFwDeuTCpbFuq/view?usp=sharing