Week 3 - Solving second order ODEs

% Matlab program to simulate the transient behaviour of simple pendulum

clear all
close all
clc

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

% Initial condition
theta_1 = [0;3];

% Time points
t_span = linspace(0,20,500);

% Solve ODE
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_1);

% Plottiong
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Time')
ylabel('Anagular displacement & Angular velocity')
legend('Angular displacement','Angular velocity')

% Animation
vv = 1;
for i = 1:length(results(:,1));
    % initial conditions
    x0 = 0;
    y0 = 0;
    % final conditions
    x1 = 1*sin(results(i,1));
    y1 = -1*cos(results(i,1));

    %plotting
    figure(2)
    %base pendulum
    plot([-1 1],[0,0],'linewidth',6,'Color','m')
    hold on
    %starting of pendulum
    line([x0,x1],[y0,y1],'linewidth',2,'color','b')
    %bob of pendulum
    plot(x1,y1,'marker','o','markersize',28,'markerfacecolor','r')
    axis([-2 2 -2 2]);
    pause(0.03);
    hold off
    M(vv) = getframe(gcf);
    vv = vv+1;
end

% Creating movie from captured freames
movie(M)
videofile = VideoWriter('simple_pendulum.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

>> Google drive link for the amimation plot: https://drive.google.com/file/d/1aaKLdMyc0WYwK3sz4WcDbH7LSQ6k5PUa/view?usp=sharing

>> The main objective is to solve the ODE and the ODE represents the equation of motion of a simple pendulum with damping, the code is to simulate the pendulum motion.

1. The input variables were defined in the program as given below
b = 0.05 % damping coeficient
g = 9.81 % gravity in m/s2
l = 1 % length of the pendulum in m
m = 1 % mass of the ball in kg
theta_1 = [0;3] % column matrix
t_span = linspace(0,20,500)
2. Created a fuction 'ode_func' to solve the ODE in matlab, the function commands are given below
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
3. The function 'ode_func' is being called in the main programm
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_1);
4. The figure(1) is plotted by using 'hold on and plot' commands of the matlab,
in the plot x_label represents 'Time' and the y_label represents 'Angular displacement & Angular velocity'
5. The 'for loop' is used to create the frames of motion of a simple pendulum, an animation for the simple pendulum plotted with the help of 'plot ans hold on' commands as given below
for i = 1:length(results(:,1));
% initial conditions
x0 = 0;
y0 = 0;
% final conditions
x1 = 1*sin(results(i,1));
y1 = -1*cos(results(i,1));

%plotting
figure(2)
%base pendulum
plot([-1 1],[0,0],'linewidth',6,'Color','m')
hold on
%starting of pendulum
line([x0,x1],[y0,y1],'linewidth',2,'color','b')
%bob of pendulum
plot(x1,y1,'marker','o','markersize',28,'markerfacecolor','r')
axis([-2 2 -2 2]);
pause(0.03);
hold off
M(vv) = getframe(gcf);
vv = vv+1;
end
6. The animation plot is saved as a movie animation with the help of below commands
movie(M)
videofile = VideoWriter('simple_pendulum.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)