Week 3 - Solving second order ODEs

AIM

Write a program in Matlab\Octave that will simulate the pendulum motion by solving the given second-order differential equation.

OBJECTIVE

• The objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

• `(d^2θ) /(d t) + b/m (d θ)/( d t) + g/L sin θ = 0`

• 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 the ordinary differential equation.

• To plot the displacement and velocity damping with respect to time.

• To make the animation of the simple pendulum.

DESCRIPTION & GOVERNING EQUATIONS

The differential equation is mostly used to define the physical phenomena into a set of mathematics equations. Differential equations are classified into ordinary differential equations and partial differential equations. Solving a 1^st order differential is quite easy and requires 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 equations and some structural phenomena are defined in the 2^nd order differential equation which needs to solve eventually to simulate 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 the periodic motion. By applying Newton’s second law motion we can obtain the equation of motion for a simple pendulum which is given as:

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

To solve the equation in MATLAB by using ‘ode45’ the equation needs to be converted from 2^nd order partial differential equation to 1^st order differential equation. Steps for converting 2nd order differential equation to first order equations are explained below

consider `θ=θ1`

The new equation of motion can be written as:

`(d^2θ1) /(d t) + b/m (d θ1)/( d t) + g/L sin θ1 = 0`

Consider,

`(dθ1)/ (dt) = θ2`

Therefore, the new equation can be written as:

`(dθ2)/(dt)+(b/m)θ2+g/Lsinθ1=0`

This leads to:

`(dθ2)/(dt) =−b/mθ2−g/L sin θ 1`

The resultant first-order equation will be:

`θ ′ 1 = θ 2 and θ ′ 2 = − b/m θ 2 − g/L sin θ 1`

MATLAB MAIN PROGRAM

%MATLAB PROGRAM FOR SIMPLE PENDULUM ANIMATION 
clear all
close all
clc

%% INPUT VALUES
l=1;     % Length of pendulum in m
m=1;     % Mass of the ball in kg
b=0.05;  % Damping coefficent
g=9.81;  % Acceleration due to gravity in m/s^2

% Time duration for simulation
t_span=linspace(0,20,100);  % time varies from 0 to 20 seconds and 100 values in betweeen

%% Initial Conditions, 
% in this solving two ODE, one for angluar Displacemnt and other is for angular velocity

theta_0=[0;3];  % initial angluar displacement of the pendulum=0 and initial angluar velocity=3rad/sec at time t=0

%% CALCULATION USING ode45 solver
%We are going to integrating the ODE for getting the solution

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

%% PLOT OF ANGULAR DISPLACEMENT & ANGULAR VELOCITY WITH RESPCT TO TIME 
figure(1) 
grid on
subplot(1,2,1)
plot(t,results(:,1),'color','r');
title('ANGULAR DISPLACEMENT vs TIME');
xlabel('TIME(s)');
ylabel('ANGULAR DISPLACEMENT(rad)');
subplot(1,2,2);
plot(t,results(:,2));
title('ANGULAR VELOCITY vs TIME');
xlabel('TIME(s)')
ylabel('ANGULAR DISPLACEMENT(rad)');
grid off
saveas(figure(1),'subplot','png')

%%POLAR COORDINATES TO CARTESIAN CORORDINATES

x1=l*sin(results(:,1));  %point of ball with respect to time
y1=-1*cos(results(:,1)); %point of ball with respect to time

%% VIDEO FOR PENDULUM MOTION

for i=1:length(t_span)
    figure(2)  % For animation 
    
    plot([0 x1(i)],[0 y1(i)],'linewidth',3,'color','y')
    hold on
    plot(x1(i),y1(i),'marker','o','linewidth',10,'markerfacecolor','g');
    plot([-1.2 1.2],[0,0],'linewidth',3,'color','g');
    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('simplependulum_4.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M)
close(videofile)
   

MATLAB FUNCTION PROGRAM

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

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 have created a for loop for which we are creating the frame motion for the simple pendulum and saved in the movie

RESULTS

GRAPH

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 decreased.

ANIMATION-YOUTUBE LINK

SIMPLE PEDNULUM ANIMATION

CONCLUSION

From the graph and from the animation we can see that we have clearly simulate the equation of motion for simple pendulum.