Week 3 - Solving second order ODEs

Simple Pendulum

Aim:

Solving a second order ordinary differential equation of simple pendulum and showing the effect of damping by an animation.

Objective:

• To convert second order ordinary differential equation to 1^st order ordinary differential equation.
• To create the function using MATLAB in-built function for solving ordinary differential equation.
• To plot the displacement and velocity damping through the time stipulated time of simulation.
• To make the animation of simple pendulum.

Description & Governing equation.

The differential equation is mostly used to define the complex physics phenomenon in the mathematical way. Differential equation is divided in two different classify as ordinary differential equation and partial differential equation. Solving a 1^st order differential is quite easy and require only one constant to define the equation, also it can be also programmed easily. The main effort is taken while solving 2^nd order differential equation. Most of the fluid equation and some structural phenomenon are defined in the 2^nd order differential equation which need to solve eventually to simulated the flow or study the structural behavior.

A simple pendulum consists of mass which is attached by the string of length L from a pivot point P. when displaced by initial force the mass gets moved back and forth performing periodic motion. By applying the Newton’s second law motion we can obtain the equation of motion for simple pendulum which is given as:

d2θdt+bmdθdt+gLsinθ=0

To solve the equation in the MATLAB by using ‘ode45’ the equation needs to be converted from 2^nd order partial differential equation to 1^st order differential equation.

Considerθ=θ1

Then new equation of motion can be written as:

d2θ1dt2+bmdθ1dt+gLsinθ1=0

Consider,

dθ1dt=θ2

Therefore, new equation can be written as:

dθ2dt+bmθ2+gLsinθ1=0

Which leads to:

dθ2dt=−bmθ2−gLsinθ1

The resultant first order equation will be:

θ′1=θ2andθ′2=−bmθ2−gLsinθ1

Code explanation: -

1. The code starts with the inputs which are provided in the problems. Then we make time domain for which the whole code is going to run. It was done by ‘linspace’ command from 0 to 20 seconds with 100 value in between them.
2. After that we provide the initial condition for the angular displacement and velocity which is column vector form, because we need to calculate the both value from the equation of motion for plotting purpose.
3. The code comprises of function which solves the second order differential equation which is taken inside the main code. The main code has another function which calculates the array of data for time and displacement.
4. After calculation plotting is done in which subplot command is used, in the first plot first column is plotted and in second subplot second column is plotted with respect to time.
5. Then we are changing the polar co-ordinate to cartesian co-ordinate because we have used ‘sin’ command instead of ‘sind’ command in function.
6. Then we have created a for loop for which we are creating the frame motion for the simple pendulum and saved in the movie.

The code is given as:

clear all
close all
clc

%Given data
b = 0.05;   %damping coefficient
g = 9.81;   %Acceleration due to gravity
l = 1;      %Length of Pendulum in meters
m = 1;      %Mass 

%Time duration for simulation 
t_span = linspace(0,20,100);
%Providing Initial Condition as Angular displacement & Velocity
y_initial = [0;3];

%Calculation using ODE45
a = @(t,y)ode_second_order_pendulum(y,b,m,l,g);
[t,y]= ode45(a,t_span,y_initial);

%Plot of Angular displacement & Velocity with respect to Time function.
figure(1);
grid on
subplot(1,2,1)
plot(t,y(:,1));
title('Displacement vs Time');
xlabel('Time in seconds');
ylabel('Angular displacement in rad');
subplot(1,2,2)
plot(t,y(:,2));
title('Angular velocity Vs Time')
xlabel('Time in seconds');
ylabel('Angular velocity in rad/sec');
grid off
saveas(figure(1),'subplot','jpg')

%Polar Co-ordinate To Cartesian Co-ordinate
y1 = -1*cos(y(:,1));
x1 = 1*sin(y(:,1));

%Video for Pendulum motion
for i = 1:100
    figure(2)
    plot([0 x1(i)],[0 y1(i)],'linewidth',3,'color','r');
    hold on
    plot(x1(i),y1(i),'marker','o','linewidth',10,'markerfacecolor','g');
    plot([-1.2 1.2],[0 0],'linewidth',3,'color','b');
    axis([-1.2 1.2 -1.2 1.2]);
    title(['Time: ', num2str(t_span(i))]);
    M(i) = getframe(gcf);
    close(2)
end

movie(M)
videofile = VideoWriter('singlependulum.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M)
close(videofile)

The Function is given as:

function ydot = ode_second_order_pendulum(y,b,m,l,g)
ydot = [y(2);-(b/m)*y(2)-(g/l)*sin(y(1))];
end

Inference of Plot:

As we can see from the plots the displacement is getting reduce as we move further in to time because of damping factor. From the plot we can clearly observe displacement is starting from 0 as initial condition. In the angular velocity vs time plot it has started from 3 rad/sec and then further get decrease.