Second order ODEs:- Equation of motion of a simple pendulum with damping.

1. AIM

• Write a program that solves the ODE for equation of motion of a simple pendulum with damping.

2. OBJECTIVE

• To understand the transient behavior of various system using matlab, example “equation of motion of a simple pendulum with damping.”
• Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

3. Governing Equation:-

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

In the above equation,

g = acceleration due to gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

4. BODY OF THE CONTENT

• Motion of Pendulum:- Above equation represents the equation of motion for Pendulum with damping.
• In Engineering application many transient behavior of system can be represent by the ODE.
• Motion of Pendulum can be represent by the Sencond order ODE.(Ordinery Differential Equation.)
• We are using “ODE45” command, to evaluate this equation.
• The main problem associate with this command is that, it cannot evaluate the second order ODE, hence we are dividing second order ODE in two different First order ODE.
• For this we have to do following assumption

Let ϴ = ϴ₁

`(dθ)/dt=(dθ_1)/dt= θ_2`

By using these substitution, we can reduce this second order ODE in first order ODE

1. `(dθ_1)/dt= θ_2`

2. `(dθ_2)/dt=-b/m dθ/dt-g/L sinθ`

Syntax:- Main Program

clear all;
close all;
clc;

%Input Parameters
b=0.05;         %Damping coefficient
l=1;            % Length of string
g=9.81;         % acceleration due to gravity
m=1;            % mass of Ball(object)

%Initial condition
theta_0=[0;-3]; %Initially Angular displacement is zero & Angular momentum is 3 radian/s

%time span for evaluation
t_span=linspace(0,20,100);  % Evaluating equation for 0 to 20 sec.

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

%Pendulum tragectory path 
for j=1:50
    u=results(j,1);
    x2_end(j)=-l*cos(u);
    y2_end(j)=-l*sin(u);
end
    %PLOTTING FOR ANGULAR DISPLACEMENT VS VELOCITY
    subplot (2,1,1)
    plot(tSol,results(:,1))
    hold on
    plot(tSol,results(:,2))
    axis([-0 20 -4 8]) 
    legend({'angular dislacement','angular velocity'})
    ylabel('Plot')
    xlabel('time') 
    title('ANGULAR DISPLACEMENT VS VELOCITY' )
    
%Loop for plotting Pendulum
for i=1:length(tSol)
    m=results(i,1);
%PENDULUM PLOTTING
    x1_end=-l*cos(m);   %Ending 'X' co-ordinates of pendulum
    y1_end=-l*sin(-m);  %Ending 'Y' co-ordinates of pendulum
    subplot (2,1,2)
    Pendulum_string = plot([-1 1],[0 0],[0 y1_end],[0 x1_end],y1_end,x1_end,'o','MarkerSize',30,'markerfacecolor','k')
    hold on
    Central_Line=plot([0 0],[0 -1.5], 'LineStyle','--')
    Pendulum_path=plot(y2_end,x2_end,'.')
    axis([-1.5 1.5 -1.4 0.3]) 
    title('Pendulum motion')
    pause(0.1)
    hold off
    %Making animation
   M(i)=getframe(gcf);
       if i<length(tSol)
       delete(Pendulum_string)
       delete(Central_Line)
    end
end

%Saving Animation
movie(M)
videofile=VideoWriter('Pendulum_motion1.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

Function syntax:- ode_func

function [dtheta_dt]=ode_func(t,theta,b,g,l,m)
k1=b/m
k2=g/l
theta1=theta(1);
theta2=theta(2);
dtheta1_dt=theta2;
dtheta2_dt=-k1*theta2-k2*sin(theta1);
dtheta_dt=[dtheta1_dt;dtheta2_dt];
end

Output:-

Syntax Explanation:-

Function:-

• Program “ode_func” is basically function which defines output as two first order ODE, which we derive earlier.
• We are doing some more substitution to make equation easy to type in program(like k1 & k2)
• In 6^th and 7^th line we are defining to variables dtheta1_dt & dtheta2_dt for the above two equations.
• Then we are saving these two output as a single array, to evaluate using “ODE45” command.
• In this way we are getting two outputs by solving these ODE, first is ‘dtheta1_dt’ which is angular displacement & second ‘dtheta2_dt’ which is angular velocity.

Main Program Explanation:-

clear all;
close all;
clc;
 
%Input Parameters
b=0.05;         %Damping coefficient
l=1;            % Length of string
g=9.81;         % acceleration due to gravity
m=1;            % mass of Ball(object)

• In These lines we are defining input parameters for pendulum, i.e Damping coefficient, Length of string, acceleration due to gravity & Mass of object

Setting initial conditions

%Initial condition
theta_0=[0;-3]; %Initially Angular displacement is zero & Angular momentum is 3 radian/s

%time span for evaluation
t_span=linspace(0,20,100);  % Evaluating equation for 0 to 20 sec.

• Now to evaluate ODE, we are defining Initial condition for Angular Displacement to 0, angular velocity to 3 radian/s.
• In 14^th line we are defining time span to evaluate ODE, 0 to 20 s.

Solving ODE

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

• In 18^th line we are solving these two ODE using “ODE45” command.
• The second argument after ODE45 command is a function which is to be evaluate, in this case it is “ode_func”
• Then we are setting initial condition for dependent & independent variable respectively.
• In output array we have define two variables, first tSol provide solution for time span, and second “results” contains two coloms, first for Angular Displacement & second for angular velocity, as we defines in function.

Plotting Results:-

• We have to Plot two figures, first ANGULAR DISPLACEMENT VS VELOCITY graph, & second is simulation for motion of pendulum.

%PLOTTING FOR ANGULAR DISPLACEMENT VS VELOCITY
    subplot (2,1,1)
    plot(tSol,results(:,1))
    hold on
    plot(tSol,results(:,2))
    axis([-0 20 -4 8]) 
    legend({'angular dislacement','angular velocity'})
    ylabel('Plot')
    xlabel('time') 
    title('ANGULAR DISPLACEMENT VS VELOCITY' )

• For first graph, we have use results which we are getting solution of ODE, and plot is single figure.(from line 26 to 35)
• At line 28 we are plotting graph for “Angular displacement” & using hold on command in 29 line we are plotting graph for “Angular Velocity”.
• Line 31 to 35 sets, the different parameters of graph like axis setting, legends, axis labels and titles.
• For second graph we want to plot simulation for motion of pendulum.

%Loop for plotting Pendulum
for i=1:length(tSol)
    m=results(i,1);
%PENDULUM PLOTTING
    x1_end=-l*cos(m);   %Ending 'X' co-ordinates of pendulum
    y1_end=-l*sin(-m);  %Ending 'Y' co-ordinates of pendulum
    subplot (2,1,2)
    Pendulum_string = plot([-1 1],[0 0],[0 y1_end],[0 x1_end],y1_end,x1_end,'o','MarkerSize',30,'markerfacecolor','k')
    hold on
    Central_Line=plot([0 0],[0 -1.5], 'LineStyle','--')
    Pendulum_path=plot(y2_end,x2_end,'.')
    axis([-1.5 1.5 -1.4 0.3]) 
    title('Pendulum motion')
    pause(0.1)

    %Making animation
   M(i)=getframe(gcf);
       if i<length(tSol)
       delete(Pendulum_string)
       delete(Central_Line)
    end
end

• For plotting results, we use angular displacement (derived by solving ODE).
• To plot pendulum for different values of displacement, we are using “For loop”, and plotting individual position of pendulum w.r.t. time.
• For pendulum string we need two co-ordinates, first end is fixed at co-ordinate at point (0,0) and second co-ordinate is varying with each iteration (x1_end, y1_end).
• The object(ball) is fixed to second end of string at((x1_end, y1_end).)
• From line 44 to 47, we plot three component of pendulum i.e. string, Central line(reference line) & trajectory path.
• Line 48 & 49 defines properties of graph, axis values & title of graph.
• In line 50, we are pausing loop for 0.1seconds, to get the animation.
• The problem using “For loop” in syntax is that, all previous values of graph are changes in each iteration, hence we are not getting correct animation. (This problem is explane in Problem section.)
• In line 55 to 59 we are using “if ” & “delete” command to solve this issue.
• From line 61 to 66, we are saving this simulation as a video file, using “movie” command.

Errors observed during Coding:-

• Inverted Prendulum.- For positive value of String length “l” we are getting inverted pendulum graph as shown below.

This problem is solved by using negative value of length of Pendulum.

x1_end=-l*cos(m);   %Ending 'X' co-ordinates of pendulum
y1_end=-l*sin(-m);  %Ending 'Y' co-ordinates of pendulum

2.Pendulum motion Trajectory Line:- We need to show pendulum motion trajectory line as well, but by the execution of every iteration of loop previous points are vanished, hence we are getting single point instead of line for motion.

To solve this we define separate loop which stores the value of angular displacement define in variable “results” and plot in advance.

But the problem insisted because every time we execute this loop, the previously plot trajectory plot is overwrite by figure, and to stop if we use “hold on” command we are getting following error

Hence to solve this, instead of using “hold on”, we use delete command, by which we are deleting only pendulum for every iteration. & keeping trajectory path as it is.

       if i<length(tSol)
       delete(Pendulum_string)
       delete(Central_Line)

5. CONCLUSION

From observation of output graph, we conclude that

• The more we increase the damping coefficient, the less time system will take to stabilize..
• The more the initial values ( Transients), system will take more time to stabilize.

Animation for Motion of Pendulum:- https://youtu.be/Q6-A397eZ98

Referances:-

• Simple Harmonic Oscillation: https://www.youtube.com/watch?v=gnD8Se92hfk
• Motion of Pendulum:- https://www.youtube.com/watch?v=WPa5IgLgDyQ&t=384s