Solving and animating ODE of a pendulum

• Objective of this project is to solve the given ODE of second order, which represents the motion of a pendulum, and simulate its motion as time passes using MatLab.

Code for the above project:

Main program:

clear all
close all
clc

%Constant values
b=0.05;
g=9.81;
l=1;
m=1;

%initial condition
theta_0=[0;3];

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

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

figure(1);
plot(t,results(:,1));
hold on;
plot(t,results(:,2));
xlabel(\'Time\');
ylabel(\'Plot\');

%Animation
%Constant values
r=1;
x1=0;
y1=0;
ct=1;

for i=1:length(results(:,1))
    THETA=results(i,1);
    x2=r*sin(THETA);
    y2=-r*cos(THETA);
    figure(2);
    plot([-r r],[0 0],\'linewidth\',5,\'color\',\'k\');
    hold on
    plot([x1 x2],[y1 y2],\'linewidth\',3);
    plot(x2,y2,\'-o\',\'markers\',20,\'markerfacecolor\',\'r\');
    grid on;
    axis([-1.5 1.5 -1.5 .5]);
    pause(.0003);
    hold off;
    M(ct)=getframe(gcf);
    ct=ct+1;
end

movie(M)
videofile=VideoWriter(\'Pendulum_animation\',\'Uncompressed AVI\');
open(videofile)
writeVideo(videofile,M)
close(videofile)

The Function:

%ODE function
function [dtheta_dt]=ode_fun(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

• Explanation for the above code:

After opening a new MatLab file, the following steps are followed:

Step 1: The first three lines consist of unsual commands to close and clear all the windows, variables of previous MatLab files.

Step 2: The initialisation of the required constant values like \'b\',\'g\',\'l\',\'m\' to solve the ODE are done.

Step 3: The initial conditions given are initialised to a column array \'theta_0\',i.e. 0 for angular displacement, 3 for angular velocity.

Step 4: The time up to which the integration has to take place is initialised to the array \'t_span\'.It consists of 500 values from 0 to 20, which means that we get solution at each of these points.

Step 5: Now the inbuilt function \'ode45\' is used to solve the given ODE, the syntax is:

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

*Firstly, anything after \'@\' enclosed within \'()\' signifies the variables in the given ODE, here it is t which is the time and theta which is the angle thro\' which the pendulum has moved.

*Then the ODE must be provided to the \'ode45\' in solver understandable format, for this we define the function \'ode_fun\'.

*Explanation of the function \'ode_fun\':

1. The function takes \'t\',\'theta\',\'b\',\'g\',\'l\',\'m\'; as the input.

2. The array \'theta\' has two values the two values are initialised to the respective variables.

3. From the given second order ODE, we are actually splitting it into two first order ODE, for the solver could solve the ODE. So, the theta in the given ODE is taken as \'theta1\'and then the below is assumed;

\'dtheta1_dt=theta2;\'

Then, the above is substituted in the given ODE to get dtheta2_dt;

\'dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);\'

Then the above are stored in the following array,

\'dtheta_dt=[dtheta1_dt;dtheta2_dt];\'

The array is the return value of the function.

*Now that we have got the ODE in the required format, the time for which is the ODE must be solved is given as the input parameter to the \'ode45\' function, then the initial conditions are given.

*Now the function solves the ODE and stores the value of angular displacement(theta1) and angular velocity(theta2) in the \'results\' field and the time taken is stored in \'t\' field.

Step 6: As the solutions have been found, they are plotted in order to make the result more understandable. This is done by using plot command, and the plot is made between time and the results obtained.

Step 7: The aim to animate the pendulum based on the result obtained is done by using the basic algorithm for animation.

*In order to do that the length of the pendulum(\'r\'), the coordinates of the point at which the pendulum is fixed(\'x1\',\'y1\'), and the counter are initialised.

*As we want to simulate the motion of a pendulum we plot the line joining the fixed point and the other end of the pendulum for each angular displacement,for this we use a for loop. The respective formula is used to calculate the moving point for each angle.

*The bob of the pendulum is made by using markers option of the plot command.

*Then the plots are captured as a movie and then converted into a video file in the next few steps.

Step 8: The video is then uploaded in YouTube.

• Errors faced:

NIL

• Results:

1.Time Vs angular displacement and angular velocity:

The blue plot is between Time and angular displacement.

The red plot is between Time and angular velocity.

As the value of damping coefficient is smaller the effect of damping is not very high in the above case.

2.YouTube link for Pendulum animation:

https://youtu.be/6mDbAFPooIw

The above is the animation of pendulum for the given scenario.