SOLVING SECOND ORDER ODE USING MATLAB

1.AIM :

WRITE A PROGRAM TO SOLVE THE ODE REPRESENTING THE MOTION OF SIMPLE PENDULUM AND ANIMATE THE MOTION

2.THEORY AND GOVERNING EQUATIONS :

A simple pendulum consists of a bob of mass m hanging from a string of length l and fixed at a pivot point. When displaced to an initial angle and released, the pendulum will swing back and forth in simple harmonic motion. By applying Newton's secont law for rotational systems, the ODE representing the equation of motion of a simple pendulum with damping is given by

`(d^ 2 θ)/ (d t^ 2) + b/ m ⋅ (d θ)/ (d t) + g/ l ⋅ sin θ = 0`

where

g = gravity in m/s^2

l = length of the pendulum in m

m = mass of the ball in kg

b=damping coefficient

3.OBJECTIVES :

1. To solve the second order ODE.
2. To plot the angular velocity v/s time and angular displacement v/s time graph.
3. To animate the motion of the simple pendulum.

4.MATLAB PROGRAM :

Program with explanation in form of comment is given below

clear all
close all
clc

%Taking inpyt from user

%damping coefficient
b=input('Enter the value of damping coefficient : ');

%length of pendulum
l=input('Enter the length of pendulum : ');

%mass
m=input('Enter the mass : ');

%gravity
g=9.81;

%initial displacement
i_disp=input('Enter the value of initial displacement : ');

%initial velocity
i_vel=input('Enter the value of initial velocity : ');

%initial condition
theta_0=[i_disp ; i_vel]; % making a column matrix of initial displacement and initial velocity

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

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

%plotting the displacement and velocity
figure(1)
plot(t,results(:,1)) %all results of column 1 which is displacement
hold on
plot(t,results(:,2)) %all result of coulmn 2 which is velocity
ylabel('Plot')
xlabel('time')
legend('Angular displacement','Angular velocity');

%plotting for the animation of simple pendulum

%initializing counter variable ct
ct=1;

for i= 1:length(results(:,1))
    
    % co-ordinates of string attached to hook/support
    x0= 0;
    y0= 0;
    theta_1 = results(i,1);
    
    % co-ordinates of string attached to bob
    x1 = l*sin(theta_1);
    y1 = -l*cos(theta_1);
    
    % plotting
    figure(2)
    
    % plotting base of pendulum
    plot([-1,1],[0,0], 'linewidth',6,'color','k')
    hold on
    % string of pendulum
    line([x0 x1],[y0 y1],'linewidth',6,'color','b')
    hold on
    % bob of pendulum
    plot(x1,y1,'marker','o','markersize', 28,'markerfacecolor','g')
    axis([-2 2 -2 2]);
    pause(0.03)
    
    % storing the current plot figure in an array 'A' at position of value of c
    A(ct)=getframe(gcf);
    
    hold off %to overwrite the next plot in the next iteration
    
    %increasing the value of counter by 1
    ct=ct+1;
end
% creating movie from captured frames

%command to generate the movie of array 'A'
 movie(A);
 
 %creating a video named simple_pendulum of type '.avi'
 videofile = VideoWriter('simple_pendulum.avi','uncompressed AVI');
 
 %opening the videofile
 open(videofile);
 
 %write the movie generated using array 'A' to the videofile
 writeVideo(videofile,A);
 
 %closing the videofile
 close(videofile);

The code for function ode_func is

function [dtheta_dt]= ode_func(t,theta,b,g,l,m)

    theta1= theta(1); %assigning the intial displacement condition
    theta2= theta(2); %assigning the intial velocity condition
    dtheta1_dt= theta2; %different of theta1 (displacement) is velocity
    dtheta2_dt=-(b/m)*theta2 - (g/l)*sin(theta1);
    dtheta_dt= [dtheta1_dt; dtheta2_dt];

end

EXPLANATION:-

1. Initially, the input for damping coefficient, length, mass, initial displacement and velocity is taken from user.
2. The timespan over which the moiton is to be studied is set.
3. We solve the ode with help of inbuilt 'ode45' function and user defined function (here 'ode_func').
4. The syntax of function ode45 can be explained by breaking the above used function ( ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0); ) in parts. The '@(t,theta)' represents that the ODE is function of these 2 variables, the part 'ode_func(t,theta,b,g,l,m)' is calling the user defined funtion 'ode_func' which solves the second order ODE for simple pendulum and the part 't_span,theta_0);' indicates the assigning of initial condition to the 2 variables on which the function id dependent.
5. The results obtained by solving ode is stored in an array (here 'results') and the angular velocity v/s time and angular displacement v/s time graph is plotted.
6. The plots are converted into an animation using the 'movie()' command.

RESULTS:-

Output window screenshot:-

plot:-

The animation of the motion of simple pendulum can be viewed using the below link:-

https://youtu.be/oj5pagS4qmg

Errors faced:-

Not using the 'hold off' command led to printing all the plot of pendulum in same window

5.CONCLUSION :

1. The second order ODE was solved.
2. The graph for angular velocity v/s time and angular displacement v/s was time obtained.
3. The animation for the movement of simple pendulum was obtained.