Week 3 - Solving second order ODEs

Our objective here was to write a MATLAB program to solve a given ODE(ordinary differential equation) that depicts the motion a simple pendulum with damping.

The given differential equation is a second order equation that depicts the motion of a simple pendulum experiencing a damping force.Apart from the we also have to animate this motion with the help of Matlab and create an animation video for the same.

Coming to the theory of the simple pendulum.

Simple pendulum is an ideal pendulum consisting of a point mass suspended by a weightless inextensible perfectly flexible thread and free to vibrate without friction —distinguished from physical pendulum.It is another mechanical system that moves in an oscillatory motion.

As we see in the above figure that the endulum shown here has a bob(that has a mass) a massless string or a rod that is connected to the pivot point.Here there is another quantity that is theta that denotes the value of angular displacement in radians.Now in order to solve the given ODE using a Matlab program we need to provide some pre defined values as the inputs that are provided in the question itself which are,

Expressions,

b- Damping co efficient

g- gravity in m/s^2

m-mass of the bob in kg

l-length of the pendulam in m

**Code explanation: Now moving on to the code explanation ,as mentioned earlier we take the inputs first that also include the initial values of angular displacement and angular velocities in the form of two column arrays.Where array 1 contains the angular displacement and arry 2 contains the angular velocity values.The time span is also provided for deciding the time interval so as to specify from what initial to final time the solutions are found.

l=1;                       %Length of the pendulum in metres
m=1;                       %Mass of the bob in kg
b=0.05;                    %Damping factor
g=9.81;                    %Acceleration due to gravity in m/sec^2

%Initial conditions for angular velocity and angular displacement
theta_0=[0;3];             %Notice semicolon given to store values as column matrix

%Initializing time span to find out nature of motion 
t_span=linspace(0,20,200);

2.Now we come to the step of solving the ODE using a predefined solver in Matlab "ode45(t,y(@....))". The ode45 solves a higher order differential equation by dividing it into further single order differential equation.In this case as the equation involved here is a second order ode the solver converts it into 2 single order differential equation.It solves the ode using 4th order runge kutta method.And the second order ODE is converted into 2 single order equations as shown below.

%%Using solver ODE45 to solve the second order differential equation of the
%%simple pendulum
[tsol,thetasol]=ode45(@(t,theta) ode_func(t,theta,b,m,g,l),t_span,theta_0);

Here if we take a closer look at what is happening when compiler executes ode465 is it first sets the parameters on which the equation is solved i.e. theta(dependent variable)and time(the independent variable) then it moves on to the function calling that goes to the funtion ode_func where the value of angular displacement and velocity are respectively calculated and returns the value dtheta 1_dt and dtheta2_dt to the solver the next thing is to decide till where the solution is desired that is done by t_span and the initial conditions i.e theta_0.The code for function ode_func is given below.

function [dtheta_dt]=ode_func(t,theta,b,m,g,l);

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

The values that are returned from the function that are the values of angular displacement and angular velocity which are stored in the array thetasol.

3.After that comes the plotting of the graph that depicts the solutions of the equation which is simply done by using plot command.

%Plotting graph for angular displacement and velocity in case of damping
%wrt time
figure(1)
plot(t_span,thetasol(:,1)) %Graph of angular displacement values(in column 1 of array theta sol)
hold on
plot(t_span,thetasol(:,2)) %Graph of angular velocity values(in column 2 of array theta sol)
xlabel('Time');
ylabel('Displacement and Angular Velocity');
legend('Displacement and Angular velocity wrt time');

4.Now in order to make an animation video we need to know the values of the angular displacement at each point of time but the values of the angular displacement are in radians which we need to change into rectangular co-ordinate system in order for plotting in x-y axes.We run a for loop for the above process and plot the respective x1 y1 values keeping the value of x0=0 and y0=0 that act as the fixed pivot point for the pendulum we depict the motion of the pendulum in the second third quadrant of the axes thus taking the y1 values as negative in order to fix the axes.

%%Now to depict the motion of a pendulum as an animation
c=1;                        %Initialized a counter variable used for frame capture
for i= 1:length(thetasol);  %Using a for loop to find out values for x1,y1(co-ordinates of the pendulum bob
    x0=0;                   %Fixing the pivot point at 0,0 on the axes
    y0=0; 
    x1=l*sin(thetasol(i));  %As the values obtained from theta were in radian (polar co-rodinate system
    y1=-l*cos(thetasol(i)); %We converted them into rectangular co-ordinate system

 %Fixing a line depicting the pivot point   
figure(2)
plot([-0.5 0.5],[0 0],'linewidth',2);%Pivot point line
grid on

hold on
plot([x0 x1],[y0 y1]);
plot([x0 x1],[y0 y1],'Color','r','linewidth',2);               %Ploting the line using values of x1,y1(here depicts the string joining the bob to the pivot point
plot(x1,y1,'Marker','o','markersize',20,'MarkerFaceColor','b');%Marking the pendulum bob using the circular marker
hold off
axis([-1.5 1.5 -1.5 0.5]);                                     %Fixing the values on the axes so as to depict the to n fro motion of the pendulum
title('Simple Pendulum')
pause(0.03);
M(c)=getframe(gcf);
c=c+1;                      %Counter increment with one
end  %End of loop

the marker command is used in order to indicate the bob and use a counter variable which will be helpful in frame capture when we create the animation.

5.Code for the animation video.

%Movie animation depicting the pendulum motion
movie(M)
videofile=VideoWriter('Simple Pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

So this was how we used matlab to solve the given ODE and created an animation of the motion characteristics along with the graph plots.

All the respective graphs are attached below along with the link to the animation video.Motion of a simple pendulum

The complete code for the program.

clear all;
close all;
clc
l=1;                       %Length of the pendulum in metres
m=1;                       %Mass of the bob in kg
b=0.05;                    %Damping factor
g=9.81;                    %Acceleration due to gravity in m/sec^2

%Initial conditions for angular velocity and angular displacement
theta_0=[0;3];             %Notice semicolon given to store values as column matrix

%Initializing time span to find out nature of motion 
t_span=linspace(0,20,200);

%%Using solver ODE45 to solve the second order differential equation of the
%%simple pendulum
[tsol,thetasol]=ode45(@(t,theta) ode_func(t,theta,b,m,g,l),t_span,theta_0);

%Plotting graph for angular displacement and velocity in case of damping
%wrt time
figure(1)
plot(t_span,thetasol(:,1)) %Graph of angular displacement values(in column 1 of array theta sol)
hold on
plot(t_span,thetasol(:,2)) %Graph of angular velocity values(in column 2 of array theta sol)
xlabel('Time');
ylabel('Displacement and Angular Velocity');
legend('Displacement and Angular velocity wrt time');
%%Now to depict the motion of a pendulum as an animation
c=1;                        %Initialized a counter variable used for frame capture
for i= 1:length(thetasol);  %Using a for loop to find out values for x1,y1(co-ordinates of the pendulum bob
    x0=0;                   %Fixing the pivot point at 0,0 on the axes
    y0=0; 
    x1=l*sin(thetasol(i));  %As the values obtained from theta were in radian (polar co-rodinate system
    y1=-l*cos(thetasol(i)); %We converted them into rectangular co-ordinate system

 %Fixing a line depicting the pivot point   
figure(2)
plot([-0.5 0.5],[0 0],'linewidth',2);%Pivot point line
grid on

hold on
plot([x0 x1],[y0 y1]);
plot([x0 x1],[y0 y1],'Color','r','linewidth',2);               %Ploting the line using values of x1,y1(here depicts the string joining the bob to the pivot point
plot(x1,y1,'Marker','o','markersize',20,'MarkerFaceColor','b');%Marking the pendulum bob using the circular marker
hold off
axis([-1.5 1.5 -1.5 0.5]);                                     %Fixing the values on the axes so as to depict the to n fro motion of the pendulum
title('Simple Pendulum')
pause(0.03);
M(c)=getframe(gcf);
c=c+1;                      %Counter increment with one
end  %End of loop

%Movie animation depicting the pendulum motion
movie(M)
videofile=VideoWriter('Simple Pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)