AIM

To write code in MATLAB that can solve second order ODE and animate the motion of simple pendulum.

INTRODUCTION

1.       ORDINARY DIFFERENTIAL EQUATIONS:

In mathematics, an ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions. The term ordinary is used in contrast with the term partial differential equation which may be with respect to more than one independent variable. [1]

An ordinary differential equation (frequently called an "ODE," "diff eq," or "diffy Q" is an equality involving a function and its derivatives. An ODE of order n is an equation of the form:

A simple example is Newton's second law of motion — the relationship between the displacement x and the time t of an object under the force F, is given by the differential equation: [1]

2.     THEORY OF PENDULUM:

A pendulum is a body suspended from a fixed support so that it swings freely back and forth under the influence of gravity. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging it back and forth. The mathematics of pendulums are in general quite complicated. Simplifying assumptions can be made, which in the case of a simple pendulum allow the equations of motion to be solved analytically for small-angle oscillations. [3]

            

Figure 1. Description of FBD of pendulum.

 

3.      GOVERNING EQUATIONS:

A simple pendulum consists of a ball (point-mass) m hanging from a (massless) string of length L and fixed at a pivot point P. When displaced to an initial angle and released, the pendulum will swing back and forth with periodic motion. By applying Newton's second law for rotational systems, the equation of motion for the pendulum may be obtained:

and rearranged as: [4]

 

This is analogous to the equation of equilibrium of vibration where the damping is neglected. But, in the real scenarios, the damping is a particularly important aspect that cannot be neglected. Therefore, introducing a damping coefficient ‘b’. the equation is then modified as follows:

 

The equation (4) is a second order ordinary differential equation. This is further broken into series of first order ordinary differential equations as follows:

Substituting equations (5), (6) and (7) in the equation (4), we get:

In Matrix format, the above equation is rewritten as:

 

OBJECTIVES

• To write a code in MATLAB that solves ODE pertaining to simple pendulum.
• To animate the damped swing of the simple pendulum.

 

 

PROGRAM – CODE WRITTEN IN MATLAB

PROBLEM STATEMENT:

A simple pendulum obeys the ordinary differential equation (4). It is required to generate a plot on how the angular displacement and angular velocity of the simple pendulum vary with time. The animation of the pendulum swing is to be written in MATLAB.

Technical data of pendulum is as follows:

Table 1. Pendulum parameters

CODE SIMULATION:

clear all

close all

clc

 

b = 0.05;%damping coefficient

g = 9.81;%acceleration due to gravity

l = 1;%length of the pendulum

m = 1;%mass of the pendulum

 

%initial conditions

theta_0 = [0;3];%array indicating the initial value of diaplacement is 0 and initial velocity is 5m/s

 

%timeplot

t_span = linspace(0,20,1000);%time begins at 0s to 10s with 500values in between

 

%calling ODE to solve ODE

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

 

%plotting displacement versus velocity curve

figure(1);

plot(t,results(:,1),'linewidth',1,'color','[0.2 0 0]');%plotting angular displacement in brown

hold on;

plot(t,results(:,2),'linewidth',1,'color','[1 0.5 0]');%plotting angular velocity in orange

xlabel('Angular Displacement (\theta) and Angular Velocity (\omega)');

ylabel('Time (s)');

legend('Angular Displacement (\theta)','Angular Velocity (\omega)');

title('Variation of Angular Displacement and Angular Velocity wrt Time');

 

 

ct =1;

for i = 1:length(results)

  

    angle = results(i);

   

    x0 = 0;

    y0 = 0;

   

    x1 = l*sind(angle);

    y1 = -l*cosd(angle);

   

    figure(2);

   

    plot([x0 x1],[y0 y1],'linewidth',2,'color','k');

    hold on;

    plot([-0.03 0.03],[0 0],'linewidth',10,'color','k');

    plot(x1,y1,'marker','o','markersize',30,'markerfacecolor','g','markeredgecolor','k');%markerfacecolor

    axis([-0.07 0.07 -1.5 0.5]);%definition of the lenght of the axes

    time = results(i,1)/results(i,2);

    txt = sprintf('Time = %f seconds',t(i));

    text(-0.02,0.25,txt);

    title('HARMONIC MOTION OF THE PENDULUM');

    pause(0.03);%used to pause the code

    M(ct) = getframe(gcf);

    ct = ct +1;

    hold off;

    end

 

 

movie(M) %movie is created with all the frames in M

videofile = VideoWriter('simple_pendulum_assignment.avi','Uncompressed AVI') %contains all the frames

open(videofile)%command to open video file

writeVideo(videofile,M)%to write every captured frame on videofile

close(videofile)

 

 

EXPLANATION:

• The code begins with the introduction of the technical data required to describe the pendulum, as tabulated in the Table.1.
• An array ‘theta_0’ is introduced that holds the values of initial angular displacement and velocity at time t=0.
• ‘time_span’ is introduced using linspace with time step from 0 second to 20 seconds with 1000 intervals in between.
• A new command named ‘ode45’ is used to solve the differential equation as per the problem statement.

ODE45 FUNCTION: Ode45 is a function in MATLAB that solves ordinary or partial differential equations using a numerical method, Runge Kutta method. It is used to solve non-stiff differential equations — medium order method.

Syntax:

[t,y] = ode45(odefun,tspan,y0)

 where tspan = [t0 tf].

integrates the system of differential equations y′=f(t,y) from t0 to tf with initial conditions y0. Each row in the solution array y corresponds to a value returned in column vector t.

All MATLAB® ODE solvers can solve systems of equations of the form y′=f(t,y), or problems that involve a mass matrix, M(t,y)y′=f(t,y). The solvers all use similar syntaxes. The ode23s solver only can solve problems with a mass matrix if the mass matrix is constant. ode15s and ode23t can solve problems with a mass matrix that is singular, known as differential-algebraic equations (DAEs). Specify the mass matrix using the Mass option of odeset.

ode45 is a versatile ODE solver and is the first solver you should try for most problems. However, if the problem is stiff or requires high accuracy, then there are other ODE solvers that might be better suited to the problem. [5]

• The following code is explained as follows:

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

In this, @(t,theta) indicates defining the independent variable ‘t’ and dependent variable ‘theta’ which are the variables of the ordinary differential equation.

The function ‘ode_fun()’ is a function that is created to solve the equation (9) that points to the differential equation in problem statement but in matrix format. It is used to define input parameters.

The function is described as follows:

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


theta1 = theta(1); %angular displacement

theta2 = theta(2);

dtheta1_dt = theta2; %angular velocity

dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1); %angular acceleration

dtheta_dt = [dtheta1_dt ; dtheta2_dt];

 

end

 The inputs to this ode_func are the values of time span, theta, damping coefficient (b), acceleration due to gravity (g), length of the pendulum (l), mass of the bob (m).

‘t_span’ differentiates the differential equation in the specified time intervals i.e. determines how far can the operation be performed. ‘theta_0’ indicates the initial condition that the effect is performed on set of parameters.

• A new figure frame i.e. figure(1) is created to plot the variation of angular displacement and angular velocity with respect to time.
• A legend is created to identify different parameters clearly. The plot representing the variation of angular displacement and angular velocity is plotted with respect to time span.
• A counter variable ct is introduced to formulate the animation required in for loop.
• The limits of the for loop begins with 1 and stretches along the length of the results.
• A variable ‘angle = results(i)’ is introduced to define the values of the path that the pendulum takes while swinging.
• Another figure command is used to capture the motion of the pendulum. i.e. figure(2).
• The initial coordinates for the pendulum’s cord is (0,0) and extends to (x1,y1) and it is colored in black. It is created using ‘plot’ command.
• The values of (x1,y1) encompass the values that lie in the third and fourth quadrant of the cartesian plane.
• A horizontal fixed end is plotted with a line parallel to X-axis in black color with linewidth of 10 using plot command.
• The pendulum bob is created using marker ‘o’ which is filled with green using ‘markerfacecolor’ and edged with black using ‘markeredgecolor’. The size of the marker is maintained at 30 using ‘markersize’.
• Axis command is used to define the length of the axes as:

X-axis :-0.07  to 0.07          Y-axis: -1.5 to 0.1

• A timer is placed in the animation that counts 20seconds which stands as a clear standard to identify the swing of the pendulum. The following code is used to set up a timer:

time = results(i,1)/results(i,2);

txt = sprintf('Time = %f seconds',t(i));

text(-0.02,0.25,txt);

time is evaluated as the ratio od angular displacement to angular velocity that are held by results(1) and results(2) respectively.

Then a variable ‘txt’ is introduced to mention the time as float (data type) as there is an increment on ‘i’.

Sprint command is used to format data into string or character vector.

The location of the text “Timer = …. Seconds” is at (-0.02,0.25) on the cartesian plane described earlier while declaring the axes and is done using the command ‘text’

• The plot is titled ‘'HARMONIC MOTION OF THE PENDULUM’ using ‘title’ command. Pause (0.03) identifies that the code pauses for 0.03 seconds for every change in the value of angle.
• In this program, an array named ‘M’ is introduced to store every frame regarding the orientation of the pendulum.
• Getframe(gcf) is assigned to the array M which makes sure that each figured is captured as a movie frame.
• ct=ct+1 is used to make sure that the results are not repetitive, and the loop advances by 1 unit.
• Then, ‘movie’ command is used on the array M that stores all the movie frames thereby creating a movie.
• This movie is stored under the name ‘simple_pendulum_assignment.avi' as an uncompressed AVI with the help of ‘videowriter’ command.
• This data is stored under the name ‘videofile’.
• Then a command named ‘open()’ is used to open the videofile.
• Following that, a command named writeVideo() is used to write every captured frame from the array M onto videofile.
• ‘close’ command is used on the videofile to close it.

 

OBSERVATION AND CONCLUSION:

• While creating the animation, in for loop:

for i = 1:length(results)

angle = t_span(i);

The value assigned to ‘angle’ was t_span intilly. This resulted in tracing the swing of the pendulum on in one direction as shown in the figure (2). Hold off command was withdrawn to understand the path clearly.

Figure 2.  Error(1)

• The plot depicting the 'Variation of Angular Displacement and Angular Velocity with respect to Time' is shown below:

Figure 3. Variation of Angular Displacement and Angular Velocity with respect to Time (t=20 seconds)

 

Angular displacement is shown in brown colour whilst angular velocity is shown in orange. The initial value of the angular displacement on the graph is ‘0’ and that of angular velocity is ‘3’ as mentioned by “theta_0” in the code. The effect of damping is clearly observed from the above figure. The values of the displacement and velocity decrease by the end of time interval. If the time-period is further increased the damping is clearly observed.

Figure 4. Variation of Angular Displacement and Angular Velocity with respect to Time (t=90 seconds)

Figure 5. Variation of Angular Displacement and Angular Velocity with respect to Time (t=180 seconds)

•  The different values of the parameters can be observed in the workspace as shown in the figure 6. The ‘results’ constitute two columns that are shown in a tabular form as shown in the figure 7.

       

Figure 6. Workspace in MATLAB

         

Figure 7. Tabulation of ‘results’

• After running the code, the command window displays the characteristic features of the video. It gives the information about the general properties that include the location of the video in the system and, the type of file in which the video is saved. The video properties specify the color channels, frame rates, video bits per pixel, video formats and so on. This information is shown in the figure 8.

Figure 8. Properties displayed in Command Window

RESULT

The simulated swing of the pendulum in the form of a movie is uploaded in the link given below:

https://youtu.be/gCnucodNqiE

CONCLUSION

By modifying the above code, the energy contours can be identified. The complex problems of compound pendula can also be simulated. Furthermore, to challenge the results, different types of numerical methods can be employed to identify the accuracy of the solution.

BIBLIOGRAPHY

1. https://en.wikipedia.org/wiki/Ordinary_differential_equation
2. https://mathworld.wolfram.com/OrdinaryDifferentialEquation.html
3. https://en.wikipedia.org/wiki/Pendulum_(mathematics)
4. https://www.acs.psu.edu/drussell/Demos/Pendulum/Pendulum.html#:~:text=By%20applying%20Newton's%20secont%20law,enough%2C%20so%20the%20small%20angle
5. https://in.mathworks.com/help/matlab/ref/ode45.html