Solving second order ODEs

Solving Second Order ODEs in MATLAB

1.AIM : To simulate the transient behaviour (motion) of a simple pendulum in MATLAB.

2.OBJECTIVE :

• To write a program in MATLAB to solve ODE (ordinary differential equation) for equation of motion of a simple pendulum with damping.
• To simulate the motion of the pendulum between time period of 0 – 20 sec.
• To create the animation of the motion of pendulum by getting its position as a function of time.

3.THEORY :

3.1 Simple Pendulum : - A simple pendulum is one which can be considered to be a point mass suspended from a string or rod of negligible mass. It is a resonant system with a single resonant frequency.

3.2 Pendulum Motion : - The motion of a simple pendulum is like ‘simple harmonic motion’.

fig.1: Animation of the motion of a Simple Pendulum

3.3 Period of Simple Pendulum : - A point mass hanging on a massless string is an idealized example of a simple pendulum. When displaced from its equilibrium point, the restoring force which brings it back to the center.

3.4 Pendulum Geometry : - The period of a simple pendulum for small amplitudes 𝚹 is dependent only on the pendulum length and gravity.

fig.2: Simple Pendulum (free body diagram)

3.5 Pendulum equation : -

• In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum.
• The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a 2^nd order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

4. CALCULATION:

The ODE which represents the equation of motion of simple pendulum with damping is –

In the above equation,

g = gravity in m/s²

L = length of the pendulum in m

m = mass of the ball in kg

b = damping coefficient.

Now,

Solving 2^nd Order differential equation –:

Let’s say,

`θ = θ_1`

Then,

`"dθ" /dt = (dθ_1)/dt = θ_2` --------------- ①

`(d^2 θ)/(dt^2 ) = (d^2 θ_1)/(dt^2 ) = d/dt ((dθ_1)/dt) = (dθ_2 )/dt` --------------- ②

Now, going back to our original 2^nd order ODE equation,

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

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

`rArr` `(d^2 θ_2)/(dt^2 ) + "b" /"m" (dθ_2)/dt + "g" /"L" sin⁡θ_1 =0`

So, in original words to solve this 2^nd order ODE; what we have done is,

we have created two ODEs –

1st ODE	2nd ODE
`(dθ_1)/dt = θ_2`	`(dθ_2)/dt = -"b" /"m" θ_2 - "g" /"L" sin⁡θ_1`

So, this is what I meant as breaking down a 2^nd order ODE into a series of 1^st order ODEs, and that is why we did this transformation.

• So, the idea is if we have a 3^rd order ODE, then we can break it down into a system of three 1^st order ODEs.

Then, we are going to say that –

`d/dt [ [θ_1], [θ_2] ] = [ [θ_2], [-"b" /"m" θ_2 - "g" /"L" sin⁡θ_1] ]`

'`θ` '

Above one is a solution vector. [We are calling it as a solution vector because it contains two quantities].

So, now when we are trying to solve this equation with programming; it can be any of the ODE solver through MATLAB, we need to provide the equation on this particular format. And that is what we are going to do at next –

5. BODY OF CONTENT :

Function Code used –

% A 'Function program' which will be called in the main program for solvin
% the ODEs

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

Main Code –

% A program for solving second order ODEs and to simulate its motion

clear all
close all
clc

% Inputs
b = 0.05;                                                                    % Damping coefficient
g = 9.81;                                                                    % Acceleration due to gravity (m/s^2)
l = 1;                                                                       % Length of string (m)
m = 1;                                                                       % Mass of bob (kg)

% initial condition
theta_0 = [0 ; 3];                                   % Initializing the condition as '0' for displacement and '3' for velocity

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

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

% Plotting graph of displacement and velocity w.r.t. time
figure(1)
plot(t,results(:,1))                                                        % Plotting time Vs all rows in 1st column
hold on
plot(t,results(:,2))                                                        % Plotting time Vs all rows in 2nd column
xlabel('Time t(sec)')
ylabel('Plot of Displacement and Velocity')
title('Damped Harmonic motion of a Simple Pendulum')
legend('Displacement','Velocity')
grid on
hold off


%Plot for different positions of Pendulum for creating an Animation
% Iteration using 'for loop' to create animation of motion of the pendulum

loop_counter =1;
for i = 1:length(results(:,1))
    x_fixed = 0;                                                            % fixed end co-ordinates
    y_fixed = 0;
    
    x_freeend = l*sin(results(i,1));                                        % free_end co-ordinates
    y_freeend = -l*cos(results(i,1));                                       % Got the equations from the geometery of the pendulum
    
    % Plotting the pendulum for motion
    figure(2)
    plot([-1 1],[0 0],'linewidth',8,'color','black')                                   % Plot for Fixed line
    hold on
    plot([0 0],[0 0],'o','markerSize',7,'markerfacecolor','yellow')                    % Plot for hinge point in fixed line
    hold on
    plot([x_fixed x_freeend],[y_fixed y_freeend],'linewidth',3,'color','green')        % Plot for pendulum string as line with green color
    hold on
    plot(x_freeend, y_freeend,'marker','o','markerSize',20,'markerFacecolor','red')    % Plot for mass bob fixed to string
    grid on
    axis([-1.5 1.5 -1.5 1.5])                                                          % Setting the aaxes limits in plot
    title('Animation of transient behaviour of a Simple Pendulum')
    pause(0.03)
    hold off
    
    M(loop_counter) = getframe(gcf);                                   % Capturing figure and plot in frames and saving in Variable M
    loop_counter = loop_counter + 1;
    
end


% Creating animation by compiling the different frames and saving the movie
% in file
movie(M)                                                              % Creating movie for all the saved frames of plots in Variable M
videofile = VideoWriter('Simple Pendulum.avi','uncompressed avi');
open(videofile)
writeVideo(videofile,M)                                               % Writing the video file to movie frames saved to Variable M
close(videofile)

CODEs OVERVIEW :

• Initial inputs are taken according to the provided data as follows –
  1. L = 1 metre
  2. m = 1 kg
  3. b = 0.05
  4. g = 9.81 m/s²
  5. time period = 0-20 sec
  6. angular displacement at time (t=0) = 0
  7. angular velocity at time (t=0) = 3 rad/sec

• An array of initial theta was created to define initial angular displacement and initial angular velocity.
• An array of time span was created to store 500 equally spaced intervals between time period of 0 – 20 sec.
• To solve the 2^nd order ODE, a function named ode_function which stores the above matrix data in it.
• The given 2^nd order ODE was solved using ode45 command followed by proper syntax of the calling function ode_function and the result obtained was used to plot the graph of variation of angular displacement and angular velocity w.r.t time.
• Now for creating the animation of obtained harmonic motion of pendulum, ‘for loop’ was used and a loop counter was initiated. In this execution, the one end of the string of pendulum was fixed and the other end moves according to the solution of the ODE solved and bob was made by using ‘marker’ command and was plotted for each time interval of the time span.
• All the frames were compiled using ‘getframe’ command and ‘movie’ command was used to create the animation in MATLAB.

6. RESULTS :

The plot generated is -

fig.3: Plot for displacement and velocity Vs time for damping (i.e the results with the time period defined in the code respectively).

• The above plot shows that the motion of pendulum id decreasing gradually. This is due to the damping effect which exponentially decreases the amplitude and velocity of the motion of pendulum.

Animation Video –

7. CONCLUSION :

From this, we get to understand the equation of motion with damping of a simple pendulum and also created a program for solving ODE for pendulum using ode45 equation solver for the transient behaviour of the system and animated the motion of the pendulum when the displacement changes with respect to time.

8. ERRORS :

8.1 While performing animation, the pendulum started moving above instead of hanging downward from the fixed wall.

This happens due to error in code at –

Solution: - solved after assigning the correct +ve , -ve sign todefine the position of pendulum correctly,

8.2 The animation was not stable.

In the created animation, the motion of pendulum was not clear visual due the motion area completely filled with a large number of strings.

This happens due to error in code at –

Solution: - Solved after adding the ‘hold off’ command where requires to make the animation of motion of pendulum to be stable and traceable as,