Week 3 - Solving second order ODEs

clear all
close all
clc

b = 0.05; %damping coefficient
g = 9.81; %acceleration due to gravity (m/s^2)
l = 1; %length (m)
m = 1; %mass (kg)

%initial conditions
theta_0 = [0;3];

%time plots
t_span = linspace(0,20,1000);

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

%Plotting Angular displacement and Angular Velocity w.r.t time
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('time')
ylabel('Angular displacement (rad) & Angular Velocity (1/s)')
legend('Angular Displacement','Angular Velocity')
title('Motion of Pendulum')

ct=1;

%for loop to get coordinatesof pendulum for each value of displacement
for i = 1:length(results(:,1))
    ang_disp=results(i,1);
    %initial position of pendulum
    x0=0;
    y0=0;
    %new position of pendulum
    x1= l*sin(ang_disp);
    y1= -l*cos(ang_disp);
    
    %plotting positions of pendulum
    figure(2)
    plot([x0 x1],[y0 y1],'linewidth',1.5);
    hold on
    plot([-1 1],[0 0],'linewidth',3);
    plot(x1,y1,'o','MarkerSize',25,'MarkerFaceColor',[0 0.36 0.36]);
    hold off
    axis([-1.5 1.5 -1.5 0.3])
    title('Motion of simple pendulum')
    M(ct)=getframe(gcf);
    ct=ct+1;
end

%creating animation
movie(M);
videofile=VideoWriter('pendulum_animation.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile)

    

Line 1: Clears all variables from the workspace

Line 2: Closes all open matlab figure windows

Line 3: Clears command window

Lines 5, 6, 7 & 8: Assigning values for damping coefficient (b), acceleration due to gravity (g), length of pendulum (l) and mass (m)

Line 11: Declaring initial conditions for angular displacement and angular velocity

Line 14: Declaring time span for the pendulum motion. Here linspace command is used to get 1000 different values of t_span between 0 and 20 seconds

Line 16: Solves the second order differential equation for the motion of the pendulum and integrates the values from t=0 to t=20 (t_span values) with initial conditions given by theta_0. Each row in the output array 'results' correspons to a value of 'time_span'

NOTE: In the ode45 function, 2 anonymous variables, t and theta, are declared initially which stores the new values of t and theta for every iteration of the value of t_span

Line 19: Plotting output array 't' against the first column of output array 'results' (the first column of the output array 'results' are the values of Angular displacement at given 't')

Line 20: Retains plots

Line 21: Plotting output array 't' against the second column of output array 'results' (the second column of the output array 'results' are the values of Angular velocity at given 't')

Line 22: Labelling x axis as 'time'

Line 23: Labelling y axis as 'Angular displacement (rad) and Angular Velociy (1/s)'

Line 24: Adds legends to the plot for better understanding of the plot

Line 27: Initiates the variable 'ct'

Line 30: Initiates for loops for i values between 1 and the length of first coulmn of the results array

Line 31: Assigns the value of results(i,1) (each value of anglular displacement) to the variable 'time'

Line 33 & 34: Declares initial coordinates of the pendulum

Line 36 & 37: Calculating the final coordinates of pendulum for given angular displacement

Line 40: Creates new figure window

Line 41: Plots the string of the pendulum with a line width of 1.5

Line 42: Retains plot in current axes

Line 43: Plots the base of the pendulum with line width of 3

Line 44: Plots the pendulum ball with a size of 25

Line 45: hold command is turned off which means no new plot will be retained

Line 46: Sets axis limits

Line 47: Gives title to the plot

Line 48: Captures the data from the current frame and stores it in M(ct)

Line 49: Increments value of the variable ct by 1

Line 50: Ends for loop

Line 53: Plays the frames in array M

Line 54: Creates video file with the name 'pendulum_animation' as an uncompressed AVI file

Line 55: Opens the file associated with 'videofile' for writing

Line 56: Writes data from the vector M into the 'videofile'

Line 57: Closes the video file

Link for animation: https://www.youtube.com/watch?v=NT-T8T5N7PA

The function given within the ode45 command is as follows:

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

end

Line 1: Creates function with name 'ode_func' which accepts input variables 't,theta,b,g,l,m'

Line 3: Assigns the first value of the array 'theta' (which is the initial angular displacement) to the variable 'theta1'

Line 4: Assigns the secomd value of the array 'theta' (which is the initial angular velocity) to the variable 'theta2'

Line 5: Assigns theta2 to the variable 'dtheta1_dt' which is the final angular displacement

Line 6: Calculates dtheta2_dt value which is the value for final angular velocity

Line 7: Forms the vector 'dtheta_dt' which has the values of final angular displacement and final angular velocity in it

Line 9: Ends the function

NOTE: The plot attached explains the variation of both angular displacement and angular velocity with respect to time