Week 3 - Solving second order ODEs

AIM-To create a code in MATLAB,to solve the differential equation of motion of a pendulum, and to plot the displacement against time on a graph. The plotted images are then combined to create an animation of the pendulum motion.

EQUATIONS USED- The given second order differential equation represents the equation of motion of the pendulum. This second order equation is broken down into 2 separate first order differential equations. This is done by using simple mathematical relations. These 2 equations are then solved to obtain the values of theta(position) and time.

`(d^2theta)/dt^2``+b/m(d(theta))/dt+g/lsin(theta)`=0

where,

b=damping coefficient

g=acceleration due to gravity (`m/s^2`)

l=length of the pendulum (m)

m=mass of the ball (kg)

We consider 2 variables,`theta_1 and theta_2`. We derive the following 2 equations, by using simple mathematical relations. -

1.`(d(theta_1))/dt=theta_2`

2.`(d(theta_2))/dt=-(b/m)*theta_2-(g/l)*sin(theta_1)`

Thus, we have split the 2nd Order differential equation, into 2 separate equations.

WORKFLOW/STEPWISE EXPLANATION- We first write the obtained differential equations,in the form of a separate function,which can be called as required.The 2 1st order equations are stored as variables and then returned as the elements of 'dtheta_dt' array.This is what is done in the 'ode_func' function.

The inputs to this function are the constant parameters in the equation,which are b,g,l,m,whose values are entered by the user as required. Every differential equation requires an intial condition for the variables involved, in order to get the right solution. The initial conditions for displacement and velocity of the pendulum are stores in the 'theta_0' vector. Since there are 2 equations to be solved,there must be 2 initial conditions,one for angular displacement,and one for angular velocity.The required time points for integration are generated by the linspace() function.

The equations are solved using the 'ode45' function in MATLAB.As part of the syntax, we first have to mention the variables as an anonymous function, or function handle.This is done by using the '@' symbol,to indicate anonymous function.First the independent variable is indicated(here,time 't'),followed by the dependent variable(angular displacement,'theta'). After this, the function containing the 2 equations('ode_func'),is passed as an argument to ode45,followed by the time points('time'),and initial conditions vector('theta_0'). The output of ode45 is stored in 2 variables 't' and 'results'. ode45 always returns the output as a column vector. Here, the first column(results(:,1)) ,contains the first result,angular displacement, and the second column (results(:,2)),contains the second result,angular velocity.

The obtained values of angular displacement and velocity are plotted against time, and the graph is displayed.

In order to create the pendulum animation, we need all the displacement points at different instances of time.The first column of results,results(:,1) is extracted and stored in the vector 'pos'.A for loop is run through the length of 'pos'. The 'x' and 'y' position of the pendulum is obtained by taking the sin and cos component of the length of the pendulum.The initial position of the pendulum is assumed as 0,0.The line with points (0,0) and (xpos,ypos) (position of the pendulum) is plotted.A marker is used to represent the 'bob' or mass of the pendulum.The figure handle is then passed to the getframe() function and stored in an array 'M',using a counter variable.This is repeated for all the displacement values,i.e,for the length of 'pos' vector.

A MATLAB movie is generated using the movie() function. A videoWriter object is created and stored in videofile. It is of type,'Uncompressed AVI'. Into this object, we write the MATLAB movie, to create a '.avi' file, and save the animation.

MAIN PROGRAM-

%%Program to solve 2nd Order ODE of a pendulum
clear all 
close all 
clc 

%Constant Inputs 
b=0.05;
g=9.81;
l=1;
m=1;

%Initial Condition 
theta_0=[0;3];%Initial Displacement is 0m and initial Velocity is 5m/s

%Time Points 
time=linspace(0,20,500);

%Solving the ODE
[t,results]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m),time,theta_0);

%Plotting the results(Displacement and Veloicty)
figure(1)
plot(t,results(:,1),'linewidth',2.5);
hold on 
plot(t,results(:,2),'linewidth',2.5);
xlabel('Time(s)');
ylabel('Angular Displacement and  Angular Velocity');
legend('Angular Displacement','Angular Velocity');
title('Angular Displacement and Angular Velocity v/s Time');
hold off

%Animation of pendulum motion
pos=results(:,1);
figure(2)
count=1;
for i=1:length(pos)
    xpos=l*sin(pos(count));
    ypos=-l*cos(pos(count));
    plot([0,xpos],[0,ypos],'linewidth',3,'color','b');
    hold on
    plot(xpos,ypos,'marker','o','markersize',10,'markerfacecolor','r');
    axis([-2 2 -2 2]);
    xlabel('Displacement(x)');
    ylabel('Displacement(y)');
    title('Pendulum Animation');
    hold off
    M(count)=getframe(gcf);
    count=count+1;
end
movie(M);
videofile=VideoWriter('Pendulum2.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);

ODE Function-

%%Function representing the system of ODE's to be solved
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

OUTPUT-

The obtained plot of Angular Displacement and Angular Velocity, as a function of Time-

Link for the Animation -

https://drive.google.com/file/d/1JS8dCuFX77k7aHgeeQYpB_0OsPQtqqmr/view?usp=sharing

CONCLUSION- Thus,using the ode solvers such ode45() in MATLAB, different types of differential equations can be sovled, and their results can be analysed to study the behaviour and motion of different mechanical and physical systems.