Week 3 - Solving second order ODEs

https://drive.google.com/file/d/1gE2gv6FfOTgI6_ndVtg-9V62I2AMqKOm/view?usp=sharing

1. Main issue faced was with plot on Y axis. Finally realised Y axis is L*cos(dtheta/dt)
2. Unable to upload the animation file in the box below so shared the G-Drive link.

%Matlab program to simulate a simple Pendulum
clear all
close all
clc

b= 0.05; %Damping 
g= 9.81; %Gravitational force in m/s^2
l= 1;    %Length of the string in m
mass= 1;    %Mass of the object in Kg

%Initial Condition
theta_0 = [0;3]; %0 represents angular displacement starts at t=0 and 3 represents angular velocity at 3 rad/sec

%time points
t_span = linspace(0,20,500); % Creating an array to store time interval of 0-20 seconds with 500 values in it

%solve ode
[t,results] = ode45(@(t,theta) ode_func(t, theta, b,g,l,mass),t_span, theta_0);

%Created a counter variable to capture every frame and finally stich them
%into an animation
counter = 1;
%Started a for loop to extract position of pendulum w.r.t time
for i=1:length(t)
    %The loop will run for 500 times as the value of t is 500
    %Result plot has been made considering Displacement in the X axis and
    %change of length of string in the Y axis
    plot([0 results(i,2)], [4 cos(results(i,1))], '-o','MarkerSize',30,'linewidth',2);
    %Here plot([x0 x1],[y0 y1]) has been entered.
    %Displacement w.r.t to time is calculated in 1st column of results and L*cos(theta) will be value for Y
    %value for theta is calculated in function so cos of results give us
    %the final value
    %2nd column of results contain damping angular velocity w.r.t time
    xlabel('Displacement');
    %Based on the values we can creata an axis suitable for animation
    axis ([-4 4 0 5])
    pause(0.01)
    %using the counter variable to capture frames in an array
    M(counter) = getframe(gcf)
    counter = counter+1;
end
%creating an animation from the captured files
movie(M)
videofile = VideoWriter('ode_animation.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)


%Function to solve ODE to further use in simulating pendulum
function [dtheta_dt] = ode_func(t,theta,b,g,l,m)

%based on the formula creating variables and calculating value of theta
%w.r.t time and will further use it to pinpoint coordinated for Y axis
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