Solving second order ODEs

AIM: To write a program in MATAB to simulate the motion of a pendulum.

THEORY: 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. ^[1]

Assumptions

Various assumptions taken for the program are:

• Friction from both air resistance and the system is negligible.
• The pendulum swings in a perfect plane.
• The arm of the pendulum cannot bend or stretch/compress.
• The arm is massless.
• Gravity is a constant 9.8 m/s²

Derivation of the Equation

Let,

m = Mass at the swinging end of the pendulum (kg)

g = Acceleration due to gravity (m/s²)

L = Length from the swivel point to the centre of mass (m)

θ = Angle between the string position to the string position at rest (rad)

t = Time (s)

T = Time period of the pendulum (s)

From Newton’s second Law,

F = ma

Acceleration due to gravity will be a function of θ. At,

θ =π/2 → |a| = g

θ = 0 −→ a = 0

Considering the relation of acceleration and θ we arrive at

a = −g sin θ (1)

Arc length (arc_L) of the pendulum can be thought of as the “position” of the system

arc_L = Lθ

Acceleration of the system will now be

Plug in equation (1) to get the simple harmonic motion of a pendulum shown in equation (2)

Now we will solve equation (2) to get T (period) reduce the second order differential equation to a first order

Integrate to get the differential equation

With the initial conditions of θ’(0) = 0 and θ(0) = θ₀ we can solve for C

Put our C into the equation

Now take the square root of both sides while ignoring the negative because we are solving for time and either time will be the same distance from our t0. It is just a matter of forwards or backwards in time. On the left hand of our equation lies the rate of change of the angle with respect to time, but we are going to solve for the period, so we need the time with respect to the angle, because of this we are going to inverse the entire equation and integrate from 0 to θ0. We will now multiply the whole thing by four to get the period. The change in time to get from 0 to θ0 is only one fourth of the entire cycle of the pendulum. This gives us our new equation of

This is where MATLAB can prove to be a great asset because equation (3) cannot be solved in terms of elementary functions.

PROGRAME:

1) Main Program

clear all 
close all
clc

% Value of the variables
b=0.05;      
g=9.81;     
l=1;         
m=1;         

% Initial conditions
theta_0=[0;3];

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

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

% Loop fo calculation of position of the Pendulum
ct=1;
for i=1:length(results)
    x1=l*sin(results(i,1));
    y1=-l*cos(results(i,1));
    plot([-0.75,0.75],[0,0],'linewidth',3);
    line([0,x1],[0,y1],'linewidth',3,'color','black')
    axis([-1.5 1.5 -1.5 0.5]);
    grid on;
    hold on;
    plot(x1,y1,'o','markers',15,'markerfacecolor','r');
    pause(0.01);
    M(ct)=getframe(gcf);
    hold off;
    ct=ct+1;
end

% Generation of animation
movie(M)
videofile=VideoWriter('pendulum_animation.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

2) User-defined function (ode_func)

% A user-definaed function for the calculation of values of displacement
% and velocity with respect to time.
function [dtheta_dt]=ode_func(t,theta,b,g,l,m)

theta1=theta(1);
theta2=theta(2);

% Displacement
dtheta_dt1=theta2;
% Velocity
dtheta_dt2=-(b/m)*theta2-(g/l)*sin(theta1);
[dtheta_dt]=[dtheta_dt1;dtheta_dt2];
% [dtheta_dt] returns the final calculated value of which first column is
% displacement and second column is velocity.
end

OUTPUT:

1) Output of Workspace

2) Output of animation

The output of animation is uploaded on YouTube with link:

https://youtu.be/n5lEvdBLinM

This animation was generated in avi format, but converted to mp4 format using FormatFactory software to reduce the size of the file.

EXPLANATION:

• clear all, close all, clc are written for clearing the workspace, command window and also to close any plots.
• Input values are taken for parameters like length of string as ‘l’, mass of bob as ‘m’, damping coefficient as ‘b’, and initial condition as ‘theta_0’ which is equal to [0,3] which denotes displacement as 0 and velocity as 3 m/s².
• Time interval is taken as ‘t_span’. linspace() is used here to generate 500 values between 0 and 20.
• Every iteration of the loop will be saved as a frame for generation of the animation, hence for this a counter is used as ‘ct’, which will change its value for every iteration.
• Loops will be used for calculating the position of bob. Loop runs for 500 times since ‘length(results)’ gives 500 values. For every iteration, position of bob is calculated and then plotted on the graph using plot(). axis() is used for giving the values to the axis of the graph. pause() is used for giving pause to every iteration of the loops, so that frames can be saved for animation. ‘hold on’ command is used to hold all three plots together. The first plot is for horizontal line, which act as the point where the string of bob is mounted. Another plot shows the position of the bob. Last plot show the bob itself in red colour. ‘hold off’ command is used to deactivate ‘hold on ‘ command. In this manner position of bob with respect to time is calculated.
• Animation is created by generating frames for each plots. getframe() is used for generating frames for each iteration. Frames are saved in ‘M’. ‘M’ is an array of 1 row. Its column value is decided by the no of iteration which is 500. Each frame is saved in new column of ‘M’, hence counter value (ct) needs to be increased every time. That why ct+1 is used, which increases its value by 1 after every iteration.
• After the for loop is completed, the generated frames are compiled together to make an animation and saved in ‘videofile’. First ‘videofile’ is created using VideoWrite(), then the file is opened using open(), then all frames are compiled together and written on it using write(), then the file is closed using close().
• The user-defined function ‘ode_func’ is used to calculate the final values of the displacement and velocity with respect to time.

Errors

1) Error for misspelling ‘VideoWriter’ command

This error was corrected by changing the name from ‘theta’ to ‘VideoWriter’.

CONCLUSION

This program helps in understanding the position of bob in simple pendulum.

Animation can be generated smoother by providing higher number of loops.