Week 3 - Solving second order ODEs

clear all
close all
clc
% Input data such as damping cofficient, length, mass of ball

b=0.05; g=9.81; l=1; m=1;

theta_0 =[0;3]; % Intial condition i.e theta = 0, angular velocity = 3 rad/sec
t_span = linspace(0,20,500); % Taking range of 500 values between 0 and 20 seconds

[t, results]=ode45(@(t,theta) ode_func(t, theta, b,g,l,m),t_span,theta_0); % Using in-built MatLab ode solver

% Plotting the results
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('Plot')
xlabel('Time')

a=1;
for i =1:length(t)
x0=0; y0=0; % defining origin
x1=l*sin(results(i,1)); % defining position of ball
y1=-1*cos(results(i,1)); % defining position of ball

figure(2)
plot([-1 1],[0 0],'linewidth',4,'color','y'); % plotting for horizontal line
hold on
plot([x0 x1],[y0 y1],'linewidth',3,'color','r');
hold on
plot(x1,y1,'o','markers',30,'markerfacecolor','g'); % For marking ball position at given time 't'
grid on
axis([-2*l 2*l -2*l 2*l]); % Axis margin
pause(0.01)
hold off;
M(a) = getframe(gcf); % getting each frame
a =a+1;
end
movie(M)
videofile= VideoWriter('second_order_ode_pendulum_motion.avi','Uncompressed AVI'); % writing frame in opened file
open(videofile)
writeVideo(videofile,M)
close(videofile)

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

Animation of Pendulum motion drive link :

https://drive.google.com/file/d/1MQ04YQpystmPjnwpMFoxX3qhVS6zGmY5/view?usp=sharing

Plot verse time graph

Explanation of code

The input data such as damping coefficient, length of pendulum and mass of the ball is provided in the code. The initial condition of angular velocity is equal to 20 rad/sec at time equal 0 seconds is given in term of a column vector. The time span of 20 seconds is taken for the motion of a pendulum and 500 different values of time between 0 and 20 second is provided using linspace command. Now, ode_function is formed that take t, theta, mass, length, ‘g’, and value of damping coefficient as input. This function breaks given the second ode into a simple differential equation of a lesser degree. The variables thetha1 and thetha2 are created for this and value of thetha2 is the first derivative of thetha1. Now using ode45 solver this given second ode is solved. This ode45 solver takes ode_function and parameter on which given ode is dependent i.e ‘t’ and ‘theta’ in this case. It also takes time duration and initial condition as an input parameter. Using this solver, the given ode is solved and output is taken in variable ‘results’. Now the graph of plot verse time is drawn.

For creating an animation of pendulum motion, a ‘for’ loop is used so that for all time between 0 to 20s, we can get the position of the pendulum. Firstly the coordinate of origin is defined and taken at (0,0) point. Now the two plot command is used to plot the connecting length of the pendulum and marking ball of pendulum respectively. The range of the axis is given in the code. Movie of animation is created using in-build ‘movie’ and ‘videoWriter’ function.