Week 3 - Solving second order ODEs

AIM:

To write a MATLAB code to solve the second order ODE representing the motion of a simple pendulum with damping.

INTRODUCTION:

A simple pendulum is a set up consisting of a bob, string (or a rod) and a pivot. The way the pendulum moves depends on the Newton's Second Law.

Basically, the pendulum swings to and fro about the mean position due to energy conservation principle i.e., kinetic energy to potential energy and vice versa.

DATA GIVEN:

L = 1 `m`

g = 9.81 `m/s^2`

b = 0.05 `(Ns)/m`

m = 1 `kg`

`theta_0` = 0

`dot theta_0` = 3 `(rad)/s`

SECOND ORDER ODE EQUATION:

PROCEDURE:

1. The input parameters are assigned to variables.

2. The initial condition for the angular diplacement and angular velocity are stored in a column vector.

3. The time duration is defined using linspace command.

4. The ODE is solved using an in-built function in MATLAB called ode45.

5. The ode45 function takes in variables involved, calculation function, time duration and the initial condition as the input. It returns the time, angular displacement and angular velocity as a result.

6. The results are plotted with the data obtained.

7. The axis limits and labels are specified.

8. The title of the plot is also specified.

9. An animation of simple pendulum's motion is coded.

10. A pivot point is defined. With the results obtained from ode45 function, the x and y coordinates of the string at each instant is stored in variables x and y using for loop.

11. A line passing through the pivot is plotted to show that the string is pivoted to a ceiling.

12. The string is plotted using the variables x and y.

13. The bob at the end of string is plotted using marker specifier in plot command.

14. The animation is paused at 0.0005s to get a visually appealing animation.

15. The getframe function is used to get each frame.

16. The movie is made by using movie command. A videofile is opened and the frames are written to that file. Once the frames are written, the videofile is closed.

17. The ode_func function called into the main program is calculating the second order ODE equation from the parameters passed into it and returning the value to the main program at the time of calling it.

CODE:

MAIN PROGRAM:

% Solving Second Order ODE for Pendulum Motion

clear all
close all
clc

% Input Parameters
L = 1;              % Length of pendulum in m
g = 9.81;           % Acceleration due to gravity in m/s^2
b = 0.05;           % Damping Coefficient
m = 1;              % Mass of sphere in kg

% Initial Condition
theta_0 = [0;3];

% Time
t_span = linspace(0,20,500);

% Solving ODE using ode45 function
[t,Result] = ode45(@(t,theta) ode_func(theta,t,b,m,g,L),t_span,theta_0);

% Plotting the Result
figure(1)

% Angular Displacement
plot(t,Result(:,1),'r');
hold on;

% Angular Velocity
plot(t,Result(:,2),'b');

% Plot parameters
legend('Angular Displacement','Angular Velocity');
xlabel('Time');
ylabel('Plot');
title('Angular Displacement and Angular Velocity vs Time');
axis([0 20 -10 10]);
hold off;


% Animation
figure(2)

% Counter Variable
ct = 1;

for i = 1:length(Result(:,1))
    % Pivot point of string
    x0 = 0;
    y0 = 0;
    
    % Defining x and y coordinates
    x1 = L*sin(Result(i,1));
    y1 = -L*cos(Result(i,1));
    
    % Line having Pivot
    plot([-2 2],[0 0],'k-','LineWidth',3);
    hold on;
    
    % Axis limit
    axis([-1.5 1.5 -1.5 0.5]); 
    
    % String
    plot([0 0],[0 -1],'k-.');
    plot([x0 x1],[y0 y1],'r-','Linewidth',2);
    
    % Bob
    plot(x1,y1,'o','MarkerSize',20,'MarkerFaceColor','b');
    
    % Pausing the animation
    pause(0.0005)
    
    % Storing frame
    M(ct) = getframe(gcf);
    hold off;
    
    ct = ct + 1;
end 

% Making Movie    
movie(M)
videofile = VideoWriter('Pendulum ODE','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

ODE FUNCTION:

% ODE Function

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

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

% Equation
dtheta2_dt = -(b/m)*theta2 - (g/L)*sin(theta1);

dtheta_dt = [dtheta1_dt;dtheta2_dt];

end

RESULT:

As you can see from the plot, the angular displacement and angular velocity dies out as time goes on. This is because there is a small damping present in the system.

ANIMATION OF MOTION OF SIMPLE PENDULUM:

ERRORS:

I didn't encounter any major errors other than some typing errors here and there.

The above error is rectified by visual inspection.

CONCLUSION:

Hence the second order ODE for simple pendulum motion is solved and the results are plotted and verified.