Week 3 - Solving second order ODEs

REPORT ON MATLAB PROGRAM TO SIMULATE THE PENDULUM MOTION

AIM =  To give the explanation on Matlab Program Code for creating an animation of the pendulum motion and on Errors encountered during the programming.

BODY OF THE CONTENT=

• The Matlab program to simulate the pendulum motion.

%Matlab Program to simulate the Pendulum Motion
clear all 
close all
clc

%Input parameters (Damping co-efficient,Gravity,Length of pendulum,Mass)
b = 0.05
g = 9.81
l = 1
m = 1

%initial conditions
theta_0 = [0;3];

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

% solve ODE
[t, results] = ode45(@(t,theta) ode_func_for_challenge(t, theta, b,g,l,m),t_span, theta_0)

ct = 1
for i = 1:length(results(:,1));
x0 = 0;
y0 = 0;
x1 = -1*sin(results(i,1));
y1 = -1*cos(results(i,1));

%plotting the line that will represent the string connected to fix end and
%Mass
plot([x0 x1],[y0 y1],'linewidth',2,'color','b');

hold on

%creating mass at the end the string
plot(x1,y1,'marker','o', 'markersize',10,'markerfacecolor','r');

grid on;
axis([-1.5 1.5 -1.5 2]);
pause(0.003);

 M(ct) = getframe(gcf);
 
 hold off
 
 ct=ct+1;
end
%creating movie of all the resulting frames
movie(M)
videofile =  VideoWriter('ode_2nd_order_challenge.avi','uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

• Explanation of the program

• Line 1

“%Matlab Program to simulate the Pendulum Motion”

This line shows the title of the program code. This shows for what task the program is been written. Here in this case it is to simulate the motion of the pendulum under the given parameters.

• Line 2 to 4

“clear all

close all

clc”

This lines deletes all the previous variables , close all the previous plots and clear the command window so that I can start a complete new program.

• Line 6 to 10

“%Input parameters (Damping co-efficient,Gravity,Length of pendulum,Mass)

b = 0.05

g = 9.81

l = 1

m = 1”

The input parameters that will be required in the program are given in this lines. The parameters includes Damping co-efficient, Acceleration of the gravity, Length of the Pendulum and Mass that is attached at the end of the string.

• Line 12 and 13

“%initial conditions

theta_0 = [0;3];”

This lines shows the initial conditions of the pendulum. Here in ‘[0;3]’ ‘0’ indicates the displacement at the initial state which is 0 in this case. ‘3’ indicates the initial angular velocity around 3 rad/sec.

• Line 15 and 16

 “ %time points

t_span = linspace(0,20,500)”

In these lines the values for time are been defined such as ‘t_span’ is an time array that contains values from 0 to 20 secs with 500 value in between. To achieve those values the ‘linspace’ is been used.

• Line 18 and 19

“% solve ODE

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

This lines are written to solve the ODE by using the an inbuilt solver in matlab ‘ode45’  and calling the function named ‘ode_func_for_challenge’. Here  ‘t’ and ‘theta’ are the inputs so we can say that ‘ode45’ is the function of ‘t’ and ‘theta’. After this the function is called having input values such as Damping co-efficient, Acceleration of the gravity, Length of the Pendulum and Mass that is attached at the end of the string. Then the time array is given in the code that shows how far the integration needs to be done. Then the initial conditions corresponding to the starting of time is added to the code. Here the output of this line is the integration time and results which shown as ‘[t,results]’.

• Line 21

“ct = 1”

In order to create the animation the ‘for loop’ will be used in the program for that the counter(ct) is added to the program with the start value of 1.

• Line 22

“for i = 1:length(results(:,1));”

In this line the ‘for loop’ command is used to apply all the values of the results in column 1 which are related to the displacement.

• Line 23 to 26

“x0 = 0;

y0 = 0;

x1 = -1*sin(results(i,1));

y1 = -1*cos(results(i,1));”

The coordinates of the string that will connect the fix point and the Mass are added to the code in these lines. The formula to calculate x1 and y1 are designed for radian values of the angle hence the ‘sin’ and ‘cos’ commands are used here instead of ‘sind’ and ‘cosd’ which consider angular values in degree.

• Line 28 to 30

“%plotting the line that will represent the string connected to fix end and

%Mass

plot([x0 x1],[y0 y1],'linewidth',2,'color','b');”

In this line the ‘plot’ command is been used to plot the line that will represent the string connected to the fix end and Mass with the linewidth value ‘2’ and color ‘blue’.

• Line 32

“hold on”

As there will be two plot command used in the program hence we need to give a ‘hold on’ command between those so that they would not overlap each other.

• Line 34 and 35

“%creating mass at the end the string

plot(x1,y1,'marker','o', 'markersize',10,'markerfacecolor','r');”

Here, the ‘plot’ command is again used to create circular marker at the end of the string that will represent the Mass. ‘markersize’ command is used to give the size value to the marker in this case it is ‘10’ and ‘markerfacecolor’ command is used to give color to the face of the maker. Here the red color is used.

• Line 37 to 39

“grid on;

axis([-1.5 1.5 -1.5 2]);

pause(0.003);”

The ‘grid on’ command is used to put the grid on so it becomes more easy to get nay value on the plot. In next line the x and y axis limits are given and after that the pause value is added to the code.

• Line 41

“M(ct) = getframe(gcf);”

Here, by using the ‘getframe’ command the individual frames are been generated in order to create the movie and those all frames are saved in an array ‘M(ct)’ where ‘ct’ is an individual frame. So in this way the windows that we get when we run the program are combined in the array to form the animation.

• Line 43

“hold off”

In this line the ‘hold off’ command is used so that after each and every iteration the new values will be displayed on the frame and there will not be any overlapping of the plot at each and every iteration.

• Line 45

“ct=ct+1;”

In this line the value of the counter is getting incremented by 1 which is been utilized by the ‘for loop’ whenever it is been executed.

• Line 46

“end”

This syntax indicates that the for loop has ended here.

• Lines 47 and 48

“%creating movie of all the resulting frames

movie(M)”

In this line the movie is been generated by the matlab by using the ‘movie’ command. Here all the frames that were generated in array ‘M’ at line no. 41 are combined to get the movie.

• Lines 49

“videofile =  VideoWriter('ode_2nd_order_challenge.avi','uncompressed AVI')”

Here by using the ‘VideoWriter’ command the new ‘.avi’ file is generated with the profile as ‘uncompressed AVI’. In this file all the frames will placed. It names must be as same as the name of the program.

• Line 50 to 52

“open(videofile)

writeVideo(videofile,M)

close(videofile)”

In those line all the captured frames are been written into the file that were generated on line 49 by opening the file and the after its completion the file is been closed.

• The explanation of the function code

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

• Line 1

“function [dtheta_dt] = ode_func_for_challenge(t,theta,b,g,l,m)”

In the first line itself the function command is been used to create the new function having inputs such as time, theta value, Damping co-efficient, Acceleration of the gravity, Length of the Pendulum and Mass that is attached at the end of the string. In this function the 2^nd  order differential equation is been used as shown below.

`frac{del^2theta}{del t^2} + frac{b}{m}frac{del theta}{del t} + frac{g}{L} sintheta = 0`

To solve this differential equation it is expressed as two single order differential equations as shown below.

`frac{del thetaunderset(1)()}{del t} = thetaunderset(2)`

`frac{del thetaunderset(2)()}{del t} = - frac{b}{m} thetaunderset(2) () - frac{g}{L} sinthetaunderset(1)`

• Line 3 and 4

“theta1 = theta(1)

theta2 = theta(2)”

Here, in these lines the values of theta1 and theta2 are given which will be used further in the code by ‘dtheta1_dt’ and  ‘dtheta2_dt’

• Line 5 and 6

“dtheta1_dt = theta2;

dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);”

Here, the values of two single order differential equations are been defined into the code.  

• Line 7

“dtheta_dt = [dtheta1_dt; dtheta2_dt];”

Here, both the ode’s are combined into an array called ‘dtheta_dt’. Both the ode’s are provided as two rows by giving the ‘;’ sign between them.  

• Line 8

“end”

In this line the ‘end’ syntax is been given which is necessary to give as it indicates the function has ended here.

• Error encountered during programming

Error = There was an error at line 35 showing ‘value not a numeric scalar’

Measures that were taken = I had given single inverted comma to the marker size value so I removed those single inverted commas and ran the program again.

• Youtube video link of the animation 

Video link = https://youtu.be/8UbwCNulDAQ

CONCLUSION = From the above entire program it can concluded that in order to create the animation it very important to follow  the right order of the command in the program to get required output.