Week 3 - Solving second order ODEs

Aim:-

To simulate a simple pendulum by solving second-order Ordinary Differential Equation with the help of MATLAB.

Objectives:-

• Solve the second-order Ordinary Differential Equation.
• to create an animation file of a pendulum.
• to simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

Solving Second-order Ordinary Differential Equation:-

The given equation represents the equation of a simple pendulum.

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

To convert this second order Ordinary Differential Equation to first order Ordinary Differential Equation.

let,

`theta` =`theta1`

`(d theta)/dt=[d theta1]/dt=theta2----- -(1)`

`(d^2 theta)/dt^2=(d^2 theta1)/dt^2=d/dt((d theta1)/dt)=d/dt (theta2)-----(2)`

`(d^2 theta)/dt^2+b/m*(d theta)/dt+g/l*sin*(theta)=0`

`(d^2 theta1)/dt^2+b/m*(d theta1)/dt+g/l*sin (theta1)`

`(d theta2)/dt+b/m theta2+g/l*sin (theta1)=0`

`(d theta2)/dt = -b/m*theta2-g/l*sin(theta1)`

After solving this equation we can solve this equation in MATLAB, the programing to solve this equation of pendulum is given below.

MATLAB Program:-

For Function equation:-

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

Explanation:-

• This program is used to create the values dtheta_dt in matrix form so that we can use these values to calculate the values of angular velocity and angular displacement.
• Here we use the function command to find the values of dtheta _dt with the help of formula [dtheta_dt] = ode_func(t,theta,b,g,l,m) .
• We can call this code to our main code to find the values of angular displacement and angular velocity with the help of this function command.
• dtheta1_dt = theta2; dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1); dtheta_dt = [dtheta1_dt; dtheta2_dt]; this formulas can be used to get the values of dtheta_dt.
• After this program, we wrote the main program code which is given below.

For Pendulum:-

clear all
close all
clc

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

% intial condition
theta_0 =[0;3];

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

% solve ODE
[t, result] = ode45(@(t,theta) ode_func(t,theta, b,g,l,m),t_span, theta_0);

% ploting and labeling
figure(1)
plot(t,result(:,1))
hold on
plot(t,result(:,2))
xlabel('time')
ylabel('plot')
title('pendulum_displacemnt')

% animation of pendulum

ct=1;
for i =1:length(t)
    x=sin(result(i,1))*l;
    y=-cos(result(i,1))*l;
    
    figure(2)
    plot([-1,1],[0,0],'linewidth',5) 
    hold on
    plot([0,x],[0,y],'linewidth',5)
    plot(x,y,'Marker','o','Markersize',25,'Markerfacecolor','r')
    hold off
    axis([-1.5,1.5,-2,2]);
    title('pendulum')
    M(ct) = getframe(gcf);
        
    ct=ct+1;
end

% to make a movie of a pendulum

movie(M)
videofile = VideoWriter('pendulum_movement','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Explanation:-

• First, we use the clear all, close all, and clc commands to clear all the unwanted data in the workspace and the command window and close background plots.
• After that, we provide all the necessary input data to the program which is

  L=1 meter,

  m=1 kg,

  b=0.05.

  g=9.81 m/s2

• Also, we provided the initial condition of angular displacement=0, angular velocity=3 rad/sec at time t=0 to simulate the motion between 0-20 sec.
• While providing input data, we used linspace to generate an array of t_span from 0 to 20 sec of 500.
• To solve the second-order Ordinary Differential Equation we use the ode45 command where we called function equation from ode_fun from function code as we mentioned earlier.
• With the help of this command, we find the values of [t, result] which we can use to plot the movement of the pendulum with respect to its angular velocity and angular displacement.
• To plot the plot of angular velocity and angular displacement we use the plot command first we plot angular displacement using the plot command plot(t,result(:,1)) then we used hold on command to stop overwriting of this curve over another curve.
• And to plot angular velocity we used plot(t,result(:,2)) then we use xlabel and ylabel commands to give names to the x and y-axis as time and plot respectively.
• Then title command is used to name the plotted graph as pendulum displacement.
• To make an animation of the pendulum we used counter command and a for loop command in for loop I variable is used to 1:length(t) to get the time values which can be used to find the values of x and y variable.
• Then to plot a pendulum we used plot command as plot([-1,1],[0,0],'linewidth',5) with this we are able to plot the horizontal bar of a pendulum, and for the vertical bar we used plot([0,x],[0,y],'linewidth',5) within this command we used hold on command to stop overwriting the plot command hold off command is used after that.
• plot(x,y,'Marker','o','Markersize',25,'Markerfacecolor','r') is used to plot the circular shape of a pendulum with all this we successfully plot our pendulum with animation.
• To give the pendulum a fixed motion we used the axis command at a point axis([-1.5,1.5,-2,2]);
• title command is used to name the figure as a pendulum.
• M(ct) command is used after a pause command to capture the current frame every time for 'for loop'.
• The code ct=ct+1 is used for incrementing the frame every time the for loop runs.
• After that, we end the for loop using the end command.
• And in the end, to make a video of the plotting of a pendulum we used the movie command after that video file command is used where 'pendulum_movement' is the name of the video which is about to save in the device and 'Uncompressed AVI' is the file type.
• The open command is used to open the video and the WriteVideo command is used to write plots in a video format after that close command is used to close the video.
• After writing this we used the run button which is in the taskbar or use f5 to run the program.
• With the help of the above codes, we get our pendulum in video formal as an output successfully.

Result:-

In the above graph with the help of MATLAB, we plot a graph to simulate the motion of a pendulum between 0-20 sec, for angular displacement at 0, the angular velocity at 3 rad/sec at time t=0.

The animation of the movement of the pendulum is in video format which can be accessed through this link https://youtu.be/Z9dlKzdCFuA

Conclusion:-

With the help of MATLAB, we solve the second-order Ordinary Differential Equation and we created an animation file of the movement of a pendulum successfully.