Week 4.1 - Solving second order ODEs

Solving second order ODEs

OBJECTIVE:

To generate a plot for a simple pendulum and create the animation for the same.

INTRODUCTION:

A simple pendulum is an arrangement in which a weighted mass is suspended with the help of a string, which is fixed at one end and free at another end. This weighted mass free to move in space when provided with an external force. When the mass is displaced sideways from its equilibrium position it is subjected to restoring force due to gravity, which tries to pull it back to the equilibrium position. This restoring force tends to oscillate the mass about its equilibrium position moving it back and forth. This motion of the pendulum will keep going until the restoring force is less than the gravity and the frictional force due to the surrounding atmosphere.

The motion of the pendulum is a circular arc. Every motion of the pendulum will make an angle θ with the vertical axis. This angle θ defines the amplitude of the oscillatory motion in the vertical plane. greater the angle higher are the oscillations and bigger is the arc created by the bob. We have to solve for the angle θ to trace the motion path of the bob, to find the change in velocity according to the time.

Fig. 1

GOVERNING EQUATION:

The governing equation that represents the motion of the pendulum in air is

In the above equations

b = Damping coefficient

m = mass of the bob

g = acceleration due to gravity

l = length of the string

This is a second-order differential equation, which first needs to be simplified for the angle θ and can be directly solver in MATLAB for the Time and velocity plot of the motion. Later this output can be used for generating the animation of the moving pendulum.

To solve this equation, we first need to convert it into a first-order ODE. Which then can be used as an input function in the code, and by using the ode45 command we can compute the results by giving proper initial conditions and input parameters.

PROGRAM PROCEDURE:

• First, we entered our input variables that are b,m,l and g. Along with that, we had our initial conditions for displacement and velocity i.e. 0m and 3 rad/sec respectively.
• Next, we introduced the time vector, as we have to simulate the motion between 0-20 seconds we will be using the linspace command to generate the time vector.
• We have simplified the second-order ODE into 2 first-order ODE and included two variables for them θ1 and θ2. These θ1 and θ2 are the two variables that will be defining our output vector of the ODE.
• In this function, we have included the θ1 and θ2 angles, the simplified form of our second-order ODE, and the solution matrix.
• With using the command ode45 we will introduce the ode solver into our code. (t,theta) are the variables defining the ode, second_orderODE_pendulum(t, theta,b,g,l,m) is the function that we created, t-span is the time vector that will change according to time, theta_0 are the initial conditions for the function.
• This solver will be solving the ode as per the inputs and initial conditions that are provided.
• The next step will be plotting the velocity and displacement profile. The x and y-axis will have time and velocity respectively. By using simple plotting command, we can easily plot the output of the ode that is solved.
• Now to generate the animation of the plot, first, we will define a variable theta2 to the results that are generated. Then this variable will be introduced into a FOR LOOP which will later be assigned the steps to a new variable theta_0. And before FOR LOOP ct is considered as 1.
• By giving the conditions for the string i.e. x0 and y0 is the fixed point of the string and the other point will be defined by x1 and y1 which will change w.r.t time.
• For plotting the animation, first, we will be plotting the fixpoint of the animation, next will be the connecting rod that is attached to the fixed point and the ball at the other end, at the last we will be using the marker command to plot to ball.
• For generating the movie first we will be creating an array 'M' which will be storing the frames generated, ct is the counter that will be incrementing the frame count.
• VideoWriter will help us to generate the avi file for the animation. First, we will be opening the video file, writeVideo will write the videofile and using array M, finally, we will be closing the video.

MATLAB CODE:

clear all
close all
clc

%Input parameters
b = 0.05;
g = 9.81;
l = 1;
m = 0.1;

% Initial condition(displacent & velocity)
D_v = [0;3];

% Time points
t_span = linspace(0,20,600);

% Solveing  using ODE45 command
[t, results] = ode45(@(t,theta) Pendulum_2nd_order_ODE(t,theta,b,g,l,m),t_span,D_v);


plot (t,results(:,1))                       % Time(Displacement) of pendulum plot
hold on
grid on
plot (t,results(:,2))                       % Velocity of pendulum plot
xlabel("Time")
ylabel("Velocity")
legend("Displacement","Velocity")
title('Time vs Velocity')



ct = 1;
%creating animation for the pendulum
theta2 = results(:,1);
for i = 1:length(theta2)          % FOR LOOP
    theta_0 = theta2(i);
    
    % Base stanrting point 1 of of pendulum 
    x0 = 0;                      % Support of pendulum Link 1
    y0 = 0;                      % Support of pendulum Link 1
    
    % Ending point of link 1 and staring point link 2(ball) with
    % length(l) and angle of theta_0 
    x1 = l*sin(theta_0);         % Support of pendulum Link 2
    y1 = -l*cos(theta_0);        % Support of pendulum Link 2
    
    %ploting the pendulum 
    figure(2)
    plot([-1 1],[0 0],'color','g','linewidth',6);         % supporting arm plot
    hold on
    plot([x0 x1],[y0 y1],'color','b','linewidth',3);      % connecting rod plot
    hold on 
    plot (x1,y1,'o','markers',15,'markerfacecolor','r');  % End point Ball plot
    hold off
    grid on
    xlabel("Displacement")
    ylabel("Velocity")
    title("Pendulum Motion Using 2nd order ODE")
    axis([-1.5 1.5 -1.5 0.5]);                            % x and y axis values i.e. [x(start) x(end)] [y(start) y(end)]
    pause (0.02);                                         % Pause time 
    M(ct) = getframe(gcf);                                % Movie creation
    ct = ct+1;
end
movie (M)
videofile = VideoWriter('Pendulum_animation.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M)
close(videofile)



FUNCTION OF 2ND ORDER ODE PENDULUIM:

function[dtheta_dt] = Pendulum_2nd_order_ODE(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

PLOT:

Plot. 1

Plot.1 shows the Time vs Velocity plot where the velocity(v) starts from 3 m/s from y-axis and the time(t) is from 0 to 20 sec.

In plot.1 The green line shows the displacement of the pendulum.

Pendulum of 2nd order ODE:

Plot.2

Plot.2 shows the movement of the pendulum at certain angle or displacement with particular velocity.

The below link is for pendulum motion

https://youtu.be/thBzWDFhGl8

ERRORS Occured during program:

The above shown error is debug by adding the , operater before t_span command.

The above pendulum error is resolved by hold off command.

CONCLUSION:

As seen in the output plot the displacement of the pendulum was slowly getting reduced because of the damping provided. The amplitude of the velocity is also decreasing as the time passes by.

The animation generated is also supporting the results that we concluded, as the time proceeds the ball slow down and the amplitude of the swing reduces.