Soving second order ODE to simulate a simple pendulum

Aim: Differential equations are the most important concept in mathematics that links this beautiful subject to all the other subjects. A differential equation act as a medium through which real physical systems can be theorized in a mathematical model for its analysis. Literally any physical system in the world can be described as a set of differential equation(s) in the form of a mathematical model. Here Ordinary Differential Equations (ODE) are focussed on. A simple pendulum is a physical system and thus its motion is governed by a second-order ODE. In this work, the solution of this equation is performed and thus visualizing the motion of the pendulum.

Governing Equations: The motion of a simple pendulum can nearly be approximated as the solution of the following differential equation:

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

where b is the damping coefficient;

m is mass (in kg);

g is the acceleration due to gravity (in m/s^2);

L is the length of string/rod of the pendulum (in m);

`theta` defines the angular position of the pendulum from the central imaginary axis

Given Parameters:

L= 1 m; m=1 kg; b=0.05; g=9.81 m/s^2

Objective: As clear from the sections above, the work is focussed on trying to solve the Governing differential equation to study the motion of the simple pendulum. But one thing to note here is that the major scope of this project is to find the solution of the second-order differential equation. And the study of motion is secondary. MATLAB script is prepared to solve the equation and ode45 tool is used.

Matlab Code:

Firstly, defining the function in the Matlab that would later be used in the main code.

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

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

dtheta1_dt=theta2; % Self assumption to convert second order ODE to first order ODE
dtheta2_dt= -(b/m)*theta2-(g/l)*sin(theta1); % This the equation we obtain after putting the former assumption in the governing equation

dtheta_dt=[dtheta1_dt;dtheta2_dt]; % Define as a column vector
end

Now the main body code:

% Program to solve the ODE and procuring the animation of the Pendulum system

% Initializing Constants of the problem
b=0.05; % Damping constant
g=9.81; % Acceleration due to gravity (in m/s^2)
l=1; % Length of pendulum string (in m)
m=1; % Mass of bob (in kg)

% Initial Boundary Conditions- First row in the following vector represents
% initial Displacement and the second row represents initial velocity
theta_0=[0;3];

% Time points
t_span=linspace(0,20,500); % Creating 500 divisions in the 0 to 20 sec time span for more smooth simulation

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

% Plotting
ct=1;
for i=1:length(results(:,1))
plot([0 l*sin(results(i,1))],[1 l*(1-cos(results(i,1)))],'linewidth',3,'color','r') % This plots the string
axis([-1 1 -0.2 1.2]) % Setting suitable limits for axis to be visible
ylabel('Y-axis')
xlabel('X-axis')
hold on
plot(l*sin(results(i,1)),l*(1-cos(results(i,1))),'marker','o','MarkerSize',30,'MarkerFaceColor','b') % This plots the Bob at the end of the string
hold on
plot([0 0],[-0.2 1.2],'--k') % This plots the axis representation
hold on
plot([-0.4 0.4],[1 1],'linewidth',5,'color','g') % Plots the Wall representation in the figure

hold off % This erases off the data of the current iteration and prepares the screen for the next iteration

pause(0.01)

M(ct)=getframe(gcf); % Putting all the frames or individual figures in each iteration of the loop in a row vector, M
ct=ct+1;
end

movie(M)
video_file=VideoWriter('Animation_Pendulum.avi');
open(video_file)
writeVideo(video_file,M)
close(video_file)

Step-wise Explanation of the codes above:

First, explaining the main code. Function definition code would be explained later when it would be called in the main code.

First, declaring the given constants in the problem statement. These are also explained in the section 'Given Parameters'. Define a new variable as "theta_0". This variable is defined as a column vector, storing the initial conditions of the given system. Initial conditions of the system are: Angular displacement, 0 and angular velocity, 3 rad/s. Then define the given time span of the problem, that is, zero to 20 seconds. This is the time for which, the problem is to be studied/ analyzed. So every calculation in this project is applicable to this time span only. For better visualization of its motion, increase the number of divisions in this time span. To solve the defined ode (ODE is defined in the function code above), MATLAB has a built-in function known as ode45(). MATLAB has many other type of solvers. But most common among them is ode45. The working and the technique of its calculation of this function is out of the scope of this project. Now, ode45 needs an equation to solve, this equation is defined in the function code above.

So now coming to the function code:

The name of the defined function is ode_func, which gives the output as the differential, `(d(theta_2))/dt`, where `theta_2=(d(theta))/dt` and `theta=theta_1`. We are basically reducing the order of the equation to one. So that ode45 can solve it. `theta_2` is nothing but a self-assumed variable to reach the solution. This function takes input as our initially defined constants and t, `theta`, as the equation is in these two parameters. These two are given as inputs in the function, as it is necessary so as to use ode45(). After substituting the `theta_2` in the main equation, we obtain the output of the function, as a first-order ODE. Define "theta" as a column vector which stores `theta_1` and `theta_2`. After the definition of "theta", "dtheta1_dt" and "dtheta2_dt" are defined as per above explanation. After substituting `theta_2` in the main equation, `(d(theta_2))/dt=-b/m*theta_2-g/L*sin(theta)`. The output of this function is nothing but a column vector containing the differential of the vector "theta".

Now, we pass this function in ode45() along with the initial conditions defined in the problem. The output of ode45 is stored in two variables, "t" and "results". "t" indicates the time with the divisions created earlier, and "results" is 500 by 2 matrix, storing the output of angular displacement and angular velocity.

Now for plotting:

Since, velocity and displacement are synced with each other, plotting on the basis of just one can solve our purpose. Since we need a continuous simulation video as our output and not just one plot for 1 position, a loop is required, which would plot every position defined by the total number of divisions between the time span in each iteration. And simultaneously clearing the previous plot as it moves to the next iteration. So first the plot for string/rod of the pendulum is done in terms of iteration variable, which is a straight line. axis() command controls the visible range of axis in the plot and is thus, set suitably. To visualize the bob of the pendulum, it is necessary to plot that particular point separately and mark the point with colour filled bold marker. And with each subsequent plots, "hold on" command is essential as it keeps the plot intact unless "hold off" command is given. Now next plot command draws the imaginary central axis of the pendulum (mean position). And the last plot command draws the wall/boundary from which the pendulum is suspended. The last two plots are just for aesthetics. "linewidth" and "color" are given accordingly. Now "hold off" command is necessary. "pause()" is for the time specification to rest on each iteration of the loop. Now, simulation is ready! But the job is not done yet. To take this simulation out of MATLAB in the form of some animation, movie() function is used. So movie basically combines all the frames obtained from the loop iterations. This collection of the frames is done by "M" variable. It is nothing but a vector which retrieves the frames by "getframe(gcf)" function. And "ct" is just the iterating variable for "M".

Movie making:

As described above, movie function joins all the above-collected frames to make a video file. "VideoWriter()" function creates a file with the specified name inside the function. Then open() the file. "writeVideo()" appends this file with all the frames in "M". And lastly, close the file.

The obtained video from this code can be viewed from the link below:

https://www.youtube.com/watch?v=zR2u9rvfgyA

Errors faced:

I faced two errors and one thing was common in both of them, that I didn't use "hold off" command. So all the iterations kept on accumulating on the same plot figure.

First one's partial code was something like this (text form of code won't be given here as this caused error so it would not be appropriate to give the code in text format):

Here, it can be seen that the y coordinates are not updated in the form of iteration variable, so y coordinate was always zero in this case. Output, in this case, appeared something like this:

After this, I updated the y-coordinates but again didn't use the "hold off" command and the output was something like this:

Conclusion:

After developing the MATLAB script for this sample problem, one can note that the solution to the problem can be more pragmatic and realistic. It can be done by defining the damping as a function of the surroundings of the system. Here, the damping is taken to be constant, which is not the case in reality. More advanced and complex physical systems can be formulated as a mathematical model and solved by extending this approach.

So I feel, that much more work can be done in this area.