Week 3 - Solving second order ODEs

Aim :- Write a program to that solves the following ODE and perform the oscillation of simple pendulum.

                      

 

Theory :-

1. The Differential equation is a mathematical equation that relates some functions with its derivatives. In application the function represents a physical quantity and the derivatives represent its rate of change and the differential equation represents the relation between the two. The highest order derivative sets the order of the equation.
2. ODE refers to an ordinary differential equation in which the highest power of equation is one. i.e. The order of the equation is one.
3. In Engineering, ODE is used to describe the transient behaviour of a system. A simple example is a pendulum.
4. The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a Second Order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.
5. Now, we will understand the motion of simple pendulum as under;

 

The equation of motion of a simple pendulum is based on Newton second law of motion.

A simple pendulum consists of a ball (point-mass) m hanging from massless string of length L and fixed at pivot point P.

 

                                           

                                                      Image – Oscillation of Pendulum

 

The above is a motion of a simple pendulum. The motion of a pendulum can be defined as under by splitting bob motion into cosine and sine component.

 

                                                                       

                                                  Fig.:- Splitting Bob motion in sine and cosine component.

 

As a result, we can find the position of bob at any given time. When displaced to an initial angle and released, the pendulum will swing back and forth with periodic motion.

By applying Newton's second law for rotational systems, the equation of motion for the pendulum may be obtained as

which is;

The below equation is, second order differential equation of a simple damped pendulum:

 

where,

b = damping coefficient(dimensionless)

m = mass(kg)

g = gravity(m/s )

L = length of pendulum(m)

 

Let us Consider

  -------------------(1)

We know,

Put  

Therefore,

After substituting the values, we get

   -----------------(2)

Rearranging equation 1 and 2 in matrix form.

  -----------------(3)

 

Main Program :-

% Write a program for animation of pendulum motion using matlab
clear all;
close all;
clc;

% Inputs

% mass of ball in kg
m=1; 
% Length of string in m
l=1; 
% Gravity in m/s^2
g=9.81;
% Damping coefficient
b=0.05;

% Input condition
theta_0=[0;3];

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

% Results
[t, results] =ode45(@(t,theta) ode_function(t,theta,b,g,l,m),t_span,theta_0);

plot(t,results(:,1),'color','r');
hold on;
plot(t,results(:,2),'color','g');
hold off;
plot(results(:,1),results(:,2));
xlabel('Time');
ylabel('Angular velocity vs Anglular Displacement');

ctr=1;
for i=1:length(results(:,1))
    x0=0;
    y0=0;
    x1=l*sin(results(i,1));
    y1=-l*cos(results(i,1));
    plot([-1 1],[0 0],'linewidth',8,'color','b');
    hold on;
    line([x0 x1],[y0 y1],'Linewidth',4,'color','k');
    hold on;
    plot(x1,y1,'-o','markers',20,'markerfacecolor','m');
    grid on;
    axis([-1.5 1.5 -1.5 0.5]);
    pause(0.01);
    hold off;
    A(ctr)=getframe(gcf);
    ctr=ctr+1;
end

movie(A)
videofile = VideoWriter('oscillation of simple pendulum.avi','uncompressed AVI');
open(videofile);
writeVideo(videofile,A);
close(videofile);

Function Program :-

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

In the above program, we have created a user-defined function for second order ODE of pendulum oscillation;
The name of the file is same as the name of funtion we will create. 

The initial value of the displacement & Angular velcoity is stored in variable theta1 & theta2 respectively.

The function that we have created is 'ode_function that will take values in the order (t,theta.b,g,l,m)

With the help of equation 1 and equation 2, we have wrriten the equation for 'dtheta1_dt' and 'dtheta2_dt'.

In equation 3, we have rearranged equation 1 and 2 in matrix form. Here too, we will do the same of assigning the equation i
'dtheta_dt'. As the function equations are written and they are assigned to variable, we will close the function using 'end' command.

Explanation :-

1. Initially, we have provided the input values of length of string(l), mass(m), gravity(g) and damping coefficient(b).

2. Thereafter, we have provided with the initial condition for anglular displacement and angular velocity as 0 radians and 3 rad/sec respectively.

3. Further, we have used linspace command to create 500 values for time span ranging from 0 to 20 sec and store this value in the array t_span.

4. Thereafter, we have tried to solve the 2nd order Ordinary differential equation using matlab's inbuilt ode45 function.

By using the ode45 function, it provides two outputs in the name of 't' & 'results'.

Syntax of ode45 is as under;

[t, results] =ode45(@(t,theta) ode_function(t,theta,b,g,l,m),t_span,theta_0)

Here we will call our function 'ode_function' that will take our given values and return the answer in variable [t,result].

As our equation is a function of time and theta, we will mention 't_span' and 'theta_0' in the end. Because, it has the values needed to be input and time limits needed for calculation. This is a variable which we use to showcase out 'time' and 'theta' function.

Once the array of [t,results] is created, we will plot graph of 'Displacement v/s Time' and 'Angular Velocity v/s Time'.

For 'Displacement v/s Time' we have written plot(t,results(:,1)), which means it will take all the answer from first column and plot Displacement versus time.

For'Angular Velocity v/s Time' we have written plot(t,results(:,2)), which means it will take all the answer from
second column and plot Angular Velocity versus time. For better illustration, I have set colour for both curve different as red and green respectively. Hold on command is used to have a complete curve in one frame.

Now to create the animation, we need to have basic idea of trignometry;

Fig:Position for (x1,y1) co-ordinates.

We can find the value of bob at different location of 'theta1' that is our angular displacement. Thats the reason, we will use our 1st column of results, as it has our result of angular displacement.
From above figure r = L;
x = x1;
y = y1;
and the pendulum is hung at origin (x0,y0) =(0,0).

Here, we will use a for loop for 'i' starting from one to length of our array results storing displacement values. Basically this, will iterate for 500 variables.

The hanging position x0,yo does not change. But for x1 and y1, our displacement value will change at every iteration and so will our x1 and y1 position.

The equations are;
x1=l*sin(results(i,1));
y1=-l*sin(results(i,1));
For every iteration, 'i' will change but column will remain same. This brings about a change in x1 and y1 value during every iteration.
Note that the sign before cosine element is negative, because we want our simulation taking place down, not other way round.

Using the code 'plot([-1,1],[0,0])' we will plot the hanging line and mention Linewidth as 8 mm and colour of line as blue.

Further, we have used line command as under to locate the xo,y0,x1 & y1 co-ordinates to visualize as a string.

line([x0 x1],[y0 y1],'Linewidth',4,'color','k');

Using the code 'plot(x1,y1,'-o','markers',20,'markerfacecolor','m')' we have actually made our pendulum bob
using marker and have set some specification to marker. The position of pendulum will change with string, as a result we have plotted it as x1 and y1.

Axis function {axis([-1.5 1.5 -1.5 0.5])} is used show that we can get our actual co-ordinates necessary to plot the figures and ultimately make animation.

The pause function actually pauses the frame for the specified amount of time. 'pause(0.01)' takes a pause of 0.01 sec in ever frame and it ultimately reflects in our video.

Here, once all frame is captured, we will make a video out of it using movie function. Movie will take the input as frames stored in A. We will make a variable videofile, where we use function VideoWriter. Here, we store the name of video in 'avi' format in uncompressed version. Then, we actually start video making using open function, where the name of the video that will be created will be 'oscillation of simple pendulum'. Then we use writeVideo function to create the video of the name 'oscillation of simple pendulum' in uncompressed version with the frames stored in A. Once that is done, we will use close function of videofile.

 

 

Ouput Plot of Angular Velocity vs Angular Displacement & Time :-

 

Animation link :- 

https://youtu.be/93XmHsw4_HM