Week 3 - Solving second order ODEs

AIM : To write a code in matlab in order to solve an ODE that represents the equation of motion of a simple pendulum with damping.

THEORY : Second order ODE involved:

Now,we have to convert this second order ODE to 2 first order ODE'S. For that, weproceed as follows:

let θ=θ1

dθ/dt=dθ1/dt=θ2 (where θ is the angular displacement and θ2 is the angular velocity.)

Now on further substitution,we get: dθ2/dt= -(b/m+g*sinθ1/L).

To carry out all the calculations and computations, we use the command :'ode45'

Ode45 is a function in matlab that solves ordinary or partial differential equations using a numerical method.

Syntax of ode 45 is:

[t,x]=ode45(@f(y),t,y0). (t and y0 are the initial conditions)

L=length of string=1 metre

m= mass of bob=1 kg

b= damping coefficient=0.05

g= acceleration due to gravity=9.81 m/s2

CODE:

1.Function Code:

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

2.Main code:

clear all
close all 
clc

%input arguments
b=0.05;
g=9.81;
l=1;
m=1;

%initial condition
theta_0=[0;5];

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

%calling the function
[t,results]= ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

figure(1)
%plotting the graph for motion of the pendulum
plot(t,results)
xlabel('Time')
ylabel('Displacement,Velocity')
grid on


ct=1;
x=results(:,1);
%plotting the moving string
for i=1:length(x)
  x0=0;
  y0=0;
  x1=l*sin(x(i))
  y1=-l*cos(x(i))
  
%plotting the bob  
figure(2)
plot([-1,1],[0,0],'linewidth',3)
hold on
plot([x0 x1] ,[y0 y1])
hold on
plot(x1,y1,'o','MarkerSize',30)
axis([-1.5 1.5 -1.5 1.5])
hold off
pause(0.01)
M(ct)=getframe(gcf);
ct=ct+1;

end

%moving pendulum movie

movie(M)
videofile =  VideoWriter('moving_pendulum','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

PROCEDURE :

1.Before starting the main code, we start off by defining the function code or the function that we have to call in the main code, i.e. the ODE. Here, theta1 is the angular displacement of the pendulum and theta2 is the angular velocity. The result is provided as displacement as rows and velocity in columns.

2. Now, starting with the main code and give the input arguments such as 'b','g','l' and 'm'. Then the initial conditions of the displacement are defined.

3.Then, we define the time interval for which the readings are being taken, here it is for 20 seconds with 500 intervals.

3. After that , we use the syntax: [t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0)

to calculate the corresponding values of displacement and velocity at different time points.

4.After that, we plot all the outcomes of the given commands on a graph as well as in a form of an animation which is saved in the form of an 'Uncompressed AVI' .

RESULT : The graph for the variation of displacement and velocity with time:

The link to the animation video for the moving pendulum with damping is: https://youtu.be/7LKtcD8U-_M