Study of Ordinary Differential Equation in matlab.

Aim: To study the Ordinary Differential Equation for pendulum in matlab.

Formulae:

Solving second order ODE:

Following is the procedure for solving second order ODE.

`frac{d^2theta}{dt^2}+frac{b}{m}frac{d theta}{dt}+frac{g}{L}sintheta=00`

Equation (01)

`theta=theta1`

Therefore,

`frac{d theta}{dt}=frac{d theta1}{dt}=theta2`

Equation (a)

`frac{d^2theta}{dt^2}=frac{d}{dt}(frac{d theta}{dt})=frac{d}{dt}theta2`

Equation (b)

Now, substituting (a) and (b) in (01)

Therefore,

`frac{d}{dt}(theta2)+frac{b}{m}(theta2)+frac{g}{L}sintheta1=00`

`frac{d theta2}{dt}=-frac{b}{m}theta2-frac{g}{L}sintheta1`

Objective:

For study the Ordinary Differential Equation of the pendulum.

To draw:

1. Displacement vs. Time plot.
2. Velocity vs. Time plot.
3. Create the animation of pendulum in matlab.

01.Introduction:

What is Ordinary Differential Equation?

An ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions.

An ODE of order n is an equation of the form:

Where,

y is a function of x, y′.

`frac{dy}{dx}`is the first derivatiive with respect to x

`frac{d^ny}{dx^n}`is the nth derivative with respect to x.

Order of Differential Equation:

The order of the differential equation is the order of the highest order derivative present in the equation. Here some of the examples for different orders of the differential equation are given.

• `frac{dy}{dx}=3x+2`. The order of the equation is 1
• `(frac{d^2y}{dx^2})+y=kt`. The order is 2
• `(frac{dy}{dt})+1=00`. The order is 1.

First Order Differential Equation:

You can see in the first example, it is a first order differntial equation which has degree equal to 1. All the linear equations in the form of derivatives are in the first order. It has only the first derivative such as dy/dx, where x and y are the two variables and is represented as:

frac{dy}{dt}=f(x,y)=y'`) = y’

Second-Order Differential Equation:

The equation which includes second order derivative is the second-order differential equation. It is represented as;

`frac{d}{dx}(frac{dy}{dx})=frac{d^2y}{dx^2}=f''(x)=y''`

Degree of Differential Equation:

The degree of differential equation is the power of the highest order derivative, where the original equation is represented in the form of a polynomial equation in derivatives such as y’, y”, y”’, and so on.

Suppose (d2y/dx2) + 2 (dy/dx) + y = 0 is a differential equation, so the degree of this equation here is 1. See some more examples here:

• `frac{dy}{dx}+1=00`, degree is 1
• `(y''')^3+3y''+6y'-12=00`, degree is 3

Applicatons:

Let us see some differential equations applications in real-time.

1) Differential equations describe various exponential growths and decays.

2) They are also used to describe the change in return on investment over time.

3) They are used in the field of medical science for modelling cancer growth or the spread of disease in the body.

4) Movement of electricity can also be described with the help of it.

5) They help economists in finding optimum investment strategies.

02. Problem Statement:

Solve the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

`frac{d^2theta}{dt^2}+frac{b}{m}frac{d theta}{dt}+frac{g}{L}sintheta=00`

In the above equation,

g = gravity in m/s^2,

L = Length of pendulum in m,

m = mass of the ball in kg,

b = damping coefficient,

Write a program in matlab that will simulate the pendulum motion.

Use,

L = 1 m,

m = 1 kg,

b = 0.5,

g = 9.81 m/s^2.

Simulate the motion between 0.20 seconds for angular displacement = 0,

angular velocity = 3 rad/sec at time = 0.

Animate the motion of the pendulum.

Solution:
Now we will write the program for the motion of the pendulum.
01.	% Study of the ODE for the motion of the pendulum. 
02.	 
03.	clear all
04.	close all
05.	clc
06.	 
07.	% inputs
08.	b = 0.5;
09.	g = 9.81;
10.	l = 1;
11.	m = 1;
12.	 
13.	% Initial condition
14.	% Here we have two ODEs Displacement and Velocity.
15.	theta_T = [0;3];
16.	 
17.	% time points
18.	t_span = linspace(0,20,200);
19.	 
20.	% Solving ODE. To solve ODE we need to call a function.
21.	[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_T);
22.	 
23.	plot(t,results(:,1))
24.	hold on
25.	plot(t,results(:,2))
26.	 
27.	ylabel('plot')
28.	xlabel('time')
29.	 
30.	ct = 1;
31.	for i = 1:length(t_span)
32.	    
33.	    THETA(i) = results(i,2);
34.	    
35.	    X0 = 0;
36.	    Y0 = 0;
37.	    
38.	    X1 = l*sind(THETA(i));
39.	    Y1 = -l*cosd(THETA(i));
40.	    figure(2)
41.	    
42.	    plot([X0 X1],[Y0 Y1],'linewidth',2,'color','g');
43.	    hold on
44.	    plot(X1,Y1,'marker','o','markersize',11,'markerfacecolor','r')
45.	    title('pendulum animation');
46.	    axis([-0.1 0.1 -1.2 0])
47.	    hold off
48.	    
49.	    M(ct)=getframe(gcf);
50.	    ct = ct + 1;
51.	end
52.	 
53.	        movie(M)
54.	        videofile=VideoWriter('C_pendulum_ani.avi','uncompressed avi')
55.	        open(videofile)
56.	        writeVideo(videofile,M)
57.	        close(videofile)

Figure01: Displacement and Velocity vs. time plot for the pendulum.

Detailed explanation of the program:

1. Name of Program, “Study of the ODE for the motion of the pendulum”.
3. Clears your workspace.
4. Closes all the figures.
5. Clears Command window.
7. Inputs
8. Damping coefficient (b) = 0.5
9. Acceleration due to gravity (g): 9.81 m/s^2
10. Length of the pendulum (l): 1 m
11. Mass of the pendulum ball (m): 1 kg
13. Defining the initial conditions.
14. Here we have two ODEs first is velocity and second is displacement.
15. theta_T = [0;3];

[At initial condition at any imaginary point in a plane pendulum remains at one position before releasing it. So in this condition the value of displacement is considered as zero whereas for the velocity the value is taken as 3 m/s considering the gravitational force acting on it. Here variable ‘theta_T’ is the column vector which represents the values of displacement and the velocity in first and second column respectively. Since it’s a column vector the values of both displacement and velocity are separated by using semi-colon (;)]

17. Defining the time points.
18. t_span = linspace(0,20,200);

[The variable ‘t_span’ holds the value of the time span for which the pendulum will be in motion. Here we took the value as 20 seconds and we defined it by using the ‘linspace’ command in order to take the required number (200) of values to make the animation of the pictures taken at the respective intervals of time.]

20. Solving ODE, for that we will have to call a function ‘ode_func’.
21. [t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_T);

[Matlab has inbuilt solver for the ODEs, here we are using ‘ode45’ which solves the ODEs of the first order only by using Range Kutta’s method. But we have to solve the equation of second order to obtain the results therefore, at first we have to convert the second order ODE into first order ODE. For that we will have to make a separate function program (ode_func) to do that. Here the variable [t,results] defines the inbuilt solver (ode45) which solves for the given values of t, theta and as the function (ode_func) is programmed.]

23. plot(t,results(:,1))

[Plotting the results of the column 1 of elements of data stored in the variable (t,results).]

24. hold on
25. plot(t,results(:,2))

[Plotting the results of the column 2 of elements of data stored in the variable (t,results).]

27. Labeling Y-axis.
28. Labeling X-axis.
30. Counter (ct) = 1
32. for i = 1:length(t_span)

[Using ‘for loop’ with variable ‘i’ holding values from 1 to over a length of data stored in variable ‘t_span’. Which simply means commands which are going to be executed in for loop will get executed only 1 to over a length of data stored in variable ‘t_span’.]

34. THETA(i) = results(i,2);

[Defining variable THETA(i) which will hold the values of the elements of data from column 2 of variable ‘results’ for the value of (i,2).]

36. X0 = 0;
37. Y0 = 0;

[Defining coordinates (X0,Y0) of fix end of the pendulum.]

39. X1 = l*sind(THETA(i));
40. Y1 = -l*cosd(THETA(i));

[Defining coordinates (X1,Y1) of free end of the pendulum.]

41. figure(2)

[Labeling the figure as figure(2).]

42. plot([X0 X1],[Y0 Y1],'linewidth',2,'color','g');

[Plotting pendulum by using ‘plot’ command and defining the coordinates of the pendulum.]

43. hold on
44. plot(X1,Y1,'marker','o','markersize',11,'markerfacecolor','r')

[Defining ‘marker’ command to draw the ball of the pendulum for animation, at free end.]

45. title('pendulum animation');

[Title of the animation, “pendulum animation”.]

46. axis([-0.1 0.1 -1.2 0])

[Defining the axis of the animation plot.]

47. hold off
49. M(ct)=getframe(gcf);

[Defining the variable ‘M’ which will hold the frame for each time the for loop is executed to create the animation of all frames altogether.]

50. ct = ct + 1;

[Defining the value of the counter for the next successive iteration.]

51. end

[End of for loop.]

53. movie(M)
54. videofile=VideoWriter('C_pendulum_ani.avi','uncompressed avi')
55. open(videofile)
56. writeVideo(videofile,M)
57. close(videofile)

[line (53-57) represents the format for creating the movie of the frames obtained after the execution of the for loop.]

Following errors were occurred during programming:

1. Error:

Wrong input coordinates for the axis command of the animation.

Figure02: Wrong values were entered for the coordinates of the axis command.

The values of the coordinates of the axis command must be entered from minimum to maximum for the array of X and Y.

Correction:

The values of the coordinates of the axis command were entered from minimum to maximum for both the axis of X and Y.

1. Error:

Wrong syntax of the plot command was used.

Figure03: Error in plot command.

Correction:

The syntax for plotting consists of arrays ‘X’ first and then arrays of ‘Y’.

Syntax: plot([x0 x1],[y0 y1],’linewidth’,value,’color’,’r’);

We used values of X0 and Y0 coordinates at first and the values of X1 and Y1 coordinates. Therefore the values were corrected as per syntax and arrays of the ‘x’ were entered at first and then the arrays of ‘y’ were entered.

Conclusion:

For a given problem we studied the programming of the pendulum.

As shown in figure01, the red curve represents the angular velocity of the pendulum for the given period of time.

And blue curve represents the angular displacement of the pendulum for the given period of time.

Click the link below to see, The behaviour (motion) of the pendulum for a given parameters.

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