Week 3 - Solving second order ODEs

AIM:

Solve 2nd order ODE using MATLAB.

OBJECTIVES:

• To solve 2nd order Simple Pendulum ODE and plot amplitude vs. time graph.
• To simulate the motion of Simple Pendulum.

GOVERNING EQUATIONS:

The ODE of the motion of simple pendulum with damping

`(d^2 θ)/(dt^2 )+ b/m dθ/dt+ g/L sin⁡θ=0` ----eq. 1

m = mass at the swinging end of the pendulum (kilograms)

g = acceleration due to gravity (meter/second2)

L = length from the swivel point to the center of mass (meters)

θ = angle between the string position to the string position at rest (radians)

b = damping co-efficient.

To solve this, we need to reduce it into lower order ODE by substitution. thus,

let `θ = θ1`

So, we can say that `(dθ)/dt= (dθ_1)/dt= θ_2` ---eq. 2

`(d^2 θ)/(dt^2 )=(d^2 θ_1)/(dt^2 )=(dθ_2)/dt` --- eq. 3

Putting (2) and (3) in (1), we get :

`(d^2 θ_1)/(dt^2 )+b/m (dθ_1)/dt+ g/L sin⁡〖θ_1 〗=0` that is

`(dθ_2)/dt+b/m θ_2 + g/L sin⁡〖θ_1 〗=0`

Since it is reduced to lower order, by rearranging terms, we can say that,

`(dθ_2)/dt =-b/m θ_2 - g/L sin⁡〖θ_1 〗`

Above reduced ODE is used to create function for main programme using ode45 function.

%Ode_func Simple pendulum

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

LINK: https://web.math.pmf.unizg.hr/nastava/nmf/Njihalo.pdf

BODY:

Programming is done in Matlab.

%Programe_that_solve_ODE_of_motion_of_simple_pendulum

%"clear all" is basically used to clear/delete the variables created in Workspace.
%"close all" is used to close all figures.
%"clc" - clears the command window
clear all
close all
clc

%Inputs
b = 0.1; %Damping co-efficient
g = 9.81; %Accleration due to gravity 'm/s^2'
l = 1; %Lengeth of the pendulum 'm'
m = 0.2; %Mass of ball 'kg'

%Initial condition
theta_0 = [0,3];

%Time points
t_span = linspace(0,20,1000);
%%linspace(i,j,n) is used to create 'n' values between i and j

%Slove ODE
%'ode45' is inbulit matlab function which allow us to calculate Ordinary or
%partial differential euqation
%with the help of function, time span and initial condition.
[t,results] = ode45(@(t,theta) ode_func_simple_pendulum(t,theta,b,g,l,m),t_span,theta_0);

%Plotting
figure(1)
%Figure command open figure page for plots.
plot(t,results(:,1))
grid on
hold on
plot(t,results(:,2))
xlabel('Time')
ylabel('Amplitude')
legend({'Displacement (rad)','Angular speed (rad/sec)'});

ct = 1;
for i = 1:length(results(:,1))
    x0 = 0;
    y0 = 0;
    x1 = l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    figure(2)
    plot([x0,x1],[y0,y1],'linewidth',3,'Color','b');
    %"plot(X,Y)" creates a 2-D line plot of the data in Y versus the corresponding values in X.
    hold on
    %hold on retains plots in the current axes so that new plots added to the axes do not delete existing plots.
   
    plot(x1,y1,'-o','markers',10,'markerfacecolor','r');
    %marker symbol in the line-specification input argument, such as plot(x,y,'-O')
    axis([-1.5 1.5 -1.5 1.5]);
    %Is used for axis(limits) specifies the limits for the current axes
    hold off;
    %hold off sets the hold state to off so that new plots added to the axes clear existing plots and reset all axes properties.
    pause(0.0003);
    M(ct) = getframe(gcf);
    %Getframe captures the current axes as it appears on the screen as a movie frame
    ct = ct+1;
end
videofile = VideoWriter('Simple_Pendulum.avi','Uncompressed AVI');
%VideoWriter object to create a video file from an array
%VideoWriter(filename,profile) additionally applies a set of properties tailored to a specific file format (such as 'MPEG-4' or 'Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
%writeVideo(v,img) writes data from an array to the video file associated with v.

%You must call open(v) before calling writeVideo.
close(videofile)

Code explaination:
1. Set physics parameters of the pendulum like mass, length, gravity and damping coefficient.
2. Set the initial conditions and timespan for which you want to run the simulation.
3. To use the ode solver, use ode45 function with main inputs as (t,theta,b,g,l,m) while you use ode_func (the one we created earlier) with inputs as (t,theta,l,b,g,m) and t_spa
4. These results are stored in an array of [t,results]
5. Since results contain 2 components, angle displacement and velocity, to plot both the graphs, use the (:,1) and (:,2) property of MATLAB and plot those graphs.
6. To create animation, Use a for loop so that we can get X and Y co-ordinates using trigonometry, Here, x0 and y0 are the hinge points. Since they are (1,1) we have subtr
7. Use the pause () function to give a pause for every looping of values.
8. Use movie making commands as per code to create a animation video.

Amplitude vs. Time

From the graph we can understand clearly that when there is damping in the motion, with time amplitude decrease.

Pendulum

Link for pendulum animation without damping.

https://youtu.be/Q9dYu55BcyE

CONCLUSION:

We have successfully solved the second order ODE of Simple Pendulum.