Week 3 - Solving second order ODEs

%inputs
clear all 
close all
clc

%length
L=1 ;
%mass
m=1 ;
%damping coefficient
b=0.05;
%gravity
g=9.81 

%initial 
theta_0=[0;3];
t_span = linspace(0,20,500);

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

%plot
plot(t, results(:,1))
hold on
plot(t, results(:,2))
xlabel('time')
ylabel('plot')

%create animation
ct =1 ;
theta_1 = results(:,1)

for i=1:length(results(:,1))
    theta_2 = theta_1(i);
    x0=0
    y0=0
    
    %pendulum side coordinates
    x1= 1*sin(theta_2);
    y1= -1*cos(theta_2);

    %plotting
    figure(2)
    plot([-1,1],[0,0],'LineWidth',4,'Color','r');
    hold on
    line([x0,x1],[y0,y1],'linewidth',3,'Color','g');
    hold on
    plot(x1,y1,'Marker','o','markersize', 20,'markerfaceColor','b');
    axis([-2 2 -2 2]);
    pause(0.03)
    hold off
    M(ct)= getframe(gcf);
    ct=ct+1
end

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

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)*sind(theta1)
         dtheta_dt = [dtheta1_dt; dtheta2_dt]
end

Graph-

Code explanation

determine initial conditions

use function code for dtheta_dt and using ode45

plotting the results of 20 seconds

determine start and end length of pendulum in the for loop for creating the animation

plotting the pendulum

using getframe to capture each frame

creating the animation in movie form.

youtube link-https://youtu.be/2PydxoBi6QI

errors faced

defined wrong variable name in function file, then corrected

space between holdon