Solving second order ODEs in MATLAB

SIMPLE PENDULUM MOTION WITH DAMPING

AIM :

To write a program to solve ordinary differential equation (ode) for equation of motion of a simple pendulum with damping and to create an animation for the motion of the pendulum.

THEORY :

• The simple pendulum is formed of a light, stiff, inextensible rod of length l with a bob of mass m. Its position with respect to time t can be described merely by the angle θ (measured against a reference line, usually taken as the vertical line straight down).
• The fact that its position can be described by using only one variable means that the simple pendulum has only one degree of freedom.
• An example of a damped simple harmonic motion is a simple pendulum.

GOVERNING EQUATIONS :

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation(ode) that describes the position of the "ball" w.r.t time.

The 2^nd order ode below represents the equation of motion of a simple pendulum with damping -

• `(d^2θ)/dt^2+b/m (dθ)/dt+g/Lsinθ=0`

where, b=damping coefficient.

m = mass of the ball in kg,

g = gravity in `m/s^2`,

L = length of the pendulum in m.

In order to solve the equation of motion for pendulum and as the ode is of higher order we break the 2^nd order ode into a series of 1^st order ode.

The final ode equations to solve is represented as :

• `(dθ1)/dt=θ2` _ _ _ _ _ _ _ _ _ _ _ _ _ _(equation-1)

• `(dθ2)/dt=-b/mθ2-g/Lsinθ1`_ _ _ _ _ _ _(equation-2)

OBJECTIVES OF THE PROJECT:

1. To create a program that solves the ODE for the equation of motion of a simple pendulum with damping and to describe the transient behavior of the system.
2. Get the position of the pendulum as a function of time and use that data to animate the motion of the pendulum and create a movie of the animation respectively.

• FUNCTION CODE -

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

• MAIN CODE -

%% Code for solving simple pendulum equation using ode45
clear all
close all
clc

% Inputs -
b = 0.05;     % Damping
m = 0.1;      % Mass in kg
g = 9.81 ;    % Gravity
l = 1 ;       % Length

% Initial Condition 
theta_0 = [0 ; 5];         % Intializing the condition as 0 for displacement and 5 for velocity

% Time interval/ Time points
t_span =  linspace(0,10,400); 

% Solve ode equation using ode45 equation solver
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m), t_span, theta_0);   %Calling the function saved ode_func to solve ode for pendulum

% Plot the displacement and velocity with respect to time
figure(1)
plot(t,results(:,1))                 % Plotting time vs all rows in 1st column
hold on
plot(t,results(:,2))                 % Plotting time vs all rows in 2nd column
xlabel('Time(t)')
ylabel('Plot of Displacement and Velocity')
legend('DISPLACEMENT','VELOCITY')

% Plot for different positions of pendulum for creating an Animation 
counter_ct = 1;
THETA_displace = results(:,1);
for i = 1:length(THETA_displace)
    THETA_0 = THETA_displace(i);
    
    % Assigning the co-ordinates -
    x0 = 0;
    y0 = 0;
    
    x1 =  l*sin(THETA_0);               % Get the equation from the geometery of the pendulum
    y1 = -l*cos(THETA_0);
    
    % Plotting the pendulum for motion
    figure(2)
    plot([-1.2 1.2],[0 0], 'linewidth',4,'color','b')             % Plot for Fixed line 
    hold on
    plot([0 0],[0 0],'o','markerfacecolor','k','markersize',8);   % Plot for hinge point in fixed line
    hold on
    plot([x0 x1],[y0 y1], 'linewidth',3.5,'color','g')            % Plot for pendulum string as line with green color
    hold on
    plot(x1,y1,'o','markersize',12,'markerfacecolor','r')         % Plot for mass bob fixed to string
    hold off
    axis([-2 2 -2 1])                                             % Setting the axis limits in plot
    pause(0.01)   
    
    E(counter_ct) = getframe(gcf);                                % Capturing figure and plot in frames and saving in Variable E
    counter_ct = counter_ct + 1 ;
end

% Creating Movie and saving it in file
movie(E)                                                    % Creating movie for all the saved frames of plots in Variable E
videofile = VideoWriter('Simple Pendulum.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,E)                                     % Writing the video file to movie frames saved to Variable E 
close(videofile)

• OUTPUTS :

• The above output shows the plot for displacement and velocity vs time for damping i.e the results with the time period defined in the code respectively.

• The above output shows the general data and the properties for the saved video animation created as movie in a ' AVI ' format.

STEPWISE EXPLANATION OF THE CODE :

1. We start the program by clearing the window and previous plots by using basic commands clear all ; close all ; clc .
2. Initially, after that we give the input parameters for the simple pendulum’s equation of motion with damping like length of pendulum, damping coefficient, gravity and the mass of ball.
3. Then set the initial conditions for displacement and velocity of pendulum. Here, we initialize the displacement as 0 and velocity as 5 m/s^2.
4. Then next the time interval or time span is defined for pendulum using linspace command for time array from 0 to 10 with 400 points in between.
5. The ode for pendulum is solved using the inbuilt command ode45 which is used for solving ordinary differential equation in matlab. The syntax is represented as, [t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m), t_span, theta_0)
6. Here, we start with the command ode45 then, give the main dependent and independent parameters for ode equation for pendulum i.e. t and theta next call the function defined as ode_func and giving all its saved input parameters like t, theta , b , g , l , & m.
7. Then mention the t_span as defined before and give the initial condition theta_0 for solving ode respectively.
8. Plot command has been used to plot the rate of change of displacement with respect to time and velocity with respect to time where results shows the 2 columns 1^st for displacement and other for velocity.
9. The counter is used and the for loop for creating a plots or figures for animation of pendulum motion.
10. In for loop we had defined the co-ordinates x0,y0 initial which is fixed and x1,y1 for pendulum mass attached for motion using the formula concluded from the geometry of the figure.
11. Then plotted the fixed line for pendulum where the string is attached and then plot the mass attached to string using command plot command with hold on for all the plots in a figure and also set the axis limits.
12. Then saved the frames for movie for different figures using getframe command and gcf for the current figure in a variable E.
13. Then increment the counter by 1 and ended the for loop.
14. Finally, created the movie using movie command for the saved frames in variable E respectively.
15. Use the video file to write the file name with the format and save it.
16. Then, open the file, write the file for Movie E and closed it using command close(videofile) respectively.
17. Finally, after running the program / code we get the respective plot for time and results and the animation showing simple pendulum motion.

Errors faced while programming :

1. Misplaced position for mass in plot due to syntax error in plot command which is use of square bracket for plotting 2 variable x1 y1 for bob.

i.e. plot([x1 y1],'o','markersize',12,'markerfacecolor','r')

b. Error in creating a movie due to typo for command movie, so matlab cannot find the instruction to run further.

Screenshots showing errors in command window and in Plot :

• Error for plotting mass of bob

• Error for creating a movie

Steps taken to overcome the Errors Faced –

• Corrected the syntax for plot and eliminating the square brackets for variables x1 & y1 for creating a mass or bob attached to string.

i.e. plot(x1, y1,'o','markersize',12,'markerfacecolor','r')

• Use the right syntax command movie instead of using capital M as Movie which is wrong.

CONCLUSION :

Analyse and understood the equation of motion with damping of a simple pendulum and created a program for solving ode for pendulum using ode45 equation solver for the transient behaviour of the system and animated the motion of the pendulum when the displacement changes with respect to time. In final output graph for displacement and velocity vs time we learned that as time increases then the displacement of the pendulum decreases so thus, the velocity for the motion of pendulum respectively.

Reference link for the animation

https://youtu.be/_JdOhDVF_Pk