SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

LAB#1:
INTRODUCTION TO MATLAB
(BASIC CONCEPTS)

The learning objectives are:

• Starting and Ending a MATLAB session.

• Learn about the major components of MATLAB environment.
• Learn to use Help Browser and help Command.
• Know about different types of files used in MATLAB.
• Know the various platforms where MATLAB can be used.
• Some other useful general commands.

Introduction:
MATLAB stands for MATRIX LABORATORY. A special purpose computer program optimized to
perform Engineering and Scientific calculations or it is a Technical programming language. It
is started with simple matrix manipulation and grown with the capability of solving any
technical problem. It is a proprietary programming language developed by Math-works. It
provides a very extensive library of pre-defined functions to make technical programming task
easier and more efficient. It is superior on other languages like C, C++, FORTRAN.

MATLAB Capabilities:
 Numerical computation
 Data Access
 Data Analysis and Visualization
 Programming and Algorithm development
 Application development and deployment

Data Analysis and Visualization:

 Preprocessing
 Fast and accurate analysis (Correlation, Fourier analysis and Filtering, Basic Statistics
and Curve fitting)
 Built-in functions for Matrix operations (Cholesky factorization, LU Factorization)
 Special Math functions (Airy function, Bessel function, Beta and Gamma function,
Laplace transforms, Fourier transforms, Interpolation, Polynomials and ODE solvers)
 Built-in Graphics for Engineering and Science (2D and 3D Plots, Volume Visualization)

MATLAB MANUAL 1
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

MATLAB Environment:

The major components of MATLAB environment are as follows:

 Command window
 Command history
 Workspace
 Current folder
 Figure window
 Edit window

1. Command Window:
Whenever MATLAB is involved the main window called command window is activated. The
command window displays the command prompt “>>” and a cursor where commands are
entered and are executed instantaneously on pressing the “enter” key of the keyboard.
Commands can be run by typing them in the command window.

2. Command History Window:

It consists of a list of all the commands that are entered at the command window. It consists
of commands of previous session also. These commands remain in the list untile they are
deleted. Any command may be executed by selecting and double clicking it with mouse. A
program file may be created by selecting a set of commands and right clicking the mouse. On
right clicking the mouse, pop-up menu is displayed. A program file containing the selected
commands can be created by choosing “Create [Link]” option from menu.
Similarly, commands may be deleted by selecting the command and right clicking the mouse
and selecting “Delete selection” from pop-up menu.

MATLAB MANUAL 2
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

3. Workspace:
A workspace is a collection of all the variables that have been generated so far in the current
MATLAB session and shows their data type and size. The workspace information can also be
obtained by typing command at the command prompt. The “who” and “whos” command will
generate list of variables currently in the workspace. Workspace stores generated variables
temporarily. The workspace will be cleared when we closed the MATLAB. All the commands
executed from command window and all the script files executed from the command window
share common workspace, so they can share all the variables. Using these variables, a
number of operations can be done such as plotting by selecting a variable and right clicking
the mouse and selecting the desired option from pop-up menu.

4. Current Folder:
It contained all executed programs and save MATLAB files. In the current folder window, all
the files and folder present in the current directory will be listed. To run any file it must either
be in the current directory or on the search path. A quick way to view or change the current
directory or its content is by using the current directory field in MATLAB toolbar. For this
click on “DESKTOP” on the menu bar and check the boxes to show window.

5. Edit Window:
An edit window is used to create a new program file or to modify existing files. In the window,
programs can be written, edited and save. The programs written using a MATALB editor are
automatically assigned an extension “.m” by the editor and are known as M. files. A new M.
file can be created by selecting an appropriate option from the desktop “file” menu. An
existing M. file can be opened by selecting “file” from menu bar and then click “Open” from
desktop menu. The window editor is open automatically when a new M. file is created or and
existing M. file is opened.

6. Figure Window:
A figure window is a separate window with default white background and is used to display
MATLAB graphics. The result of all the graphics command executed or displayed in the figure
window. There can be any number of figure windows depending upon the system memory. To
open figure window write “figure” in command window and press enter key. A figure window is
displayed.

MATLAB Symbols with their Operations:

Symbols Operations
+ Used for addition e.g. >>2+3
- Used for subtraction e.g. >>3-2
* Used for multiplication e.g. >>2*x
MATLAB MANUAL 3
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

/ Used for division e.g. >> 2/5

^ Used for exponentiation e.g. >>2^4 →24
sqrt To find square root e.g. sqrt(4)
% Used to write comments or any statements that we do not want to execute.
pi Used for π (3.141592……)
To terminate a line and suppress the result to display on command window
; e.g. >>x=2;
It will not be executed
Used to define a range
e.g.
we want to write counting from 1 to 100 with gap 1, then >>[Link]
It will display all numbers from 1 to 100.
: Generally,
Initial point : desired gap : final point
>>[Link] %TABLE OF 2
It will display.
2 4 6 8 …...20
It is typed in the beginning of a line in edit window. It create cell (highlighted
%%
rectangle) which can be executed separately / individually by using icon.

Some useful MATLAB Commands:

Commands Description
clc To clear screen of command window
clf To clear figure window
who List the variables currently in the workspace
whos Same as “who” command but gives more information such as type and size
what Display all saved file on command window
clear To clear workspace
clear all To clear variables in workspace
close all To close all windows
plot To plot a graph
fprintf To print the output on the screen / command window
quit/exit To close the MATLAB session

Docking and Undocking:

Docking means attaching the window to MATLAB environment. This is done by clicking the
icon on the top right corner of every window. Undocking means detaching the window from
MATLAB environment.

Abort:
MATLAB MANUAL 4
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

In order to abort a command in MATLAB hold down the control key and press “C” to generate
a local abort MATLAB.

LAB#2:
CONSTANTS, VARIABLES AND EXPRESSIONS

The learning objectives are:

 Knowledge of MATLAB character set and Data types.

 Study of constants and variables and related useful functions.
 Study of functions used for conversion of Co-ordinate systems.
 Study of arithmetic, relational and logical operators.
 Knowledge of hierarchy of operations and built-in functions.
 Use of assignment statement.

Introduction:
MATLAB MANUAL 5
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

MATLAB is a powerful computing language. Programming languages process data consisting

of numbers, characters and strings and provide output, known as useful information. Data
processing is carried out by executing a sequence of instructions called a computer program.
This group of instructions or program is written by the user, using certain symbols and words
according to syntax rules of a particular language. MATLAB also has its vocabulary and
grammar. In this Lab, concepts of character set, constants, variables, data types, operators
and expressions are discussed.

Character Set:
It is a set of characters that form a part of the statements written using the programming
language. Characters can be broadly classified into four categories:
 Alphabets
 Numerals
 Special characters
 White space characters
MATLAB recognizes all 26 alphabets of English language. It is case-sensitive that is, upper-
case letters and lower-case letters are recognized as different characters. Numerals 0 to 9 are
recognized in MATLAB special characters.
White space characters like tab, blank, new line, vertical tab and line-feed are also included in
the character set.

Data Types:
Data types are in form of array. This array is a minimum single element in size and can grow
to an n-dimensional array of any size. Two dimensional version of these arrays are called
Matrices.
1. Character Data Types:
Alphabets or words e.g. MATLAB is a character array. Array of strings (alphabets) can be
used to store number of strings, provided all of them have the same length that is 16 bits
long.

2. Numeric Data Types:

Single data type e.g. int8, int16, int32 and int64 as well as double data type e.g. uint8,
uint16, uint32, and uint64 includes signed and unsigned integer arrays that are 8, 16, 32,
64 bits in length respectively.

Constants and Variables:

Constants:
Constants refer to find values which donot change during the execution of a program.
Constants can be of two types:

1. Numeric Constants:
MATLAB MANUAL 6
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Numeric constants are of three types. That are Integer constants, Real constants and
Complex constants. MATLAB does not differentiate between Integer and Real constants
during execution.

Commands Related To Complex Numbers:

Let x is defined as: x=a+ bi
Sr. # Operations Description
1 real ( x) Give real part of x
2 imag( x ) Give imaginary part of x
3 conj( x ) Produce conjugate of x , i.e x=a−bi
4 |( x)| Use to obtain magnitude of complex number (|x|=√ a2 +b2 )
5 angle (x ) −1 y
Use to obtain angle of complex numbers in radians (θ=tan )
x

By default MATLAB accepts and return angle in radians.

Conversions of Coordinate System:

Conversion Command
[theta,r]=cart2pol(x,y)
Cartesian To Polar e.g. >>x=1; y=1;
>>[theta,r]=cart2pol(1,1)
( x , y ) →(r , θ)
>>theta=0.7854
r =1.4142
[x,y]=pol2cart(theta,r)
Polar To Cartesian e.g. >>theta=1, r=5;
(r , θ)→ ( x , y ) >>[x,y] = pol2cart(1,5)
>>x=2.7015 y=4.2074

2. Character Constant:
Character constants are of following types
 Single Character constants
 String Constants
 Escape Sequence Constants

Single Character Constants

A single character constant contains a single character enclosed with a pair of single quote
marks. Some examples are valid single character constants are ‘y’,’g’,’1’,’_’,’ ’. It may be noted
that space is also a character.

String Constants:

MATLAB MANUAL 7
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

A string constant is a sequence of characters enclosed in pair of single quotes. The characters
may be alphabets, numbers, special characters or blank space. e.g. ‘GOOD’ ‘1947’ ‘BEST
WISHES’. A word containing a single quote is specified by using two single quote together e.g.
isn’t written as ‘isn’t’.

Escape Sequence Constants:

Escape Sequence character constants are used in output functions. For example
\b Backspace
\t Tab
\n New line
\r Carriage return
\f Form feed
\\ Print back slash
\’ Print apostrophe
%% Print % or create a cell

Variables:
Variables and constants form an integral form of program in any language. The different types
of data are identified by their names for case in reference. Programming languages require a
programmer to declare these variables and the type of data in advance at the beginning of
program, whereas no such declaration is required in MATLAB.
 Blank spaces cannot be included in MATLAB variables names.
 MATLAB is case-sensitive. Lower-case and upper-case latters represent two different
variables names, for example “GOLD” and “gold” are different variable names.
 Words that are reserved for the key constructs of a programming language are called
“KEYWORDS” and cannot be used as variable names. For example: for, while, switch
and functions are keywords.

Special Constants & Variables:

Pi 3.141519
i/j √−1
Inf ∞
Ans Default output variables

Built-in Functions:
The functions which are already available in MATLAB are called built-in functions. A complete
list of all MATLAB functions can be seen through Help browser. General Syntax for using
these functions is:
Function_name (expression)

MATLAB MANUAL 8
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Built-in function Description

exp(x) Calculatee x
log(x) Calculate Natural logarithm
factorial(x) Calculate Factorial of a number
sin(x) Compute the sine of ‘x’ in radians
cos(x) Compute the cosine of ‘x’ in radians
tan(x) Compute the tangent of ‘x’ in radians
csc(x) Compute the cosecant of ‘x’ in radians
sec(x) Compute the secant of ‘x’ in radians
cot(x) Compute the cotangent of ‘x’ in radians
sind(x) Compute angles in degree
cosd(x)
tand(x)
cscd(x)
secd(x)
cotd(x)
asin(x) Compute the inverse sine of ‘x’ where ‘x’ must between -1 and 1. The

−π π
function returns an angle in radians between and
2 2
acos(x) Compute the inverse cosine of ‘x’ where ‘x’ must between -1 and 1. The
function returns an angle in radians between 0 and π
atan(x) Compute the inverse tangent of ‘x’. The function returns an angle in

−π π
radians between and
2 2
atan2(y,x) y
Compute the inverse tangent of value the function returns an angle
x
in radians that will be between –π and π depending on signs of ‘x’ & ‘y’.
sinh(x) e x −e− x
Compute the hyperbolic sine of ’x’ which is equal to
2
cosh(x) e x + e−x
Compute the hyperbolic cosine of ’x’ which is equal to
2
tanh(x) e −e− x
x
Compute the hyperbolic tangent of ’x’ which is equal to x − x
e +e
asinh(x) Compute the inverse hyperbolic sine of ’x’ which is equal to ln [x + √ x 2+1]
acosh(x) Compute the inverse hyperbolic cosine of ’x’ which is equal to

ln [x + √ x 2−1]
atanh(x) Compute the inverse hyperbolic tangent of ’x’ which is
1+ x
equal to ln
√ 1−x
for |x|≤1

Element by element operation:

MATLAB MANUAL 9
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

We use element by element operation, when we have expressions in fractional forms.

Operations Operators Algebraic form MATLAB form

Element by element .* p ×q p.*q
multiplication
Element by element ./ p/q p./q
division
Element by element .^ pq p.^q
exponentiation

ILLUSTRATIVE PROGRAMS:

Task#1:
Write a program to calculate area of circle having radius 3cm?

AreaOf ¿˚ π r 2

Program:
r = 3;
area = pi * (r ^ 2)

Output:
area =
28.26

Task#2:
Given two sides a=3.2 & b=4.6 of a triangle and angle 60o between these two sides, find the
length of the third side and the area of the triangle?

C=√ a2 +b 2−2 ab cos(θ)

1
AreaOf Triangle= ab sin ( θ )
2
Program:
a = 3.2;
b = 4.6;
theta = 60;
theta = theta * pi / 180;
c = sqrt ( a ^ 2 + b ^ 2 – 2 * a * b * cos (theta) )
area = 0.5 * a * b * sin (theta)

Output:

MATLAB MANUAL 10
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

c=
4.0841
area =
6.3739

Task#3:
Write a program to convert temperature given in oC say 35.4oC to oF?
o
F= 9/5 oC+32

Program:
degcenti=35.4;
degfah=1.8*degcenti+32

Output:

degfah =

95.7200

Task4#:
Write MATLAB statement and calculate sum of series for x=1.5?

x2 x4 x6 x8
s=1− + − +
2! 4 ! 6 ! 8 !

Program:
x=1.5;
x1=1;
x2=(x*x)/(factorial(2));
x3=(x^4)/(factorial(4));
x4=(x^6)/(factorial(6));
x5=(x^8)/(factorial(8));
S=x1-x2+x3-x4+x5

Output:
S=

0.0708

Task#5:
Evaluate the following assignment statements.

1- q=2 log 10 x +cos π +| y 2−z 2|+ √ 5 yz for x =2, y=4∧z=3

2 1
2- ln ( x + x+ 1 ) where x= ∧x=1
2

MATLAB MANUAL 11
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

x π π
3- where x= ∧x=
2
(x +1)(sin x ) 4 2

Program:
%%
x=2;
y=4;
z=3;
pi=3.14;
q=2*log(10)*x+cos(pi)+abs((x^2)-(z^2))+sqrt(5*y*z)

%%
x=1;
a=log((x^2)+x+1)

%%
x=1/2;
b=log((x^2)+x+1)

%%
x=180/4;
c=x/(((x^2)+1)*sin(x))

%%
x=180/2;
d=x/(((x^2)+1)*sin(x))

Output:

q=
20.9563
a=
1.0986
b=
0.5596
c=
0.0261
d=
0.0124

Task#6:
Evaluate the given number for ‘x’ from 1 to 2 in steps of 0.1.

1 x3
y= + 4
x x +5 x sin x

Program:

MATLAB MANUAL 12
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

x=1:0.1:2;
y=(1./x)+((x.^3))./((x.^4)+(5.*x.*sin(x)))

Output:
y=

Columns 1 through 10
1.1920 1.1182 1.0587 1.0102 0.9698 0.9357 0.9065 0.8810 0.8583 0.8378

Column 11
0.8188

Task#7:
Write assignment statements to evaluate the following equations:

(A)Volume=π r 2 h for r=2 , h=3

2 πNT
(B) Power= for N=1000 , T =10
60

VI cos θ
(C) Eff iciency= for V =250 , I =5 , R=2 , cos θ=0.8 , ωc =20
VI cos θ+ I 2 R+ ωc

Program:

%%
pi=3.14;
r=2;
h=3;
V=pi*(r^2)*h

%%
pi=3.14;
N=1000;
T=10;
P=(2*pi*N*T)/60

%%
V=250;
I=5;
cos(theta)=0.8;
R=2;
Wc=20;
E=(V*I*cos(theta))/((V*I*cos(theta))+(I^2)*R+Wc)

Output:

MATLAB MANUAL 13
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

V=
37.6800

P=
1.0467e+003

E=
71

Task#8:
For an electrical circuit with an inductance l=0.01mh and resistance r=180 ohms the damped
natural frequency of oscillations is:

1 R2
Frequency=
√ −
LC 4 C2

Write a program to calculate the frequency for different values of ‘c’ varying from 0.1 to 1 in
step of 0.1.

Program:
L=0.01;
R=180;
C=0.1:0.1:1;
F=sqrt((1./L.*C)+(R.^2)/4.*(C.^2))

Output:
F=
Columns 1 through 9
9.5394 18.5472 27.5500 36.5513 45.5522 54.5527 63.5531 72.5534 81.5537

Column 10
90.5539

MATLAB MANUAL 14
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

LAB#3:
VECTORS AND MATRICES IN MATLAB
The learning objectives are:

 Knowledge of scalar, vectors and matrices and basic operations such as product and
transpose
 Useful commands/functions related to vectors and matrices.
 Entering data to a matrix, line continuation, matrix indices, dimension statement.
 Study of different matrix operations.
 Knowledge of special matrices.
 Study of different commands/functions used with matrices .

Introduction:
The matrix notation usually simplifies the complex mathematical expressions /equations and
makes solution of problems much easier to handle and manipulate. In MATLAB, a matrix is a
rectangular array of real or complex numbers. Matrices with only one row or with only one
column are called row and column vectors, respectively. A matrix having only one element is
called a scalar.

Although other higher programming languages work with one number at a time, in MATLAB it
is possible to work with complete matrix simultaneously. This feature is very important as it
removes the un-necessary loops and repetition of same statements. The program, therefore,
becomes short, concise and easily understandable. In MATLAB, matrix is chosen as a basic
data element. All variables when used as a single data element are treated as single element
matrix that is a matrix with one row and one column.

Arrays:
An array is a list of numbers arranged in rows and/or columns. A one-dimensional array is a
row or a column of numbers and a two dimensional array has a set of numbers arranged in
rows and columns. An array operation is performed element by element.

Scalars and Vectors:

MATLAB MANUAL 15
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

A scalar is a 1 ×1 matrix containing a single element only. A column vector is an m× 1 matrix

that has m-number of rows but a single column only. A row vector is a 1 ×n matrix which has
n-number of columns and has only one row.

Assigning data to elements of a Scalar/Vector:

Row vector:
In a row vector, the elements are entered with a space or a comma between the elements
inside the square brackets e.g. x=[a ,b , c , d ] or x=[a b c d ]

Column vector:
In a column vector, the elements are entered with a semicolon between the elements inside
the square brackets e.g. x=[a; b ; c ; d ]

Scalar:
A scalar does not need any square brackets, for example a scalar x=10 can be entered as:
x=10 ;

Creation of evenly spaced Row-vectors:

A row vector with evenly spaced, real elements can be created by using the following
command:

Variable name=K : L: M

K is the initial value of the row vector

L is the step-size or increment
M is the final value of row vector
Care should be taken that the initial and final values must be appropriate. If a step-size is not
specified, +1 is taken as default value of the step-size.

Commands for Evenly Spaced Row-Vectors:

linspace:

A row vector can also be created by using linepsace command. The syntax is given as follows:

linspace ( x a , x b , n)

It generates a linearly space vector of length ‘n’ that is ‘n’ equally spaced points between x a
and x b. e.g. A=linspace(0,10,5)

When executed, produces a vector having 5 equally spaced elements between the limits 0 and
10 given as follows: A=[0 2.5 5.07.5 10.0]

MATLAB MANUAL 16
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Built-In Functions for Handling Arrays:

Some of the built-in functions available in MATLAB for managing and handling arrays are
listed in table below:

Function Description Example

>> A=[5 9 2 4];

Returns the number of elements in
length(A) >>length(A)
the Vector A.
>>ans = 4
>>size(A)
size(A)
Returns a row vector [m,n], where m >>ans = 1 4
or
and n are the size mxn of the array A (means 1 row and 4
[m,n]=size(A)
columns)
>>V=[3 2 1];
>>A=diag(V)
When ‘V’ is a vector, creates a square
>>A =
diag(V) matrix with the elements of ‘V’ in
3 0 0
diagonals.
0 2 0
0 0 1

It lists the elements of the vector ‘x’ in >>x=[5 -3 10 -10]

sort(x,’mole’)
ascending or descending order. >>a=sort(x)
sort(x,’ascend’)
If we only write sort(x), it will sort the >>a=[-10 -3 5 10]
sort(x,’descend’
array in ascending order >>sort(x,’descend’)
)
automatically. >>b=[10 5 -3 -10]

Built-In Functions for Analyzing Arrays:

Some of the many built-in functions available in MATLAB for analyzing arrays are listed
below:

Functions Description Example

mean(A) If ‘A’ is vector, returns the mean value of >>A= [3 7 2 16];
elements >>mean(A)
>>ans = 14
c=max(A) If ‘A’ is a vector, C is the largest element in >>A= [3 7 2 16 9 5
‘A’ 18 13 0 4];
>>c=max(A)
>>c=18
[d,n]=max(A) If ‘A’ is a vector ‘d’ is the largest element in >>[d,n]=max(A)
‘A’, ‘n’ is the position of the element (the >>d=18
first if several have the max value) n=7

MATLAB MANUAL 17
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

min(A) The same as max(A), but for the smallest >>A=[3 7 2 16];
[d,n]=min(A) element >>min(A)
The same as [d,n]= max(A) but for the >>ans=2
smallest element
sum(A) If ‘A’ is a vector, return the sum of the >>A=[3 7 2 16];
element of the vector >>sum(A)
>>ans=28
median(A) If ‘A’ is a vector, return the median value >>A=[3 7 2 16];
of the element of vector >>median(A)
>>ans=5
std(A) If ‘A’ is a vector return the standerd >>A=[3 7 2 16];
derivation of elements of vector >>std(A)
>>ans=6.3770
det(A) Return the determinant of a square matrix >>A=[3 7 2 16];
‘A’ >>det(A)
>>ans=-2
dot(A,B) Calculates the scalar (dot) product of two >>A=[5 6 7];
vectors A & B. The vector can each be row >B=[4 3 2];
or column vectors >>dot(a,b)
>>ans=52
cross(A,B) Calculates the cross product of two >>A=[5 6 7];
vectors A & B (AxB). The two vectors must >B=[4 3 2];
have three elements. >>cross(a,b)
>>ans=-9 18 -9

Matrix:
A matrix is a two dimensional array which has n-number of rows and columns. A matrix is
entered row-wise with consecutive elements of a row separated by a space or a comma, and
the columns are separated by semi-colons or carriage returns. The entire matrix is enclosed
with in square brackets. The elements of the matrix may be real numbers or complex
numbers. e.g.

To enter the matrix

A= 1 3 −4
[
0 −2 8 ]
In MATLAB input command is:

A=[1 3 −4 ; 0 −2 8]

Similarly for B= [−53ix ln 2 x +7 sin 3 y

5−13i ]
In MATLAB: B=[ −5∗x log ( 2∗x ) +7∗sin( 3+ y); 3∗i 5−(13∗i) ]

Matrix Subscripts/Indices:
MATLAB MANUAL 18
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

The elements of a matrix can be specified/ accessed by specifying their respective row and
column numbers. The first index refers to the row numbers and second index refers to
column number. The elements in ith row and jth column of a matrix A is denoted by A(i , j).
For example A ( 2,3 )refers to the element in the second row and the third column of the matrix
A.

Colon for a Matrix:

The table below gives the use of a colon in addressing arrays in a matrix.

Command Description
A(: ,n) Refers to the elements in all the rows of the column ‘n’ of the matrix
‘A’.
A(m, :) Refers to the elements in all the columns of the row ‘m’ of the matrix
‘A’.
A(: ,m :n) Refers to the element in all the rows between column ‘m’ and ‘n’ of
the matrix ’A’.
A(m:n , :) Refers to the elements in all the columns between row ‘m’ and ‘n’ of
matrix ‘A’
A(m:n , p : q) Refers to the elements in rows ‘m’ through ‘n’ and columns ‘p’
through ‘q’ of matrix ‘A’.

Element By Element Operations:

Element by element operations can only be done with arrays of the same size. Element by
element multiplication, division and exponentiation of two vectors or matrices is entered in
MATLAB by typing a period in front of arithmetic operator. Table below consist of list of these
operations.

ARITHMETIC OPERATORS
Matrix Operators Arrays Operators
+ Addition + Addition
- Subtractions - Subtractions
* Multiplication .* Array Multiplication
^ Exponentiation .^ Array Exponentiation
/ Left Division ./ Array Left Division
\ Right Division .\ Array Right Division

Some Useful Commands Related To Matrices:

Command Description Example

det(A) It returns the determinant of a matrix >>A=[ 1 2 ; 0 4]
(| A|). >>det(A)
>>ans=4

MATLAB MANUAL 19
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

rank(A) It returns the rank of a given >>rank(A)

rectangular matrix. >>ans=2
(no. of non-zeros rows)
trace(A) It returns the sum of diagonal elements >>trace(A)
of rectangular matrix. >>ans=5
inv(A) It returns the inverse of non-singular >>inv(A)
matrix (A-1)
norm(A) It returns the Euclidean norm of the >>norm(A)
given rectangular matrix ‘A’ of a given >>ans=4.4954
vector ‘A’
transpose(A) It returns the transpose of the specified >>A’
/A’ rectangular matrix. >>ans=
1 0
2 4

Generation Of Special Matrices:

The following MATLAB commands and functions which generate special matrices, are often
used in engineering computation involving matrices.

1. Zeros (m,n)
It returns a rectangular m× nmatrix of all zeros.

2. Ones (m,n)
It returns a rectangular m× n matrix of all ones.

3. eye (m,n)
It returns a rectangularm× nmatrix with ones on the main diagonal and zeros
elsewhere.

Matrix Manipulation:

1. Reshaping matrices as a vector:

All the elements of a matrix A can be grouped in to a single column vector b by the
command:
b= A(:)
Matrix A will be stored column wise in vector b i.e. first column is stored first, then
second column and so on.

2. Reshaping matrices as a differently sized matrix:

If a given matrix A is a p ×q matrix, it can be reshaped in to a new matrix m× nas long
as total elements of the two matrices are same i.e.
p ×q=m ×n
The command is as follows:
B=reshape ( A , m, n)

MATLAB MANUAL 20
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

When executed, this command will reshape the matrix A in to matrix B of size m× n.
The elements of the matrix B are taken column wise from matrix A.

3. Expanding the matrix size:

When the single element or a few elements are entered in a matrix, MATLAB creates the
matrix of proper dimensions to accommodate the entered element (elements).
Remaining unspecified elements are assumed to be zero.
Commands are:

 A ( i. j ) =aij (for single element)

 B ( i. a j :b j ) =¿] (for ith row and specified columns)
 C ( ai : bi , j )=¿] (for specified rows and jth column)

4. Appending a Row/Column to a Matrix:

Sometimes it is required that a column or a row to be added to a given matrix.

1. Appending Column:
A column can be appended using following command:
A=[ A x ] where x=[ a ; b ; c ]
2. Appending Row:
A row can be appended using following command:
A=[ A ; y ] where y=[ a b c ]
Semicolon is used to append row. In appending rows and columns, take care of
dimensions of matrix.

Deleting a Row/Column of a Matrix:

Rows and Columns of a matrix can be deleted by setting the corresponding row or column
vector equal to null vector i.e. pair of empty square brackets [ ], without any element within
brackets.

Commands are as follows:

 A ( i, : )=[ ]

It deletes all columns of ith rows

 A ( : , j )=[ ]

It deletes all rows of j th column.

 A ( i, j )=[ ]

This type of command will give error, because single element cannot be deleted from a matrix.

 B(:, aij :bij )= [ ]

MATLAB MANUAL 21
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

This will delete the a thto b thcolumn to all rows.

 C ( aij :b ij , : )=[ ]

This will delete the a th to b throws to all columns.

Random Numbers Generation:

There are many physical processes and engineering applications that require the use of
random numbers in the development of the solution.

MATLAB has two commands that can be used to assign random numbers to variables

 The rand command

 The randn command

The rand Command

Command Description Example
rand Generates a single random numbers >>rand
between 0 and 1. >>ans=0.9501
rand(1,n) Generates an ‘n’ elements row vector of >>a=rand(1,3)
random number between 0 and 1. >>a=
0.4565
0.0185
0.8214
rand(n) Generates an nxn matrix with random >>b=rand(3)
numbers between 0 and 1. >>b=
0.73 0.93 0.89
0.17 0.91 0.05
0.40 0.41 0.35
rand(m,n) Generates an mxn matrix with random >>c=rand(2,3)
numbers between 0 and 1. >>c=
0.20 0.60 0.19
0.19 0.27 0.01
randn(n) Generates an nxn matrix containing
normal random numbers with mean 0
and variance of 1.

randn(m,n) Generates an mxn matrix containing

normal random numbers with mean 0
and variance of 1.

Using range within matrix:

We can add equally spaced elements as much as we want in rows. This will display equally
spaced data in a matrix. Care should be taken while selecting ranges. If dimensions of each

MATLAB MANUAL 22
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

row do not match, it will not be executed and displays error.

Command is as follows:

A=[ k 1 :l1 :m 1 ; k 2 :l 2 : m2 ; k 3 :l 3 :m 3 ; … ]

Or

A=¿

Keep ‘n’ same for all rows. It will generate same number of elements in each row. The initial
and final points can be taken any value but initial point must be less than the final point.

Polynomials in Matrix Form:

 Poly(A):

It returns the coefficients of the characteristics polynomial ‘p’ of the matrix ‘A’ i.e. the
coefficient of the determinant of matrix [ xI−A ], where I is the identity matrix. The syntax is
given as:
p= poly (A )
When executed gives the co-efficient of the polynomial ‘p’ in descending power of ‘x’ as:

p=¿

This applies that the characteristic polynomial is x 2−5 x+ 4

 eig(A):

It returns the eigen values of matrix ‘A’ where ‘A’ is any given matrix.

 [ V , x ] =¿ eig(A):
It returns the eigen vectors ‘V’ and eigen values shown on the main diagonal of the

matrix ‘x’.

ILLUSTRATIVE PROGRAMS:

Task#1:
2

[]
Given the row vector P = [ 1−2 5 ] and a column Q = 4 . Solve the following:
8
a) P*Q
b) Q*P
c) PT*QT

MATLAB MANUAL 23
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

d) Sum and mean of all the components of vectors P & Q.

e) Number of elements in vectors P & Q
f) Maximum & Minimum values in vectors P & Q.
g) Product of elements of vectors P & Q.
h) Arrange the elements of vectors P & Q in ascending order.

Program:
%%
p=[1, -2, 5]; % row vector defined
q=[2; 4; 8]; % column vector defined
p*q % product of vectors
q*p % product of vectors
p'*q' % product of transpose of vectors

sum(p)
sum(q)

mean(p)
mean(q)

length (p) % counts number of elements in a vector

length (q)

max(p)
max(q)

min(p)
min(q)

prod(p) % product of elements in a vector

prod(q)

sort(p) % arranging elements in ascending order

sort(q)

Output:
ans =

1 3

ans =

3 1

ans =

MATLAB MANUAL 24
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

34

ans =

2 -4 10

4 -8 20

8 -16 40

ans =

2 4 8

-4 -8 -16

10 20 40

ans =

ans =

14

ans =

1.3333

ans =

4.6667

ans =

ans =

ans =

ans =

ans =

-2

MATLAB MANUAL 25
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

ans =

ans =

-10

ans =

64

ans =

-2 1 5

ans =

Task#2:
2 −4 6 −8
Given the matrix A= 1
[ 3 5
]
7 Write MATLAB statement to obtain:
2 12 30 56
a) All the elements of all rows but first column.
b) All the elements of first row but all columns.
c) Elements in the 2nd row & 3rd column.

Program:

A=[2 -4 6 -8; 1 3 5 7; 2 12 30 56];

a=A(:,1) % part (a)

b=A(1,:) % part (b)
c=A(2,3) % part (c)

Output:
a=

MATLAB MANUAL 26
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

b=

2 -4 6 -8

c=

Task#3:
Generate the following 3*3 matrices of common use with the help of MATLAB functions:
a) Unity or Identity matrix.
b) Null matrix.
c) Matrix with all elements equal to one where P=3 & Q=4.

Program:

eye (3,3) % identity matrix of 3x3

zeros(3,3) % null matrix of 3x3
ones(3,4) % matrix of ones of 3x4

Output:

ans =
1 0 0
0 1 0
0 0 1

ans =
0 0 0
0 0 0
0 0 0
ans =
1 1 1 1
1 1 1 1
1 1 1 1

Task#4:
For the following matrices, find:
a) Determinants.
b) The inverse of each, if they exist & the corresponding product AA -1, BB-1

MATLAB MANUAL 27
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

 A= [ 28 −14 ]
3 0 −2

[
B= −1 1 2
0 1 −2 ]
1 2 3

[
C= 0 −1 1
1 0 1 ]
0 2 0

[
D= 3 −3 4
0 2 6 ]
Program:

A=[2,4; 8,-1];
B=[3,0,-2; -1,1,2; 0,1,-2];
C=[1,2,3; 0,-1,1; 1,0,1];
D=[0,2,0; 3,-3,4; 0,2,6];
det(A)
det(B)
det(C)
det(D)
inv(A)
inv(B)
inv(C)
inv(D)
A*inv(A)
B*inv(B)
C*inv(C)
D*inv(D)

Output:
ans =
-34

ans =
-10

ans =
4

ans =
-36

ans =

MATLAB MANUAL 28
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

0.0294 0.1176
0.2353 -0.0588

ans =
0.4000 0.2000 -0.2000
0.2000 0.6000 0.4000
0.1000 0.3000 -0.3000

ans =

-0.2500 -0.5000 1.2500

0.2500 -0.5000 -0.2500
0.2500 0.5000 -0.2500

ans =
0.7222 0.3333 -0.2222
0.5000 0 0
-0.1667 0 0.1667

ans =
1 0
0 1

ans =
1 0 0
0 1 0
0 0 1

ans =
1 0 0
0 1 0
0 0 1

ans =
1.0000 0 0
-0.0000 1.0000 0
0 0 1.0000

Task#5:
Determine the ranks of the following matrices:
4 3

[ ]
a) L = 7 2
4 0

MATLAB MANUAL 29
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

b) M = [ 32 5 0 9
3 6 4 ]
1 2 3

[ ]
c) N= 4 5 6
7 8 9

Program:

L=[4 3; 7 2; 4 0];
M=[3 5 0 9; 2 3 6 4];
N=[1 2 3; 4 5 6; 7 8 9];
a=rank(L)
b=rank(M)
c=rank(N)

Output:
a=
2
b=
2
c=
2

Task#6:
Obtain the following products:
a) AB
b) BA
c) ATA
d) BTB
For the matrices given below:
2 4
A=[−1 1 3
4 5 6 ] ,
[ ]
B= 6 8
10 12
Program:
A=[-1 2 3; 4 5 6];
B=[2 4; 6 8; 10 12];
size(A)
size(B)
A*B
B*A
A'*A

MATLAB MANUAL 30
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

B'*B

Output:
ans =
2 3

ans =
3 2

ans =
40 48
98 128

ans =
14 24 30
26 52 66
38 80 102

ans =
17 18 21
18 29 36
21 36 45

ans =
140 176
176 224

Task#7:
Evaluate using array operations:
a) A+B
b) A-B
c) A^3 & B^4
d) A/B
5 −8 10 −2

[
A= 0 4
2 4
3 −4
6 8 ]
3 −3 −6 5

[
B= 8
0 4
0 5 4
3 2 ]
Program:

A=[5 -8 10 -2; 0 4 3 -4; 2 4 6 8];

MATLAB MANUAL 31
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

B=[3 -3 -6 5; 8 0 5 4; 0 4 3 2];
A+B
A-B
A.^3
B.^4
A./B

Output:
ans =
8 -11 4 3
8 4 8 0
2 8 9 10
ans =
2 -5 16 -7
-8 4 -2 -8
2 0 3 6

ans =
125 -512 1000 -8
0 64 27 -64
8 64 216 512

ans =
81 81 1296 625
4096 0 625 256
0 256 81 16

ans =
1.6667 2.6667 -1.6667 -0.4000
0 Inf 0.6000 -1.0000
Inf 1.0000 2.0000 4.0000

Task#7:
3 4 5 6
Given 3*4 matrix Q¿ 6
[ ]
7 8 9 , using MATLAB commands:
9 10 11 12
a) Delete the 1st row of matrix Q.
b) Delete the 1st & 2nd column of all the rows of matrix Q.
c) Extract the 2*2 sub-matrix from Q.
d) Replace the elements Q (2,2) with 800.

Program:

MATLAB MANUAL 32
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

%%
Q1=[3 4 5 6;6 7 8 9;9 10 11 12];
Q1(1,:)=[ ] % Delete the 1st row of matrix Q.

%%
Q2=[3 4 5 6;6 7 8 9;9 10 11 12];
Q2(:,1:2)=[ ] % Delete the 1st row 1st & 2nd column of all the rows of matrix

%%
Q3=[3 4 5 6;6 7 8 9;9 10 11 12];
Q3(2:3,2:3) % Extract the 2*2 sub-matrix from Q.

%%
Q4=[3 4 5 6;6 7 8 9;9 10 11 12];
Q4(2,2)=800 % Replace the elements Q (2,2) with 800.

Output:

Q1 =
6 7 8 9
9 10 11 12

Q2 =
5 6
8 9
11 12

Q3 =
7 8
10 11

Q4 =
3 4 5 6
6 800 8 9
9 10 11 12

Task#8:

3 4 5 1
Given the matrix P:

a) Reshape this matrix as a:

[
P= 5 6 7 2
7 8 9 3 ]
MATLAB MANUAL 33
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

 Column vector.
 (4*3) matrix.
 (6*2) matrix

[]
b) Given a column vector, x = 1 & a row vector, y = [ 0 2 3 4 ] , append these
2
column/row vector to the matrix P, given above.
Program:
%%
p=[3 4 5 1; 5 6 7 2; 7 8 9 3];
a=p(:)
b= reshape(p,4,3)
c= reshape(p,6,2)

%%
p=[3 4 5 1; 5 6 7 2; 7 8 9 3];
x=[0;1;2];
y=[0,2,3,4];
p1=[p x]

%%
p=[3 4 5 1; 5 6 7 2; 7 8 9 3];
y=[0,2,3,4];
p2=[p;y]

Output:
a=
3
5
7
4
6
8
5
7
9
1
2
3
b=
3 6 9
5 8 1

MATLAB MANUAL 34
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

7 5 2
4 7 3
c=
3 5
5 7
7 9
4 1
6 2
8 3
p1=
3 4 5 1 0
5 6 7 2 1
7 8 9 3 2
p2 =
3 4 5 1
5 6 7 2
7 8 9 3
0 2 3 4

Task#9:
Consider matrix A is given by:
2 5 6

[
A= 3 4 10
7 11 8 ]
Show the use of the following MATLAB command:
a) Diagonal (A), A is a matrix.
b) Diagonal (A, K), A is a matrix.
c) Polynomial (A).
d) Eigen (A).

Program:
A=[2 5 6; 3 4 10; 7 11 8];
diag(A)
diag(A,2)
poly (A)
eig(A)

Output:
ans =
MATLAB MANUAL 35
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

2
4
8
ans =
6
ans =
1.0000 -14.0000 -111.0000 -104.0000
ans =
19.8545
-1.1022
-4.7523

Task#10:
Express the following sets of algebraic equations in matrix form, Ay=b
x + y−z=2
−x +3 y−z=2
3 x+ 5 y −2 z =2
For this set of equations;
a) Find the inverse of the matrix, if it exists.
b) Obtain the solution of the variables x, y & z.
c) Find the Eigen values and Eigen vectors of matrix A.
d) Find the trace of matrix A.
e) Find the transpose of matrix A.

Program:
A=[1 1 1;-1 3 -1;4 -4 0];
B=[4;4;0];
X=inv(A)*B
[v x]=eig(A)
trace(A)
transpose(A) % second method for transpose is A’

Output:

X=
2
2
0

MATLAB MANUAL 36
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

v=
-0.0000 0.7071 -0.3487
0.7071 0.7071 0.1162
-0.7071 0.0000 0.9300
x=
4.0000 0 0
0 2.0000 0
0 0 -2.0000
ans =
4
ans =
1 -1 4
1 3 -4
1 -1 0

Task#11:
Create a 4*4 matrix of random numbers, multiply all the elements of matrix by 10 & then
round off all the elements to integers using appropriate commands.

Program:

a=rand(4)
b=10*a
round_off=round (b)

Output:
a=
0.8147 0.6324 0.9575 0.9572
0.9058 0.0975 0.9649 0.4854
0.1270 0.2785 0.1576 0.8003
0.9134 0.5469 0.9706 0.1419

b=
8.1472 6.3236 9.5751 9.5717
9.0579 0.9754 9.6489 4.8538
1.2699 2.7850 1.5761 8.0028
9.1338 5.4688 9.7059 1.4189

round_off =

MATLAB MANUAL 37
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

8 6 10 10
9 1 10 5
1 3 2 8
9 5 10 1

LAB#4:
POLYNOMIALS IN MATLAB

The learning objectives are:

 Writing a polynomial in MATLAB

 Perform different operations on polynomials
 Generate a polynomial from given roots

MATLAB MANUAL 38
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

 Obtain a characteristics polynomial for a given Matrix

 Compute the derivative and integration of a polynomial
 Fit a polynomial curve for a given set of point
 Evaluate polynomial with matrix argument

Introduction:
A polynomial is an expression in which a finite number of constants and variables are
combined using addition subtraction multiplication and non-negative whole number
exponents (raise to power). Thus the following expressions are all considered polynomials.

 x 2+ 2 x−7
 x 4 −7 x 3
 x
The MATLAB software provides built-in functions for standard operation on polynomials such
as arithmetic operation like multiplication and division etc., evaluation of roots,
differentiation, and integration and so on. In addition there are functions for more advance
applications such as curve fitting. Polynomials have a very significant role to play in
engineering system analysis. The polynomial functions are stored in a polyfun directory in
MATLAB.

Entering a Polynomial:
In MATLAB, a polynomial is represented by a row vector. Polynomials are entered with
variable term arranged in descending order of its power. Only the coefficient of variable terms
are entered i.e. variables are not entered. The left most term is the highest power term of the
variable and the right most term is the constant term with variables raise to the power zero.
The total number of term is always one more than the order of the polynomial. The coefficient
of polynomial is entered as an element of a row vector and can be separated from each other
by blank space or comma. For example:

p ( x ) =x 4 +3 x3 −15 x 2−2 x +9

This can entered in MATLAB as a row vector ‘P’ as column.

p=[1 3 −15−2 9] or p=[1,3 ,−15 ,−2,9 ]

MATLAB can interrupt a vector of length (n+1) as an nth order polynomial. Thus, if there are
missing terms, zeroes must be entered in the appropriate places in the vector.

e.g. y=x 3 +1

In MATLAB: y=[10 0 1] or y=[1,0,0,1]

Two zeroes have been included in a vector to account for missing powers of x, i.e. x 2 and x .

Arithmetic Operations on Polynomials:

Operations Syntax Description

MATLAB MANUAL 39
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

X and y are two vectors for two

Addition of two polynomials and z vector is their resultant.
Z=x+ y
polynomials The dimension (size) of both polynomials
vectors must be same for addition.

For subtraction also, the size/dimensions

of both vectors must be same. If the
Subtraction of Two
Z=x− y number of terms of two polynomials is not
Polynomials
equal, then zero is padded as coefficient to
make dimensions of vector same.

X and y are the vectors of coefficients of

polynomials to be multiplied and z
contains the coefficient of the resultant
Multiplication of Two conv ( x , y) polynomials.
Polynomials Or z=conv ( x , y) If there are more than two polynomials,
then the product is obtained by
performing multiplication of two of the
polynomials at one time.

Z = Dividend vector,
Division of Two y = Divisor vector,
[ x , r ]=deconv (z , y )
Polynomials x = Vector of quotients obtained,
r = vectors of remainders obtained.

Useful Commands for Polynomials:

List Of Commands Description

polyval( p , x) Evaluate polynomials at ‘x’
roots ( p) Computes the roots of polynomial
conv ( p , q) Multiplication of two polynomials
deconv ( p , q) Division of two polynomials
poly(r ) Converts given roots to polynomial
polyder( p) Differentiate a polynomial
Integrates polynomial analytically, ‘c’ is constant
polyint ( p , c ) of integration. If it is not mentioned in command,
1 is taken by default at place of ‘c’.
polyfit (x , y , k ) Fits a polynomial to given data

ILLUSTRATIVE PROGRAMS:
Task#1:
Obtain the Values of the following Polynomials at x=-2, 4, 5 and 2+3i?
 P(x) = x2 +4x +1.

MATLAB MANUAL 40
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

 P(x) = x2 + 1
 P(x) = x2 +4x + 2.

Program:

x=[-2 , 4.5 , 2+4*i];

P1 = [1 4 1];
P2 = [1 0 1];
P3 = [1 4 2];
V1=polyval(P1,x)
V2=polyval(P2,x)
V3=polyval(P3,x)

Output:
V1 =
-3.0000 39.2500 -3.0000 +32.0000i

V2 =
5.0000 21.2500 -11.0000 +16.0000i

V3 =
-2.0000 40.2500 -2.0000 +32.0000i

Task#2:
Find the Roots of the following Polynomials?
 P1 = x2 + 3x +2.
 P2 = x4 + 4x3 + 9x +2.
 P3 = x5 + 1.
 P4 = x3 + 3x2 + 4x + 1.

Program:
P1=[1 3 2];

P2=[1 4 0 9 2];
P3=[1 0 0 0 0 1];
P4=[1 3 4 1];
r1=roots(P1)
r2=roots(P2)
r3=roots(P3)
r4=roots(P4)

Output:
r1 = -2
-1

r2 = -4.4347

MATLAB MANUAL 41
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

0.3263 + 1.4012i
0.3263 - 1.4012i
-0.2179

r3 = -1.0000
-0.3090 + 0.9511i
-0.3090 - 0.9511i
0.8090 + 0.5878i
0.8090 - 0.5878i

r4 = -1.3412 + 1.1615i
-1.3412 - 1.1615i
-0.3177

Task#3:
Multiply the following Polynomials?
 P1 = x + 3.
 P2 = x + 6.
 P3 = x2 + 4x +6.
Program:
P1=[1 3];
P2=[1 6];
P3=[1 4 6];
P4=conv(P1,P2)
P5=conv(P3,P4)

Output:
P4 =
1 9 18
P5 =
1 13 60 126 108

Task#4:
Divide Polynomial P by q?
a) P = x3 + 4x +10. Q = x3 + 3x2 + 4x + 2.
b) P = 5x2 + 8x + 2. Q = x2 + x + 2.
c) P = 2x2 + 3x + 1. Q =x2 + x – 1.

Program:

P1=[1 0 4 10];
Q1=[1 3 4 2];
[z1,r1]=deconv(P1,Q1) % z=Quotient , r1=Remainder, P/Q (P divided by Q)

P2=[5 8 2];

MATLAB MANUAL 42
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Q2=[1 1 2];
[z2,r2]=deconv(P2,Q2)

P3=[2 3 1];
Q3=[1 1 -1];
[z3,r3]=deconv(P3,Q3)

Output:

z1 =
1

r1 =
0 -3 0 8

z2 =
5

r2 =
0 3 -8

z3 =
2

r3 =
0 1 3

Task#5:
Roots of Polynomial are given below find the Corresponding Polynomials?
a) r1 = 2, 4.
b) r2 = -2, -4.
c) r3 = -2, 4.
d) r4 = -2, 4, 6.

Program:
r1= [2; 4];
P1= poly(r1)

r2 = [-2,-4];
P2 = poly(r2)

r3 = [-2;4];
P3 = poly(r3)

r4 = [-2;4;6];
P4 = poly(r4)

MATLAB MANUAL 43
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Output:
P1 =
1 -6 8

P2 =
1 6 8

P3 =
1 -2 -8

P4 =
1 -8 4 48

Task#6:
Find Polynomial of Degree 5 for a given Data, for x=12?

X: 2 5 7 9 11 13
Y: 24 156 384 580 650 780

Program:
n=5; % n shows Degree of Polynomial
x=[2 5 7 9 11 13];
y=[24 156 384 500 650 780];
P=polyfit(x,y,n)
x=12;
P12=polyval(P,x)

Output:
P = 1.0e+003 * -0.0001 0.0048 -0.0739 0.5270 -1.6181 1.6700
P12 = 741.1275

Task#7:
Evaluate the Derivative and Integral of the following Polynomials?
 P1 = x4 + 4x3 + 10x2 + 20x + 15, c = 1.
 P2 = x3 + x + 10, c = 3.
 P3 = x + 5, c = -3.

Program:
P1=[1 4 10 20 15];
c1=1;
P2=[1 0 1 10];
c2=2;
P3=[1 5];
c3=-3;

MATLAB MANUAL 44
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

D1=polyder(P1) %Derivative
D2=polyder(P2)
D3=polyder(P3)

I1=polyint(P1,c1) %Integration
I2=polyint(P2,c2)
I3=polyint(P3,c3)

Output:

D1 =
4 12 20 20
D2 =
3 0 1
D3 =
1

I1 =
0.2000 1.0000 3.3333 10.0000 15.0000 1.0000
I2 =
0.2500 0 0.5000 10.0000 2.0000
I3 =
0.5000 5.0000 -3.0000
Task#8:
Express the following system of linear equations into Matrix form and obtain the solution of
the variables x, y, z.
x + y−z=2
−x +3 y−z=2
3 x+ 5 y −2 z =2
Program:

A= [1 1 1; -1 3 -1; 4 -4 0];
B= [4; 4; 0];
X= inv(A)*B

Output:

X=
2
2
0

MATLAB MANUAL 45
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

LAB#5:
CONSTRUCTION OF TWO DIMENSIONAL GRAPHS IN
MATLAB
The learning objectives are:

 To learn the basic features of plotting a graph

 To draw multiple plots using ‘plot’, ‘hold’ and ‘line’ functions
 To learn about different style options to vary color, marker and line-style
 To provide a legend on the plot
 To draw multiple sub-graph windows using subplot functions
 To draw specialized plots such as logarithmic, polar, bar, histogram, pie, stairs and stem
plots

Introduction:
MATLAB provides very good tools for 2-dimensional plots, functions, tools for interactively
creating plots and the ability to export results to many popular graphics formats. It is possible
to customized plots by adding multiple axis, changing line colors and markers, adding
annotation and legends.

Vector data can be visualized with 2-D plotting functions such as line, area, bar and pie
charts, direction and the velocity plots, histograms, polygons and surfaces, scatter or bubble

MATLAB MANUAL 46
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

plots and animations. MATLAB provides functions for visualizing 2-D matrices. In MATLAB,
one can specify plot characteristics, such as viewing angle, perspective, lighting effect, light
source, location and transparency.

List of Commands:

1. General Graphics Command:

Function Name Syntax Description

plot plot(x , y) To draw a 2-D plot
or
plot(x,y,’style’,’linewidth’,a Style could be line color, line style or marker
) style. ‘a’ is constant for line width in cm let
say 2.
xlabel xlabel(‘string’) To label the x-axis
ylabel ylabel(‘string’) To label the y-axis
title title(‘string’) To give a title to the figure
grid grid on To show the grid line on the plot
text text(‘string’) To write illustrative text on the plot
gtext gtext(‘string’) To place text on the plot using mouse
axis axis(‘equal’) It sets equal scale on axis
axis(‘square’) It sets rectangular frame
axis(‘normal’) It resets to default values
axis(‘off’) It turn off all axis lines, tick-marks and
labels
axis tight It fits the graph in figure window.
legend legend(‘x’, ‘y’, ‘z’) To produce a boxed legend on the plot
subplot subplot(r , c , l) To subdivide graphic window into many sub-
or graphs.
subplot rcl Use this command before ‘plot’ command
r=number of rows
c=number of columns
l=location of subplot on
figure window
axis axis( [ x a x b y a y b ] ) To fix range of points on axis
pause pause(n) To fix the time in seconds for drawing
n=number of seconds multiple graphs
hold hold on Use plot command and then hold command
to plot multiple graphs on single axis
line line(x , y) Use it after plot command to draw multiple
graphs on single axis.

2. 2-D Plot Functions:

Function name Description

semi log x ( x , y ) Plot with logarithmic scale x-axis and linear scale y-axis
semi log y (x , y ) Plot with linear scale x-axis and logarithmic scale y-axis
MATLAB MANUAL 47
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

log log(x , y) Plot with logarithmic scale x-axis and y-axis

polar(θ , r ) To draw a plot on polar axis coordinates
area( x , y) Same as ‘plot’ command except area under curve is filled with
color
stairs(x , y ,' style ') To give stair case plot
contour ( x , y , ' style ' ) To give contour of the given matrix

Style Options:

1. Color Style Options:

y = yellow m = magenta c = cyan

r = red g = green b = blue
w = white k = black

2. Marker Style Options:

. = point o = circle x = x-mark

+ = plus * = star s = square

d = diamond v = triangle (down) h = hexagram

^ = triangle (up) < = triangle (left)

> = triangle (right) p = pentagram

3. Line Style Options:

For Solid line= - For Dashed line= --

For Dotted line= . For Dashed dotted line= -.

ILLUSTRATIVE PROGRAMS:

Task#01:
Plot a curve for the equation
y=x 3 +2 x 2−5 ,−10< x <10.
Use labels, title and gtext command to write this equation on the curve plotted?

Program:
x=-[Link];
y=x.^3+2*x.^2-5;
plot (x,y) % command for drawing graph
xlabel ('X-Axis'); % for labeling x-axis
ylabel ('Y-Axis'); % for labeling y-axis
title ('2D Graph'); % for giving title on top in figure window

MATLAB MANUAL 48
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

gtext ('y=x^3+2*x^2-5'); % for dropping text inside figure window on graphical area
grid on % for drawing grid lines inside graphical area

OUTPUT:

Task#02:
Write a program to illustrate the use of area command with specified limits, for plot of
function y=cos x for x=0 ¿ π with step size of 0.02?

Program:
x=0:0.02:2*pi;
y=cos(x);
area (x,y) % command for drawing graph having filled area.
xlabel ('X-Axis');
ylabel ('Y-Axis');
title ('Graphing s single plot with filled area using Area Command');
gtext ('Filled Area');

OUTPUT:

MATLAB MANUAL 49
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Task#03:
Write a program to draw the curves for functions in a single figure window using single plot
command. Illustrate the use of different styles of curves which can be drawn using plot
command?
 y 1=sin x
 y 2=sin (x−0.35)
 y 3=sin (x−0.55¿) ¿
 y 4 =sin( x−0.85)

Program:
x=0:pi/10:2*pi;
y1=sin(x);
y2=sin(x-0.35);
y3=sin(x-0.55);
y4=sin(x-0.85);

plot (x,y1,'r-o' , x,y2,'b-o' , x,y3,'c-o' , x,y4,'k-o','linewidth',2)

xlabel ('X-Axis');
ylabel ('Y-Axis');
title ('Graphing multiple plotsin a single figure window using plot command');
legend (‘y1=sin(x)’,’y2=sin(x-0.35)’,’y3=sin(x-0.55’),’y4=sin(x-0.85)’)
grid on;
axis tight

OUTPUT:

MATLAB MANUAL 50
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Task#04:
Write a program to draw the curves for function
 y=200sin x∗e−0.5 x
 y=0.8 sin10 x∗e−0.5 x
Illustrate how the “Legend” command is used to indicate the different styles used for the
different curves plotted on same graph?

Program:
x=0:0.01:100;
y1=200*exp(-0.5*x).*sin(x);
y2=0.8*exp(-0.5*x).*sin(10*x);
plot (x,y1,'b',x,y2,'r')
xlabel ('X-Axis');
ylabel ('Y-Axis');
title ('Graphing multiple plots using plot command');
legend ('y1=200*exp(-0.5*x).*sin(x)','y2=0.8*exp(-0.5*x).*sin(10*x)');
grid on
axis tight

OUTPUT:

MATLAB MANUAL 51
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Task#05:
Write a program to draw the curves for functions in a single figure window using using Hold
and Plot Command?
 y 1=sin x
 y 2=sin (x−0.35)
 y 3=sin (x−0.55¿) ¿
 y 4 =sin( x−0.85)

Program:
x=0:pi/10:2*pi;
y1=sin(x);
plot (x,y1, 'r-.o','linewidth',2)
pause (2)

hold on

y2=sin(x-0.35);
plot (x,y2,'b-.o','linewidth',2)
pause (2)

y3=sin(x-0.55);
plot (x,y3,'c-.o','linewidt',2)
pause (2)

y4=sin(x-0.85);
plot (x,y4,'k-.o','linewidth',2)
pause (2)

hold off

xlabel ('X-Axis');

MATLAB MANUAL 52
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

ylabel ('Y-Axis');
title (‘Graphing multiple curves using Hold Command’);
grid on
axis tight

OUTPUT:

Task#06:
Write a program to draw the curves for functions in a single figure window using Line
Command?
 y 1=sin x
 y 2=sin (x−0.35)
 y 3=sin (x−0.55¿) ¿
 y 4 =sin ( x−0.85 )

Program:
x=0:pi/10:2*pi;
y1=sin(x);
y2=sin(x-0.35);
y3=sin(x-0.55);
y4=sin(x-0.85);

plot (x,y1, 'r-o','linewidth',2)

pause (2)

plot(x,y2, 'b-box','linewidth',2)
pause (2)

plot(x,y3, 'g-*','linewidth',2)
pause (2)

MATLAB MANUAL 53
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

plot(x,y4, 'k-d','linewidth',2)
pause (2)

xlabel ('X-Axis');
ylabel ('Y-Axis');
title ('Graphing multiple plots using Line Command');
legend (‘y1=sin(x)’,’y2=sin(x-0.35)’,’y3=sin(x-0.55’),’y4=sin(x-0.85)');
grid on
axis tight

OUTPUT:

Task#07:
Divide the Figure window into four sub-windows and plot the following functions?
 Plot V versus I where V=4I and I=1,2,3,4 on the upper left sub-window.
 Plot y versus x where y=x 2 and x=1,2,3,4 on the upper right sub-window.
π
 For t=0 :2 π in steps of t= . Plot sin t versus t on the lower left window
60
π
 For t=0: :2 π . Plot cos t versus t on the lower right window.
30

Program:
I=[Link];
V=4*I;
subplot (2,2,1)
plot(I,V)
title('V Versus I');

x=[Link];
y=x.^2;

MATLAB MANUAL 54
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

subplot (2,2,2)
plot(x,y)
title('y Versus x');

t=0:pi/60:2*pi;
Y=sin(t);
subplot (2,2,3)
plot(t,Y)
title('t Versus Sint');

t=0:pi/60:2*pi;
Y1=cos(t);
subplot (2,2,4)
plot(t,Y1)
title('t Versus Cost');

OUTPUT:

Task#08:
Plot the following functions on the polar plot:
−π π
f (θ)=sin 4 θ , < θ< , Where θ is in radians?
2 2

Program:
theta=-pi/2:pi/40:pi/2;
f=sin(4*theta);
polar (theta,f,'r-')
title ('Graphing using Polar Function')

OUTPUT:

MATLAB MANUAL 55
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Task#09:
a) Draw the stairs plot to show the function y=x 2 where −2 ≤ x ≤2
b) Draw the stairs plot for the following data:
X=[0 2 3 4 5 7]
Y=[5 -1 8 4 7 3]

Program:
%%
X=-2:0.1:2;
Y=X.*X;
stairs(X,Y,'c','linewidth',2) % c is for cyan color
xlabel ('X-Axis');
ylabel ('Y-Axis');
title ('Graphing stairs plot using stairs Command');
grid on
%%
X=[0 2 3 4 5 7];
Y=[5 -1 8 4 7 3];
stairs(X,Y,'*-','linewidth',2)
xlabel ('X-Axis');
ylabel ('Y-Axis');
title ('Graphing stairs plot using stairs Command');
grid on

OUTPUT:

MATLAB MANUAL 56
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Graphing
Graphingstairs
stairs plot
plot using
usingstairs
stairs Command
Command
44

3.5
3.5

33

2.5
2.5

-Axxisis
22

YY-A
1.5
1.5

11

0.5
0.5

00
-2
-2 -1.5
-1.5 -1
-1 -0.5
-0.5 00 0.5
0.5 11 1.5
1.5 22
X-Axis
X-Axis

Graphing
Graphingstairs
stairs plot
plot using
usingstairs
stairs Command
Command
88

77

66

55

44
-Axxisis
YY-A

33

22

11

00

-1
-1
00 11 22 33 44 55 66 77
X-Axis
X-Axis

Task#10:
Draw the stem plot for the following data:
X=[0 1 2 3 4 5 6]
Y=[3 -1 6 -4 5 2 3]

Program:
%%
X=[0 1 2 3 4 5 6];
Y=[3 -1 6 -4 5 2 3];
stem(X,Y,'k','linewidth',2.5)
xlabel ('X-Axis');
ylabel ('Y-Axis');
title ('Graphing stairs plot using stairs Command');
grid on

MATLAB MANUAL 57
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

OUTPUT:
Graphing
Graphing stairs
stairs plot
plot using
using stairs
stairs Command
Command
66

55

44

33

22
Y-Axis
Y-Axis

11

00

-1
-1

-2
-2

-3
-3

-4
-4
00 11 22 33 44 55 66
X-Axis
X-Axis

LAB#6:
CONSTRUCTION OF THREE DIMENSIONAL GRAPHS
IN MATLAB
The learning objectives are:

 To learn the basic features of plotting a 3D graph

 To draw multiple plots using ‘plot’, ‘hold’ and ‘line’ functions
 To learn about different style options to vary color, marker and line-style in 3D graphs.

Introduction:
MATLAB provides very good tools for 3-dimensional plots, 3-D volume visualization, functions,
tools for interactively creating plots and the ability to export results to many popular graphics
formats. It is possible to customized plots by adding multiple axes, changing line colors and
markers, adding annotation and legends.
MATLAB MANUAL 58
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

MATLAB provides functions for visualizing 3-D scalar and the 3-D vector data. In MATLAB,
one can specify plot characteristics, such as viewing angle, perspective, lighting effect, light
source, location and transparency. 3-dimenstional plotting function includes surface,
contour, mesh, image plots simple and easily understandable. This makes graphics one of the
most powerful features of MATLAB.

List of 3D Graphic Commands:

Function Description
plot 3 (x , y , z ) To plot a simple 3D plot.
[ X , Y ] =meshgrid (x , y) To give array for mesh grid
mesh(x , y , z) To plot wireframe mesh plot
surf (x , y , z) To plot surface plot
To show shadow or circles beneath
meshc ( x , y , z)
graph
meshz (x , y , z) To show curtons beneath graph
surfc( x , y , z) To plot shadow beneath surface
contour 3 ( x , y , z ) To give contour over the plane
To show lines of 3D graph like
waterfall ( x , y , z ) waterfall
stem 3(x , y , z) To draw the discrete surface as stems
comet 3( x , y , z ) To see graph while drawing

ILLUSTRATIVE PROGRAMS:

Task#01:
Given is y=sin x∧z =x+5 , for 0< x <5000
Obtain a 3D plot showing the matrix is (y,z) space as a factors of ‘x’?

Program:
x=0:2*pi:5000;
y=sin(x);
z=x+5;
plot3(x,y,z)
xlabel ('X-Axis');
ylabel ('Y-Axis');
zlabel ('Z-Axis');
title ('Plot 3D Graph');
grid on;

OUTPUT:

MATLAB MANUAL 59
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Task#02:
Given is y=sin x∧z =x+5 , for 0< x <5000
Obtain a 3D plot using a Comet Command?

Program:
x=0:2*pi:5000;
y=sin(x);
z=x+5;
comet3(x,y,z) %OR comet3(x,y,z) we see graph is drawn
xlabel ('X-Axis');
ylabel ('Y-Axis');
zlabel ('Z-Axis');
title ('Use of Comet Command');
grid on
axis tight

OUTPUT:

Task#03:
sin r
Create a surface plot of the function: z=
r
where r =√ x 2 + y 2 for−8< x <8∧−8< y < 8?

MATLAB MANUAL 60
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

PROGRAM:
x=-8:0.5:8;
y=-8:0.5:8;
[X,Y]=meshgrid(x,y); % x and y changed into X,Y(now we use these X and Y)
r=sqrt(X.^2+Y.^2);
Z=sin(r)./r;
surf(X,Y,Z)

OUTPUT:

Task#04:
sin r
Create a mesh plot of the function: z=
r
2 2
where r =√ x + y for−8< x <8∧−8< y < 8?

PROGRAM:
x=-8:0.5:8;
y=-8:0.5:8;
[X,Y]=meshgrid(x,y); %x and y changed into X,Y(now we use these X and Y)
r=sqrt(X.^2+Y.^2);
Z=sin(r)./r;
mesh(X,Y,Z)

OUTPUT:

MATLAB MANUAL 61
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Task#05:
Plot the function using “meshc”, “meshz” and “waterfall” commands using the function
sin r
z= where r =√ x 2 + y 2 for−8< x <8∧−8< y < 8. Create subplots for each command.
r

PROGRAM:
x=-8:0.5:8;
y=-8:0.5:8;
[X,Y]=meshgrid(x,y);
r=sqrt(X.^2+Y.^2);
Z=sin(r)./r;

subplot (3,1,1);
meshc(X,Y,Z)
title ('meshc Graph');

subplot (3,1,2);
meshz(X,Y,Z)
title ('meshz Graph');

subplot (3,1,3);
waterfall(X,Y,Z)
title ('waterfall Graph');

OUTPUT:

MATLAB MANUAL 62
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Task#06:
Write a program to illustrate the use of stem3 function to plot the following function:
y=x cos x∧z=cosx e x/5 +1 for 0 ≤ x ≤6 π

Program:
%%
x=0:pi/50:6*pi;
y=x.*cos(x);
z=exp(x/5).*cos(x)+1;

MATLAB MANUAL 63
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

stem3(x,y,z)
xlabel ('X-Axis');
ylabel ('Y-Axis');
zlabel ('Z-Axis');
title ('Use of stem3 Command');
grid on
axis tight

OUTPUT:

Use
Useofofstem3
stem3Command
Command

40
40
30
30
20
20
Z-Axis
Z-Axis

10
10
00

-10
-10
-20
-20

10
10
15
15
00 10
10
-10
-10 55
Y-Axis 00
Y-Axis X-Axis
X-Axis

LAB#7:
CONSTRUCTION OF STATISTICAL GRAPHS IN MATLAB

The learning objectives are:

 To learn the basic features of plotting a 3D graph

 To draw multiple plots using ‘plot’, ‘hold’ and ‘line’ functions
 To learn about different style options to vary color, marker and line-style
 To draw specialized plots such as bar, histogram, pie plots

MATLAB MANUAL 64
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Introduction:
MATLAB provides very good tools for 2-dimensional and 3-dimensional plots, 3-D volume
visualization, functions, tools for interactively creating plots and the ability to export results to
many popular graphics formats. It is possible to customized plots by adding multiple axis,
changing line colors and markers, adding annotation and legends.

In MATLAB, one can specify plot characteristics, such as viewing angle, perspective, lighting
effect, light source, location and transparency. This makes graphics one of the most powerful
features of MATLAB.

List of Commands:

Function name Description

' ¯ optio n' )
( x , y , style To draw vertical bar-graph
Or Style options means colors of bins
Grouped means multiple bar graph and stacked means
( x , y ,' grouped∨stacked
¯ ')
simple bar graph
barh(x , y) To draw horizontal bar-graph
hist ( x , y) To draw histogram
pie(x , y ) To draw pie plot
stem( x , y ) To draw a stem plot
rose( x , y) To draw a histogram in polar form
3̄(x , y) To plot a 3D vertical bar graph
3̄ h(x , y) To plot a 3D horizontal bar graph
pie3 (x , y ) To plot a 3D pie plot

ILLUSTRATIVE PROGRAMS:

Task#01:
Plot the bar graph for the data given as:
x=[01 2 3 4 5 6]∧ y=[10 15 25 20 3027 19] ?

PROGRAM:
x= [0 1 2 3 4 5 6];
y = [10 15 25 20 30 27 19];
bar (x,y,'c') % c is for cyan color
xlabel ('x-axis');
ylabel ('y-axis');
title ('Graph to show bar function');

OUTPUT:

MATLAB MANUAL 65
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Task#02:
Plot a bar graph to show the comparison of average temperature in city A, B and C for the for
the months from September to February, data is given as:

MONTHS CITY A CITY B CITY C

September 31 28 24
October 29 26 22
November 28 25 20
December 27 24 16
January 26 22 17
February 29 25 20

PROGRAM:
%%
x= [31,28,24; 29,26,22; 28.25,20; 27,24,16; 26,22,17; 29,25,20];
bar(x, 'grouped') % grouped/stacked
xlabel('Months from September to February');
ylabel('Temperature in Celsius');
legend ('City A’, ‘City B’, ‘City C');
title ('Comparison of Average Temperature in City A, B and C');

OUTPUT:

MATLAB MANUAL 66
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Comparison
ComparisonofofAverage
AverageTemperature
TemperatureininCity
CityA,
A,BBand
andCC
35
35
City
CityAA
City
CityBB
30
30 City
CityCC

25
25

Celsius
TemperatureininCelsius
20
20

Temperature
15
15

10
10

55

00
11 22 33 44 55 66
Months
Monthsfrom
fromSeptember
SeptembertotoFebruary
February

Task#03:
Plot the Histogram graph with 15 bins using the random data?

PROGRAM:
x=randn (200,1);
hist(x,15)
title ('Histogram to show the 200 random values in 15 bins')

OUTPUT:

Task#04:
Illustrate with the use of ‘hist’ function with number of bins specified by a vector P?

PROGRAM::
x = randn(300,1);
P = -2:0.1:2;
hist (x,P)
title ('Histogram to show the values of 300 and no of bins specified by Vector');

MATLAB MANUAL 67
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

OUTPUT:

Task#05:
Illustrate the use of ROSE Function?

PROGRAM:
X=randn(500,1)*pi; % generate 500 random data points
P=-pi:0.5:pi; % specifies the number of bins in a vector V
rose (x,P)
title ('Polar Histogram using Rose Function');

OUTPUT:

Task#06:
Illustrate the use PIE Function to show the concentration of different industries in the region
as per following data:
INDUSTRIES INDUSTRIAL UNITS
Cement 4
Textile 8
Software 20
Chemical 2
Telecom 7
Banking 10

Show the sliced and labeled pie chart?

MATLAB MANUAL 68
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

PROGRAM:
x=[4 8 20 2 7 10];
k=[1,0,0,1,1,0];
label=({'Cement','Textile','Software','Chemicel','Telecom','Banking'});
pie(x,k,label) %For 2D PIE Chart

OUTPUT:

Task#07:
Draw 3D bar graph for the data given as:
x= [0 1 2 3 4 5 6] and y = [10 15 25 20 30 27 19]?

PROGRAM:
x= [0 1 2 3 4 5 6];
y = [10 15 25 20 30 27 19];
bar3 (x,y,'g')
xlabel ('x-axis');
ylabel ('y-axis');
title ('Graph to show bar 3D function');

OUTPUT:

MATLAB MANUAL 69
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Task#08:
Illustrate the use of pie3 function to show different professions people groups as per the
following data?

Professions No. of Persons

Manager 12
Engineer 15
Professor 13
Doctor 11
Architect 6
Designer 10

PROGRAM:

x = [12 15 13 11 6 10];
k = [0 0 1 0 0 1];
label = ({'Manager', 'Engineer','Professor','Doctor','Architect','Designer'});
pie3 (x, k, label)
title ('Persons in a Group using pie3 Function');

OUTPUT:

MATLAB MANUAL 70
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

LAB#8:
CONTROL STRUCTURES IN MATLAB

The learning objectives are:

MATLAB MANUAL 71
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

 To learn conditional and unconditional flow control structures in MATLAB.

 Looping control structures.
 Branding control structures.
 Special MATLAB statements which control program execution sequence.

Introduction:
The instruction given in a program are executed in a sequence in which they are written,
however many times program applications demands that the sequence of execution of the
instruction be changed or a group of instruction can be repeated depending upon some
conditions, control structure are provided in a computer language to make care of such needs
of the program, use for testing a condition for true or false, repeated execution of statements
out of several different available possibilities.
Like other programming language, the MAT-Lab also has different control structures for
controlling the sequence of execution of statement of a program.
 First there is a conditional execution in which a condition is checked and then action is
taken as per the result of the condition. In another type of control structure
instructions are executed repeatedly until some conditions is met. This is called
Conditional looping, no of repetitions are not known.
 A group of instructions are executed repeatedly foe specified no of times are called
unconditional looping and in this case no of repetitions are already known.

Control structures are also available for selective execution of group of instructions, to solve
such type of executions; MAT-Lab has a number of control structures. There are 2 broad
categories of these control structures.

LOOPING STRUCTURES:
They cause specific group of instructions to be repeated for a fixed no of times or until the
given condition is specified, some of looping structures are for, while among others.

BRANCHING STRUCTURES:
They select the specific blocks of the program for execution depending upon the condition
satisfied, some of them are if, if-else among others.
Beside, these there are some other commands such as break, error; continue to control the
execution of such structures.

LOOPS:
They are classified depending upon the no in which the iterations in a program are controlled,
for and while loops. In for loop the no of repetitions are already known but in while loop
instructions are executed until condition satisfied.

FOR LOOP:

MATLAB MANUAL 72
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

It is used for unconditional looping, in it a group of statements say ‘n’ in no are executed for
a fixed no of times, which are known before the loop starts and is decided by an index, its
syntax are as follows:

for index = expression

statement # 1;

statement # 2;
.........................
statement # n;
end
next statement;

The index is assigned an initial value and then a group of statements are executed till ‘end’
statement. After completion of one loop the index is incremented by specifying and again
group of statements are executed until the provided value of index not exceeds the last value.
Then loop is terminated.

NESTED FOR LOOP:

When a for loop is embedded in another for loop, the 2 loops are called nested for loops.
The necessary condition is that one loop is completely nested in 2nd loop.

WHILE LOOP:
It is used for executing a group of statements repeatedly on the basis of condition. In
this loop the group of statements associated with the while statement are executed infinitely
until the condition true, its general form are:

while (condition)
statement # 1;

statement # 2;
........................
statement # n;
end

The group of statements may include same statement which may alter the result of the
condition, when condition is false the loop is terminated.

BRANCHES CONTROL STRUCUTRE:

They allow selection and execution of some specific statements depending upon the condition
They include if, if-else and switch etc. If is the conditional statement which allow some group
of statements to be executed if the given condition is true.

IF CONTROL STRUCTURE:
In if control structure a group of statements are executed only if the condition is true, some
are given blow:

MATLAB MANUAL 73
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

 If structure.
 If-else structure.
 If-else, if-else structure.
 Nested if structure.

IF STRUCUTRE:
This is simplest form of structure, it has only 1 group of statements which is executed if the
condition is true otherwise not, its syntax is:
If (expression)
statement # 1;

statement # 2;
........................
statement # n;
end

next statement;

where expression is a logic expression and give the result either true or false, if expression
give true result then the statements are executed till end otherwise not.

IF-ELSE STRUCUTRE:

In this form control structure has 2 statements if or else one for true and other for false, its
syntax is as follows:

If expression
statement # 1;

statement # 2;
........................
statement # n;
else
statement # 1;

statement # 2;
........................
statement # n;
end

if the expression is true then if statements are executed and else statements are skipped or if
not true then else statements are executed and if statements are executed.

IF-ELSE-IF STRUCTURE:

It is most use of if structure and it allow the user to check large no of statements and
expressions. Its syntax is as follows:

If (expression)

MATLAB MANUAL 74
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

statement # 1;

statement # 2;
........................
statement # n;
else if (expression)

statement # 1;

statement # 2;
........................
statement # n;
else
statement # 1;

statement # 2;
........................
statement # n;
end

next statement;

 If the expression give true condition then all of its statements are run and and other all
are skipped or if it give false then after end statements are run and if statements are
skipped.
 There are some multiple blocks in if else if else structure, the if else condition is tested
sequentially if the previous condition are false, once the true statement is found then at
block is executed and other all statements are skipped.

NESTED IF STRUCTURE:
Like other control structures, if structure can also be nested if structure are called nested,
when1 of the “if structure” lies entirely within in the domain of other if structure, its syntax is:
If (condition)
statement # 1;

statement # 2;
If (condition)
statement # 1;

statement # 2;
end
statement # 1;
statement # 2;
end

LAB#9:
OPEN ENDED LAB
MATLAB MANUAL 75
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

The learning objectives are:

LAB#10:

MATLAB MANUAL 76
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

CODE OF BISECTION METHOD AND REGULA-FALSI

METHOD FOR SOLUTION OF TRANSCENDENTAL
EQUATIONS IN MATLAB

The learning objectives are:

 How does Root Bracketing Methods work?
 Background of Bracketing Methods.
 Algorithm.
 Convergence sheet.

1. BI-SECTION METHOD:
ALGORITHM:
Newton’s method is a popular technique for the solution of nonlinear equations, but
alternative methods exist which may be preferable in certain situations. The Bisection method
is yet another technique for finding a solution to the nonlinear equation f(x) = 0, which can be
used provided that the function f is continuous. The motivation for this technique is drawn
from Bolzano’s theorem for continuous functions.

Theorem (Bolzano):
If the function f(x) is continuous in [a, b] and f(a).f(b) < 0 (i.e. the function f has values with
different signs at a and b), then a value c belongs to (a, b) exists such that f(c) = 0.

The bisection algorithm attempts to locate the value c where the plot of f crosses over zero, by
checking whether it belongs to either of the two sub-intervals [a, xm], [xm, b], where xm is the
midpoint xm = (a + b)/2
The steps to apply the bisection method to find the root of the equation f (x )=0 are

1. Choose
xℓ and
xu as two guesses for the root such that f (x ℓ )f ( x u )<0 , or in
other words, f (x ) changes sign between
x ℓ and x u .

MATLAB MANUAL 77
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

2. Estimate the root,

x m , of the equation f (x )=0 as the mid-point between
x ℓ and
x u as
x ℓ +x u
xm =
2

3. Now check the following

a) If f (x ℓ )f ( x m )<0 , then the root lies between xℓ and
x m ; then x ℓ=x ℓ and
x u=x m .
b) If f (x ℓ )f ( x m )>0 , then the root lies between x m and x u ; then x ℓ=x m and
x u=x u
c) If f (x ℓ )f ( x m )=0 ; then the root is
x m . Stop the algorithm if this is true.
4. Find the new estimate of the root
xℓ + xu
xm =
2

Find the absolute relative approximate error as

new old
xm - xm
|∈a| = | |× 100
x mnew
where

x new
m = estimated root from present iteration

x old
m = estimated root from previous iteration

5. Compare the absolute relative approximate error

|∈a| with the pre-specified relative
∈
error tolerance s . If
|∈a|>∈s , then go to Step 3, else stop the algorithm. Note one
should also check whether the number of iterations is more than the maximum
number of iterations allowed. If so, one needs to terminate the algorithm and notify the
user about it.

MATLAB CODE:

TASK:

A mass balance for a pollutant in a well-mixed lake can be written as:

dc
v =w−Qc−kv √ c
dt

Given the parameter values as follows:

MATLAB MANUAL 78
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

v=1× 106 m3
Q=1 ×10 5 m 3 / year
w=1 ×106 gm/ year
k =0.25 m 0.5 / g0.5 / year
The equation takes the form:

25 √ c+ 10 c−100=0

Solve for the steady-state concentration, find the root using Bisection method in MATLAB.

PROGRAM:
clc
close all
clear all

f=@(x) 25*sqrt(x)+10*x-100 % Let steady-state concentration = x

f(0)
f(5)
f(10)
fplot(f,[0,10])
xlabel('Steady-state concentration, c');
ylabel('f(c)');
title('BI-SECTION METHOD');
gtext('xr')
grid on

xa=0; % Initial guess

xb=5; % Initial guess

fprintf ('i xa xb xr f(xr) RE\n');

for i=1:20
xr=(xa+xb)/2; % Bisection formula

f_xr=25*sqrt(xr)+10*xr-100;
RE=abs((xr-xb)/xr)*100;

fprintf ('%6d %15.5f %15.5f %15.5f %15.5f %18.5f \n' ,

i,xa,xb,xr,f_xr,RE);

if f(xr)*f(xa)<0 % Root Bracketing

xb=xr;
xa=xa;

else
xa=xr;
xb=xb;

end %End of if-else Statement

end %End of For-Loop Statement
fprintf('\n The root of given equation = %f \n' , xr); % to print statement

MATLAB MANUAL 79
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

OUTPUT:
f=
@(x)25*sqrt(x)+10*x-100

ans =
-100
ans =
5.9017

ans =
79.0569

i xa xb xr f(xr) RE
1 0.00000 5.00000 2.50000 -35.47153 100.0000
2 2.50000 5.00000 3.75000 -14.08771 33.33333
3 3.75000 5.00000 4.37500 -3.95875 14.28571
4 4.37500 5.00000 4.68750 1.00159 6.66667
5 4.37500 4.68750 4.53125 -1.47067 3.44828
6 4.53125 4.68750 4.60938 -0.23261 1.69492
7 4.60938 4.68750 4.64844 0.38496 0.84034
8 4.60938 4.64844 4.62891 0.07630 0.42194
9 4.60938 4.62891 4.61914 -0.07813 0.21142
10 4.61914 4.62891 4.62402 -0.00091 0.10560
11 4.62402 4.62891 4.62646 0.03769 0.05277
12 4.62402 4.62646 4.62524 0.01839 0.02639
13 4.62402 4.62524 4.62463 0.00874 0.01320
14 4.62402 4.62463 4.62433 0.00392 0.00660
15 4.62402 4.62433 4.62418 0.00150 0.00330
16 4.62402 4.62418 4.62410 0.00030 0.00165
17 4.62402 4.62410 4.62406 -0.00031 0.00082
18 4.62406 4.62410 4.62408 -0.00000 0.00041
19 4.62408 4.62410 4.62409 0.00015 0.00021
20 4.62408 4.62409 4.62409 0.00007 0.00010

The root of given equation = 4.624085

GRAPH:

MATLAB MANUAL 80
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

BI-SECTION METHOD

50
xr

0
f(c)

-50

-100
0 2 4 6 8 10
Steady-state concentration, c

2. REGULA-FALSI METHOD:

ALGORITHM:
The steps to apply the false-position method to find the root of the equation f ( x )=0 are as follows.

x
Choose x L and U as two guesses for the root such that
f ( x L ) f ( x U ) <0
, or in other words,
f ( x ) changes sign between x L and x U .

Estimate the root, x r of the equation f ( x )=0 as

x U f ( x L )−x L f ( x U )
xr=
f ( x L ) −f ( xU )

Now check the following

If
f ( x L ) f ( x r ) <0
, then the root lies between x L and x r ; then x L =x L and x U =x r .

If
f ( x L ) f ( x r ) >0
, then the root lies between x r and
x U ; then x L =xr and x U =x U .

, then the root is x r . Stop the algorithm.

f ( x L ) f ( x r )=0
If

Find the new estimate of the root

x U f ( x L )−x L f ( x U )
xr=
f ( x L ) −f ( xU )

MATLAB MANUAL 81
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Find the absolute relative approximate error as

new old
x r −x r
|∈a|=| new
|×100
xr

where
new
xr = estimated root from present iteration
old
xr = estimated root from previous iteration

Compare the absolute relative approximate error

|∈a| with the pre-specified relative error tolerance
∈s . If |∈a|>∈s , then go to step 3, else stop the algorithm. Note one should also check

whether the number of iterations is more than the maximum number of iterations allowed. If so, one
needs to terminate the algorithm and notify the user about it.

Note that the false-position and bisection algorithms are quite similar. The only difference is the formula
used to calculate the new estimate of the root x r as shown in steps #2 and #4!

MATLAB CODE:

TASK:

A mass balance for a pollutant in a well-mixed lake can be written as:

dc
v =w−Qc−kv √ c
dt

Given the parameter values as follows:

v=1× 106 m3
Q=1 ×10 5 m 3 / year
w=1 ×106 gm/ year
k =0.25 m 0.5 / g0.5 / year
The equation takes the form:

25 √ c+ 10 c−100=0

Solve for the steady-state concentration, find the root using Regula-Falsi method in MATLAB.

PROGRAM:
clc
close all
clear all

f=@(x) 25*sqrt(x)+10*x-100 % Let steady-state concentration = x

MATLAB MANUAL 82
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

f(0)
f(5)
f(10)
fplot(f,[0,10])
xlabel ('Steady-state concentration, c');
ylabel ('f(c)');
title ('REGULA-FALSI METHOD');
grid on
gtext('xr');

xa=0; % Initial guess

xb=10; % Initial guess

fprintf ('i xa xb xr f(xr) RE\n');

for i=1:15
xr= ((xa*f(xb)-xb*f(xa))./(f(xb)-f(xa))); % Regula-Falsi formula
f_xr=25*sqrt(xr)+10*xr-100;
RE=abs((xr-xb)/xr)*100;

fprintf ('%6d %15.5f %15.5f %15.5f %15.5f %18.5f \n' ,

i,xa,xb,xr,f_xr,RE);

if f(xr)*f(xa)<0 % Root Bracketing

xb=xr;
xa=xa;

else
xa=xr;
xb=xb;

end %End of if-else Statement

end %End of For-Loop Statement
fprintf('\n The root of given equation = %f \n' , xr); % to print statement

OUTPUT:
f=
@(x)25*sqrt(x)+10*x-100

ans =
-100

ans =
5.9017

ans =
79.0569

i xa xb xr f(xr) RE

MATLAB MANUAL 83
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

1 0.00000 10.0000 5.58482 14.92869 79.05694

2 0.00000 5.58482 4.85937 3.70372 14.92869
3 0.00000 4.85937 4.68582 0.97516 3.70372
4 0.00000 4.68582 4.64057 0.26068 0.97516
5 0.00000 4.64057 4.62851 0.06997 0.26068
6 0.00000 4.62851 4.62527 0.01880 0.06997
7 0.00000 4.62527 4.62440 0.00505 0.01880
8 0.00000 4.62440 4.62417 0.00136 0.00505
9 0.00000 4.62417 4.62410 0.00037 0.00136
10 0.00000 4.62410 4.62409 0.00010 0.00037
11 0.00000 4.62409 4.62408 0.00003 0.00010
12 0.00000 4.62408 4.62408 0.00001 0.00003
13 0.00000 4.62408 4.62408 0.00000 0.00001
14 0.00000 4.62408 4.62408 0.00000 0.00000
15 0.00000 4.62408 4.62408 0.00000 0.00000

The root of given equation = 4.624081

GRAPH:

REGULA-FALSI METHOD

50 xr

0
f(c)

-50

-100
0 2 4 6 8 10
Steady-state concentration, c

LAB#11:

MATLAB MANUAL 84
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

CODE OF FIXED POINT ITERATION METHOD, SECANT

METHOD AND NEWTON RAPHSON METHOD FOR SOLUTION
OF TRANSCENDENTAL EQUATIONS IN MATLAB
The learning objectives are:

 How does Open Methods work?

 Background of Open Methods.
 Algorithm.
 Convergence sheet.

3. NEWTON-RAPHSON METHOD:
ALGORITHM:

The steps of the Newton-Raphson method to find the root of an equation f ( x )=0 are
'
1. Evaluate f ( x ) symbolically
2. Use an initial guess of the root,
x i , to estimate the new value of the root, x i+1 , as
f ( xi )
x i+1 = xi −
f ' ( xi )

3. Find the absolute relative approximate error

|∈a| as
x i+1 − x i
|∈a| = | |×10 0
x i+1

4. Compare the absolute relative approximate error with the pre-specified relative error
∈
tolerance, s . If
|∈ | ∈
a > s , then go to Step 2, else stop the algorithm. Also, check if the
number of iterations has exceeded the maximum number of iterations allowed. If so, one
needs to terminate the algorithm and notify the user.

MATLAB CODE:

TASK:

A mass balance for a pollutant in a well-mixed lake can be written as:

dc
v =w−Qc−kv √ c
dt

Given the parameter values as follows:

v=1× 106 m3
Q=1 ×10 5 m 3 / year

MATLAB MANUAL 85
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

w=1 ×106 gm/ year

k =0.25 m 0.5 / g0.5 / year
The equation takes the form:

25 √ c+ 10 c−100=0

Solve for the steady-state concentration, find the root using Newton-Raphson method in MATLAB.

PROGRAM:
%%
clc
clear all
close all
f=@(x) 25*sqrt(x)+10*x-100 % Let steady-state concentration = x
Df=@(x) 12.5*x^(-1/2)+10;
f(0)
f(5)
f(10)
fplot(f,[0,10])
xlabel ('Steady-state concentration, c');
ylabel ('f(c)');
title ('NEWTON RAPHSON METHOD');
grid on
gtext('xr');

xa=10; % Initial Guess

fprintf ('i xa xr f(xr) RE\n');

for i=1:5
xr=xa - (f(xa)/Df(xa));
f_xr=25*sqrt(xr)+10*xr-100;
Df_xr=12.5*xr^(-1/2)+10;
RE=abs((xr-xa)/xr)*100;
xa=xr;
fprintf ('%6d %15.5f %15.5f %15.5f %18.5f \n' , i,xa,xr,f_xr,RE);

end %End of For-Loop Statement

fprintf('The root of given equation = %f \n' , xr);

OUTPUT:

f=
@(x) 25*sqrt(x)+10*x-100
ans =
-100
ans =
5.9017
ans =
79.0569

MATLAB MANUAL 86
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

i xa xr f(xr) RE
1 4.33399 4.33399 -4.61447 130.73415

2 4.62232 4.62232 -0.02787 6.23770

3 4.62408 4.62408 -0.00000 0.03812

4 4.62408 4.62408 0.00000 0.00000

5 4.62408 4.62408 0.00000 0.00000

The root of given equation = 4.624081

GRAPH:

NEWTON RAPHSON METHOD

50
xr

0
f(c)

-50

-100
0 2 4 6 8 10
Steady-state concentration, c

4. SECANT METHOD:
ALGORITHM:
The Newton Raphson method of solving a non-linear equation of f (x)=0 is given by the
iterative formula:
f ’ ( x)=((f ( xi) – f (xi−1))/( xi – xi−1)).

MATLAB MANUAL 87
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

One of that’s method drawback is that we find derivative of that function, with availability of
symbolic manipulator like math, MAT-Lab and they are more convenient, laborious,
intractable if function derive as part of numerical scheme. In overcome that drawback th
function derivative is reduce to:
Put in above formula:
Xi +1=xi – ¿
The above eq. Is called the Secant Method, requires 2 initial guesses and no need to bracket
the root of the eq. It is an open method or may converge and faster than the Bi-Sectional
Method but slower than Newton Raphson Method.
The Newton-Raphson method of solving a nonlinear equation f (x )=0 is given by the iterative
formula

f ( xi )
x i+1 = xi −
f ' ( xi ) (1)One of the
drawbacks of the Newton-Raphson method is that you have to evaluate the derivative of the
function. With availability of symbolic manipulators such as Maple, MathCAD,
MATHEMATICA and MATLAB, this process has become more convenient. However, it still can
be a laborious process, and even intractable if the function is derived as part of a numerical
scheme. To overcome these drawbacks, the derivative of the function, f (x ) is approximated
as

' f ( x i )−f ( x i−1 )

f ( x i )=
x i −xi −1 (2)

Substituting Equation (2) in Equation (1) gives

f ( x i )(x i −x i−1 )
x i+1 =xi −
f ( x i )−f ( xi−1 ) (3)

The above equation is called the secant method. This method now requires two initial
guesses, but unlike the bisection method, the two initial guesses do not need to bracket the
root of the equation. The secant method is an open method and may or may not converge.
However, when secant method converges, it will typically converge faster than the bisection
method. However, since the derivative is approximated as given by Equation (2), it typically
converges slower than the Newton-Raphson method.

MATLAB CODE:
TASK

A mass balance for a pollutant in a well-mixed lake can be written as:

dc
v =w−Qc−kv √ c
dt

MATLAB MANUAL 88
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Given the parameter values as follows:

v=1× 106 m3
Q=1 ×10 5 m 3 / year
w=1 ×106 gm/ year
k =0.25 m 0.5 / g0.5 / year
The equation takes the form:

25 √ c+ 10 c−100=0

Solve for the steady-state concentration, find the root using Secant method in MATLAB.

PROGRAM:
%%
%%
clc
close all
clear all

f=@(x) 25*sqrt(x)+10*x-100 % Let steady-state concentration = x

f(0)
f(5)
f(10)
fplot(f,[0,10])
xlabel ('Steady-state concentration, c');
ylabel ('f(c)');
title ('Secant METHOD');
grid on
gtext('xr');

xa=0;
xb=10;
f_xa=@(x) 25*sqrt(xa)+10*xa-100;
f_xb=@(x) 25*sqrt(xb)+10*xb-100;

fprintf ('i xa xb xr f(xr) RE\n');

for i=1:8
xr=(xa*f(xb)-xb*f(xa))/(f(xb)-f(xa));
f_xr=25*sqrt(xr)+10*xr-100;
RE=abs((xr-xa)/xr)*100;

fprintf ('%6d %15.7f %15.7f %15.7f %15.7f %18.8f \n' ,

i,xa,xb,xr,f_xr,RE);
xa=xb;
f_xa=f_xb;

xb=xr;
f_xb=f_xr;

end
fprintf ('\n The root of given equation = %f \n' , xr);

MATLAB MANUAL 89
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

OUTPUT:

f=
@(x)25*sqrt(x)+10*x-100
ans =
-100

ans =
5.9017

ans =
79.0569

i xa xb xr f(xr) RE

1 0.0000000 10.0000000 5.5848156 14.9286921 100.000000

2 10.00000 5.5848156 4.5569858 -1.0623985 119.4433017

3 5.5848156 4.5569858 4.6252716 0.0188273 20.74567944

4 4.5569858 4.6252716 4.6240825 0.0000253 1.45102891

5 4.6252716 4.6240825 4.6240809 -0.0000000 0.02574892

6 4.6240825 4.6240809 4.6240809 -0.0000000 0.00003454

7 4.6240809 4.6240809 4.6240809 0.0000000 0.00000000

8 4.6240809 4.6240809 4.6240809 0.0000000 0.00000000

The root of given equation = 4.624081

GRAPH:

SECANT METHOD

50 xr

0
f(c)

-50

-100
0 2 4 6 8 10
Steady-state concentration, c

MATLAB MANUAL 90
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

5. FIXED POINT ITERATION METHOD:

ALGORITHM:
In numerical Analysis, this method is used to find the fixed point iteration, more specifically
given function defines on the real values and a point x b in the domain of f, the fixed point
iteration is:
x n+1¿ g ¿n¿ n=0,1,2,3.... .
which gives rise to sequence x0, x1, x2 , ....., which is hoped to convert to a point x, if f is
continuous, then one can prove that the obtained ‘x’ is a fixed point of f, i.e.
f (x)= x

More generally the function f can be defined on any matrix space with values in that same
space.

MATLAB CODE:

TASK:

A mass balance for a pollutant in a well-mixed lake can be written as:

dc
v =w−Qc−kv √ c
dt

Given the parameter values as follows:

v=1× 106 m3
Q=1 ×10 5 m 3 / year
w=1 ×106 gm/ year
k =0.25 m 0.5 / g0.5 / year
The equation takes the form:

25 √ c+ 10 c−100=0

Solve for the steady-state concentration, find the root using Fixed-point iteration method in MATLAB.

PROGRAM:

%%
%%
clc
close all
clear all

f=@(x) 25*sqrt(x)+10*x-100 % Let steady-state concentration = x

g=@(x) 10-2.5*sqrt(x);
f(0)
f(5)
f(10)
fplot(f,[0,10])

MATLAB MANUAL 91
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

xlabel('Steady-state concentration, c');

ylabel('f(c)');
title('FIXED-POINT ITERATION METHOD');
grid on
gtext('xr');

xa=0; % Initial guess

fprintf('i xa xr g_xr RE \n');

for i=1:45
xr=g(xa);
g_xr=10-2.5*sqrt(xr);
RE=abs((xr-xa)/xr)*100;
fprintf('\n%3d %15.7f %15.7f %15.7f %15.7f
%18.8f\n',i,xa,xr,g_xr,RE);
xa=xr;
end
fprintf('\n The root of given equation = %f \n',xr);

OUTPUT:

f=
@(x) 25*sqrt(x)+10*x-100

ans =
-100

ans =
5.9017

ans =
79.0569

i xa xr g_xr RE
1 0.0000000 10.0000000 2.0943058 100.0000000

2 10.0000000 2.0943058 6.3820708 377.4851773

3 2.0943058 6.3820708 3.6843098 67.1845407

4 6.3820708 3.6843098 5.2013610 73.2229686

5 3.6843098 5.2013610 4.2983769 29.1664276

6 5.2013610 4.2983769 4.8168682 21.0075604

7 4.2983769 4.8168682 4.5131588 10.7640749

8 4.8168682 4.5131588 4.6889509 6.7294187

MATLAB MANUAL 92
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

9 4.5131588 4.6889509 4.5865036 3.7490704

10 4.6889509 4.5865036 4.6459690 2.2336682

11 4.5865036 4.6459690 4.6113725 1.2799356

12 4.6459690 4.6113725 4.6314734 0.7502439

13 4.6113725 4.6314734 4.6197854 0.4340064

14 4.6314734 4.6197854 4.6265785 0.2529973

15 4.6197854 4.6265785 4.6226293 0.1468258

16 4.6265785 4.6226293 4.6249248 0.0854305

17 4.6226293 4.6249248 4.6235904 0.0496329

18 4.6249248 4.6235904 4.6243661 0.0288607

19 4.6235904 4.6243661 4.6239152 0.0167735

20 4.6243661 4.6239152 4.6241773 0.0097514

21 4.6239152 4.6241773 4.6240249 0.0056681

22 4.6241773 4.6240249 4.6241135 0.0032950

23 4.6240249 4.6241135 4.6240620 0.0019153

24 4.6241135 4.6240620 4.6240919 0.0011134

25 4.6240620 4.6240919 4.6240745 0.0006472

26 4.6240919 4.6240745 4.6240846 0.0003762

27 4.6240745 4.6240846 4.6240788 0.0002187

28 4.6240846 4.6240788 4.6240822 0.0001271

29 4.6240788 4.6240822 4.6240802 0.0000739

30 4.6240822 4.6240802 4.6240814 0.0000430

31 4.6240802 4.6240814 4.6240807 0.0000250

32 4.6240814 4.6240807 4.6240811 0.0000145

33 4.6240807 4.6240811 4.6240808 0.0000084

34 4.6240811 4.6240808 4.6240810 0.0000049

35 4.6240808 4.6240810 4.6240809 0.0000029

MATLAB MANUAL 93
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

36 4.6240810 4.6240809 4.6240809 0.0000017

37 4.6240809 4.6240809 4.6240809 0.0000010

38 4.6240809 4.6240809 4.6240809 0.0000006

39 4.6240809 4.6240809 4.6240809 0.0000003

40 4.6240809 4.6240809 4.6240809 0.0000002

41 4.6240809 4.6240809 4.6240809 0.0000001

42 4.6240809 4.6240809 4.6240809 0.0000001

43 4.6240809 4.6240809 4.6240809 0.0000000

44 4.6240809 4.6240809 4.6240809 0.0000000

45 4.6240809 4.6240809 4.6240809 0.0000000

The root of given equation = 4.624081

GRAPH:

FIXED-POINT ITERATION METHOD

50
xr

0
f(c)

-50

-100
0 2 4 6 8 10
Steady-state concentration, c

MATLAB MANUAL 94
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

LAB#12:
CODE OF GAUSS-JACOBI METHOD AND GAUSS-SEIDEL
METHOD FOR SOLUTION OF SYSTEM OF LINEAR
EQUATIONS IN MATLAB
The learning objectives are:

 How does Gauss-Jacobi and Gauss-Seidel mathods work?

 Background.
 Algorithms.
 Convergence Sheet

1. JACOBI’S METHOD:
ALGORITHM:
Given a general set of n equations and n unknown’s, we have:
a11x1 + a12x2 + a13x3 + ...................... + a1nxn = c1.
a21x1 + a22x2 + a23x3 + ..................... + a2nxn = c2.
............................................................................
............................................................................
............................................................................
an1x1 + an2x2 + an3x3 + ...................... + ann xn = cn.

MATLAB MANUAL 95
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

If the Diagonal elements are non-zero each eq. Is re written for the corresponding
unknown, that is, the 1st eq. Is rewritten with ‘x 1’ on the LHS and 2nd with ‘x2’ on LHS and
so on.
i.e.
x1 = (c1 - a12x2 - a13x3 - ...................... - a1nxn ) / a11.
x2 = (c2 - a21x1 - a23x3 - ...................... - a2nxn ) / a22.
x3 = (c3 – a31x1 – a32x2 - ...................... - a3nxn ) / a33.
...............................................................................
...............................................................................
xn = (cn – an1x1 – an2x2 - ...................... - an,n-1xn-1 ) / an-1.

Generally, iterative form of method is:

x1(i+1) = (c1 - a12x2i - a13x3 i - ...................... - a1nxn i ) / a11.
x2(i+1) = (c2 - a21x1i - a23x3i - ...................... - a2nxni) / a22.
x3 (i+1) = (c3 – a31x1i – a32x2i - ...................... - a3nxni) / a33.
...............................................................................
...............................................................................
xn(i+1) = (cn – an1x1i – an2x2i - ...................... - an,n-1xn-1i) / an-1.

Also find the absolute relative error:

R.E=|((present value-Previous value)/ present value) | * 100.

MATLAB CODE:
TASK:
Solve the given system of equations using Jacobi's Method on MATLAB?
3 x 1−0.1 x 2−0.2 x 3=7.85
0.1 x 1−7 x2 −0.3 x3 =−19.3
0.3 x 1−0.2 x2 −10 x 3=71.4
PROGRAM:

%%
clc
clear all
close all

a11=3; a12=-0.1; a13=-0.2;

a21=0.1; a22=7; a23=-0.3;
a31=0.3; a32=-0.2; a33=10;

b1 = 7.85; b2=-19.3; b3=71.4;

x(1)=0; % Initial Guess

y(1)=0; % Initial Guess
z(1)=0; % Initial Guess

fprintf('i x(i+1) y(i+1) z(i+1) REx REy REz \n');

for i=1:10
x(i+1) = (b1 - a12*y(i) - a13*z(i))/a11;

MATLAB MANUAL 96
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

y(i+1) = (b2 - a21*x(i) - a23*z(i))/a22;

z(i+1) = (b3 - a31*x(i) - a32*y(i))/a33;

REx = (abs((x(i+1)-x(i))/x(i+1)))*100;
REy = (abs((y(i+1)-y(i))/y(i+1)))*100;
REz = (abs((z(i+1)-z(i))/z(i+1)))*100;

fprintf ('%6d %15.5f %15.5f %15.5f %15.5f %15.5f %15.5f \n', i,

x(i+1) , y(i+1) , z(i+1) , REx , REy, REz);
end

fprintf('\n x = % f' , x(i+1));

fprintf('\n y = % f' , y(i+1));
fprintf('\n z = % f\n' , z(i+1));

% to plot system of equations

X = [-20:20];
Y = [-20:20];
[x,y] = meshgrid (X,Y);

z1 = (b1 - a11*x - a12*y)/a13;

z2 = (b2 - a21*x - a22*y)/a23;
z3 = (b3 - a31*x - a32*y)/a33;

mesh(x,y,z1)
hold on
mesh (x,y,z2)
hold on
mesh (x,y,z3)
hold off
xlabel('x-axis');
ylabel('y-axis');
zlabel('z-axis');
title ('Plotting system of linear equations');
grid on

% Shortcut method: x=inv(A)*B;

OUTPUT:

i x(i+1) y(i+1) z(i+1) REx REy REz

1 2.61667 -2.75714 7.14000 100.00000 100.00000 100.0000
2 3.00076 -2.48852 7.00636 12.79992 10.79431 1.90745
3 3.00081 -2.49974 7.00021 0.00148 0.44863 0.08786
4 3.00002 -2.50000 6.99998 0.02613 0.01057 0.00322
5 3.00000 -2.50000 7.00000 0.00079 0.00006 0.00026
6 3.00000 -2.50000 7.00000 0.00004 0.00004 0.00001
7 3.00000 -2.50000 7.00000 0.00000 0.00000 0.00000
8 3.00000 -2.50000 7.00000 0.00000 0.00000 0.00000
9 3.00000 -2.50000 7.00000 0.00000 0.00000 0.00000
10 3.00000 -2.50000 7.00000 0.00000 0.00000 0.00000

x = 3.000000

MATLAB MANUAL 97
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

y = -2.500000
z = 7.000000

GRAPH:

2. GUASS-SEIDEL
METHOD:
ALGORITHM:
Given a general set of n equations and n unknown’s, we have:
a11x1 + a12x2 + a13x3 + ...................... + a1nxn = c1.
a21x1 + a22x2 + a23x3 + ..................... + a2nxn = c2.
............................................................................
............................................................................
an1x1 + an2x2 + an3x3 + ...................... + annxn = cn.

If the Diagonal elements are non-zero each eq. Is re written for the corresponding
unknown, that is, the 1st eq. Is rewritten with ‘x 1’ on the LHS and 2nd with ‘x2’ on LHS and
so on.
i.e.
x1 = (c1 - a12x2 - a13x3 - ...................... - a1nxn ) / a11.
x2 = (c2 - a21x1 - a23x3 - ...................... - a2nxn ) / a22.
x3 = (c3 – a31x1 – a32x2 - ...................... - a3nxn ) / a33.
...............................................................................
...............................................................................
xn = (cn – an1x1 – an2x2 - ...................... - an,n-1xn-1 ) / an-1.

Iterative form of method is:

x1(i+1) = (c1 - a12x2i - a13x3 i - ...................... - a1nxn i ) / a11.
x2(i+1) = (c2 - a21x1i+1 - a23x3i - ...................... - a2nxni) / a22.
x3 (i+1) = (c3 – a31x1i+1 – a32x2i+1 - ...................... - a3nxni) / a33.
...............................................................................
...............................................................................
xn(i+1) = (cn – an1x1i+1 – an2x2i+1 - ...................... - an,n-1xn-1i+1) / an-1.

MATLAB MANUAL 98
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Also find the absolute relative error:

R.E=|((present value-Previous value)/ present value) | * 100.

MATLAB CODE:

TASK:
Solve the given system of equations using Gauss Seidel Method on MATLAB?
3 x 1−0.1 x 2−0.2 x 3=7.85
0.1 x 1−7 x2 −0.3 x3 =−19.3
0.3 x 1−0.2 x2 −10 x 3=71.4

PROGRAM:
%%
clc
clear all
close all

a11=3; a12=-0.1; a13=-0.2;

a21=0.1; a22=7; a23=-0.3;
a31=0.3; a32=-0.2; a33=10;

b1 = 7.85; b2=-19.3; b3=71.4;

x(1)=0;
y(1)=0;
z(1)=0;
fprintf('i x(i+1) y(i+1) z(i+1) REx REy REz \n');

for i=1:6
x(i+1) = (b1 - a12*y(i) - a13*z(i))/a11;
y(i+1) = (b2 - a21*x(i+1) - a23*z(i))/a22;
z(i+1) = (b3 - a31*x(i+1) - a32*y(i+1))/a33;

REx = (abs((x(i+1)-x(i))/x(i+1)))*100;
REy = (abs((y(i+1)-y(i))/y(i+1)))*100;
REz = (abs((z(i+1)-z(i))/z(i+1)))*100;

fprintf ('%6d %15.5f %15.5f %15.5f %15.5f %15.5f %15.5f \n', i,

x(i+1) , y(i+1) , z(i+1) , REx , REy, REz);
end

fprintf('\n x = % f' , x(i+1));

fprintf('\n y = % f' , y(i+1));
fprintf('\n z = % f \n' , z(i+1));

X = [0:20];
Y = [0:20];
[x,y] = meshgrid (X,Y);

z1 = (b1 - a11*x - a12*y)/a13;

z2 = (b2 - a21*x - a22*y)/a23;
z3 = (b3 - a31*x - a32*y)/a33;

MATLAB MANUAL 99
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

mesh(x,y,z1)
hold on
mesh (x,y,z2)
hold on
mesh (x,y,z3)
hold off
xlabel('x-axis');
ylabel('y-axis');
zlabel('z-axis');
title ('Plotting system of equations');
grid on

% Shortcut method: x=inv(A)*B;

OUTPUT:

i x(i+1) y(i+1) z(i+1) REx REy REz

1 2.61667 -2.79452 7.00561 100.00000 100.00000 100.00000
2 2.99056 -2.49962 7.00029 12.50235 11.79774 0.07598
3 3.00003 -2.49999 7.00000 0.31584 0.01453 0.00416
4 3.00000 -2.50000 7.00000 0.00105 0.00048 0.00001
5 3.00000 -2.50000 7.00000 0.00001 0.00000 0.00000
6 3.00000 -2.50000 7.00000 0.00000 0.00000 0.00000

x = 3.000000
y = -2.500000
z = 7.000000

GRAPH:

MATLAB MANUAL 100

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

LAB#13:
CODE OF LAGRANGE INTERPOLATION METHOD AND
NEWTON DIVIDED DIFFERENCE METHOD FOR
INTERPOLATION IN MATLAB

The learning objectives are:

 How does Lagrange and Newton Divided Difference Interpolation Formula work?
 Background.
 Algorithms.
 Convergence Sheet

1. LAGRANGE INTERPOLATION METHOD:

ALGORITHM:
Polynomial interpolation involves finding a polynomial of order ‘n’ that passes through the
‘n+1’ data points. One of the methods used to find this polynomial is called ‘Lagrange-
Interpolation Formula’, other methods includes Newton’s divided difference polynomials
method and the direct methods.
The Lagrange-Interpolation polynomial is given by:
P n ( x ) =∑ Li(x ). f ( x i).
Where n in Pn(x) stands for the nth order polynomial that approximates the function y=f(x)
given that n+1 data points as (x0,y0),( x1,y1), ( x2,y2),.................... ( xn-1,yn-1), ( xn,yn).
and

L i(x )=∏ (x−x j/ x i−x j).

Li(x) is a weighting function that includes a product of n-1 terms and with terms of i=1
omitted.

MATLAB CODE:

MATLAB MANUAL 101

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

TASK:
Illustrate the implementation of Lagrange Interpolation method using MATLAB code for the
given set of points:

T: 0 8 16 24 32 40
O(T ): 14.621 11.843 9.870 8.418 7.305 6.413

Estimate 'O' at T=27?

PROGRAM:
%%
clc
clear all
close all

% Declaring Table
x=[0,8,16,24,32,40];
y=[14.621,11.843,9.870,8.418,7.305,6.413];
n=length(x);

% Given unknown x whose value has to be interpolated.

xa=27;
ya=0; % initializing y

for i=1:n
L=1; % Initializing Lagranges L

for j=1:n
if (j~=i)
L=L*(xa-x(j))/(x(i)-x(j));
end

end

ya=ya+L*y(i);
end

% To fit the polynomial

fx=polyfit(x,y,5)

% To plot the curve for the given data.

xnew=[0,8,16,24,xa,32,40];
ynew=[14.621,11.843,9.870,8.418,ya,7.305,6.413];
fprintf('The value of y at x=27 is: %f',ya);
plot(xnew,ynew,'k-',xa,ya,'ro','linewidth',2)
xlabel('x-values');
ylabel('y-values');
title('Lagrange Interpolation');
gtext('Interpolating point');
grid on

% To find the Error

exact=7.986;
RE=(abs((exact-ya)/exact))*100;

MATLAB MANUAL 102

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

fprintf('\n Percentage Relative Error=%f \n',RE);

%%
% shortcut method
x=[12,13,14,15,16];
y=[5,6,9,11.5,11];
xa=15;
p=polyfit(x,y,3)
ya=polyval(p,xa)
plot(x,y,'k-',xa,ya,'ro')

OUTPUT:

fx =

-0.0000 0.0000 -0.0002 0.0095 -0.4135 14.6210

The value of y at x=27 is: 7.968239

Percentage Relative Error = 0.222402

GRAPH:

Lagrange Interpolation
16

14
y -v alues

12
Interpolating point
10

6
0 10 20 30 40
x-values

2. NEWTON DIVIDED DIFFERENCE INTERPOLATION METHOD:

ALGORITHM:

For Newton’s divided difference method let us revisit the quadratic polynomial interpolation formula
f 2 ( x )=b 0 +b1 ( x−x 0 )+b2 ( x−x 0 )(x −x1 )

where

MATLAB MANUAL 103

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

b0 =f ( x 0 )

f ( x 1 )−f ( x 0 )
b1 =
x1 −x 0

f ( x 2 )−f ( x 1 ) f ( x 1 )−f ( x 0 )
−
x 2 −x 1 x 1 −x 0
b2 =
x 2 −x 0

Note that 0
b, b1 , and b2 are finite divided differences. 0 b,
b1 , and b2 are the first, second,
and third finite divided differences, respectively. We denote the first divided difference by

f [ x 0 ]=f ( x 0 )

the second divided difference by

f (x 1 )−f ( x 0 )
f [ x 1 , x0 ]=
x 1−x 0

and the third divided difference by

f [ x 2 , x 1 ]−f [ x 1 , x 0 ]
f [ x 2 , x1 , x 0 ]=
x 2 −x 0

f ( x 2 )−f (x 1 ) f ( x 1 )−f ( x0 )
−
x 2 −x 1 x 1 −x 0
=
x 2 −x 0

where
f [ x ],
0 f [ x , x ],
1 0 and
f [ x 2 , x1 ,x 0 ] are called bracketed functions of their variables
enclosed in square brackets.

Rewriting,

f 2 ( x )=f [ x0 ]+f [ x 1 , x 0 ]( x−x 0 )+f [ x 2 , x1 , x 0 ](x− x0 )( x−x 1 )

This leads us to writing the general form of the Newton’s divided difference polynomial for n+1 data
x , y , x , y ,. .. . .. , ( x n−1 , y n−1 ) , ( x n , y n )
points, ( 0 0 ) ( 1 1 ) , as

f n ( x )=b0 +b1 ( x−x 0 )+. .. .+bn ( x−x 0 )( x−x 1 ).. .( x−x n−1 )

where

b0 =f [ x0 ]

MATLAB MANUAL 104

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

b1 =f [ x 1 , x 0 ]

b2 =f [ x 2 , x 1 , x 0 ] ……

bn−1=f [ x n−1 , x n−2 , . .. . ,x 0 ]

bn =f [ xn , x n−1 ,. ... , x 0 ]
th
where the definition of the m divided difference is

bm =f [ x m ,. ......., x 0 ]

f [ x m ,. .. . .. .. , x 1 ]−f [ xm−1 , . .. .. .. . , x0 ]
=
x m−x 0

From the above definition, it can be seen that the divided differences are calculated recursively.

MATLAB CODE:

TASK:
Illustrate the implementation of Lagrange Interpolation method using MATLAB code for the
given set of points:

T: 0 8 16 24 32 40
O(T ): 14.621 11.843 9.870 8.418 7.305 6.413

Estimate 'O' at T=27?

PROGRAM:
%%
% Newton Divided Difference
clc
clear all
close all
x=[0,8,16,24,32,40];
transposeofx=x';
y=[14.621,11.843,9.870,8.418,7.305,6.413];
xa=27; % point where y is evaluated
yo=14.621;
n=length(x);
D=zeros(n,n); % Matrix of zeros of nxn where n= number of elements in x
D(:,1)=y; % All rows and 1st column of D is replaced by y

for j=2:n
for i=j:n
fprintf(' Newtons Divided difference table:\n');
D(i,j)=(D(i,j-1)-D(i-1,j-1))/(x(i)-x(i-j+1))
end
end

for i=n:-1:2

MATLAB MANUAL 105

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

for j=n:-1:2
p=1;
if i==j
p=p*(xa-(x(i)))*D(i,j);
end
end
end

ya=yo+p;
fprintf(' ya=%f\n',ya);
exact=7.96820;
fprintf(' exact value=7.96820\n');
RE=abs((exact-ya)/exact)*100;
fprintf(' Percentage Relative Error=%f \n',RE);

OUTPUT:

Newtons Divided difference table:

D=
14.6210 0 0 0 0 0
11.8430 -0.3473 0 0 0 0
9.8700 -0.2466 0.0063 0 0 0
8.4180 -0.1815 0.0041 -0.0001 0 0
7.3050 -0.1391 0.0026 -0.0001 0.0000 0
6.4130 -0.1115 0.0017 -0.0000 0.0000 -0.0000

ya=8.023250
exact value=7.96820
Percentage Relative Error=0.690871

MATLAB MANUAL 106

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

LAB#14:
CODE OF TRAPEZOIDAL RULE, SIMPSON’S 1/3 RULE AND
SIMPSON’S 3/8 RULE FOR NUMERICAL INTEGRATION IN
MATLAB

The learning objectives are:

 How does Trapezoidal rule and Simpson rule work?

 Background.
 Algorithms.
 Convergence Sheet.

1. TRAPEZOIDAL RULE :
ALGORITHM:
Just like in multiple-segment trapezoidal rule, one can subdivide the interval [ a,b ] into n
segments and apply Simpson’s 1/3 rule repeatedly over every two segments. Note that n

needs to be even. Divide interval [ a,b ] into n equal segments, so that the segment width
is given by

b−a
h=
n .

Now
b xn

 f ( x)dx   f ( x)dx
a x0

where

x0  a xn  b

b x2 x4 xn  2 xn

 f ( x)dx   f ( x)dx   f ( x)dx  ......   f ( x)dx   f ( x)dx

a x0 x2 xn  4 xn  2

Apply Trapezoidal Rule over each interval, we get:

b n−1
( b−a )
(
∫ f ( x ) dx= 2 n f ( x 0 ) + f ( x n ) +2 ∑ f ( x i )
a i=1
)
MATLAB MANUAL 107
SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

MATLAB CODE:

Task:

Integrate the given integral using Trapezoidal Rule in MATLAB. Use 10 segments with
step-size of ‘h’ to approximate the integral.
0.8
I =∫ (0.2+25 x−200 x 2¿ +675 x 3−900 x 4 + 400 x5 )dx ¿
0

Program:
%%
% Trapezoidal rule
clc
clear all
close all
f=@(x) (0.2+25*x-200*x^2+675*x^3-900*x^4+400*x^5);
fplot(f,[0,1])
xlabel('x');
ylabel('f(x)');
title('Finding area using Trapezoidal rule');
gtext('Area under curve');
grid on

a=0;
b=0.8;
n=10;
h=(b-a)/n;

y=0; % Initializing the point for loop

x=a; % Initializing 'x' point

for i=a+h:h:b-h
x=x+h; % Increment in x-values
y=y+feval(f,i); % add up the area of each trapezoid from y1 to yn-1
I1=h*y;
end

I2=(h*(feval(f,a)+feval(f,b))/2); % Area of yo and yn

I=I1+I2; % Total Area

fprintf('\n I1=%f',I1);
fprintf('\n I2=%f',I2);
fprintf('\n Trapezoidal Area=%f',I);

OUTPUT:
I1=1.597763
I2=0.017280
Trapezoidal Area=1.615043
GRAPH:

MATLAB MANUAL 108

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Finding area using Trapezoidal rule

4

f(x)
1

Area under curve

0

-1

-2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x

2. SIMPSONS 1/3 RULE:

ALGORITHM:
Just like in multiple-segment trapezoidal rule, one can subdivide the interval [ a,b ] into n
segments and apply Simpson’s 1/3 rule repeatedly over every two segments. Note that n

needs to be even. Divide interval [ a,b ] into n equal segments, so that the segment width
is given by

b−a
h=
n .
b xn

 f ( x)dx   f ( x)dx
Now a x0

Where
x0  a xn  b

b x2 x4 xn  2 xn

 f ( x)dx   f ( x)dx   f ( x)dx  ......   f ( x)dx   f ( x)dx

a x0 x2 xn  4 xn  2

Apply Simpson’s 1/3rd Rule over each interval,

b
 f ( x0 )  4 f ( x1 )  f ( x 2 )   f ( x 2 )  4 f ( x3 )  f ( x 4 ) 
 f ( x )dx  ( x
a
2  x0 ) 
 6 

 ( x4  x2 ) 
 6   ...

 f ( x n  4 )  4 f ( x n 3 )  f ( x n  2 )   f ( x n  2 )  4 f ( x n 1 )  f ( x n ) 
 ( xn2  xn4 )   ( xn  xn2 )  
 6   6 

MATLAB MANUAL 109

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Since

xi  x i  2  2h i  2, 4, ..., n

b
 f ( x0 )  4 f ( x1 )  f ( x 2 )   f ( x 2 )  4 f ( x3 )  f ( x 4 ) 
 f ( x )dx  2h 
a
6   2 h  6   ...

 f ( x n  4 )  4 f ( x n 3 )  f ( x n  2 )   f ( x n  2 )  4 f ( x n 1 )  f ( x n ) 
 2h    2 h  
 6   6 

h
  f ( x0 )  4 f ( x1 )  f ( x3 )  ...  f ( xn1 )  2 f ( x2 )  f ( x4 )  ...  f ( xn2 )  f ( xn )
3

 n 1 n 2

h
 f ( x 0 )  4  f ( xi )  2  f ( xi )  f ( x n ) 
3  i 1 i 2

 i  odd i  even 

b  
ba  n 1 n 2

 f ( x )dx  f ( x 0 )  4  f ( xi )  2  f ( xi )  f ( x n ) 
3n  i 1 i 2

 
a
i odd i even

MATLAB CODE:

Task:

Integrate the given integral using Simpsons 1/3 Rule in MATLAB. Use 10 segments
with step-size of ‘h’ to approximate the integral.
0.8
I =∫ (0.2+25 x−200 x 2¿ +675 x 3−900 x 4 + 400 x5 )dx ¿
0

Program:
%%
% Simpsons 1/3 rule
clc
clear all
close all
f=@(x) (0.2+25*x-200*x^2+675*x^3-900*x^4+400*x^5);
fplot(f,[0,1])
xlabel('x');
ylabel('f(x)');
title('Finding area using Simpsons 1/3 rule');
gtext('Area under curve');
grid on

MATLAB MANUAL 110

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

a=0;
b=0.8;
n=10;
h=(b-a)/n;
y1=0; % Initializing the point for loop 1
y2=0; % Initializing the point for loop 2
x=a; % Initializing 'x' point

for i=a+h:2*h:b-h % for-loop for odd points x1, x3, x5, ..., xn-2
x=x+h; % Increment in x-values
y1=y1+feval(f,i); % add up the area of odd points
end
I1=(4*h*y1)/3;
fprintf('\n I1=%f',I1);

for j=a+2*h:2*h:b-h % for-loop for even points x2,x4,x6,...,xn-1

y2=y2+feval(f,j); % add up the area of even points
end
I2=(2*h*y2)/3;
fprintf('\n I2=%f',I2);

I3=(h*(feval(f,a)+feval(f,b))/3); % Area of yo and yn

fprintf('\n I3=%f',I3);

I=I1+I2+I3; % Total Area

fprintf('\n Area using Simpsons 1/3 Rule=%f \n',I);

OUTPUT:

I1=1.126803
I2=0.501774
I3=0.011520
Area using Simpsons 1/3 Rule=1.640096

GRAPH:

MATLAB MANUAL 111

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Finding area using Simpsons 1/3 rule

4

f(x) 1

0 Area under curve

-1

-2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x

LAB#15:

MATLAB MANUAL 112

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

CODE OF MODIFIED EULER’S METHOD AND RUNGE-KUTTA

METHOD FOR SOLUTION OF FIRST ORDER DIFFERENTIAL
EQUATIONS IN MATLAB
The learning objectives are:

 How does Modified Euler’s method and Runge-Kutta method of order-4 work?
 Background.
 Algorithms.
 Convergence Sheet

1. MODIFIED EULERS METHOD:

ALGORITHM:
The objective in numerical methods is, as always, to achieve the most accurate and reliable
result with least effort. For integrating the Initial value problem, the effort is usually measured
by the number of times the function f (t , y ) must be evaluated in stepping from a to b. As we
will see, a simple improvement doubles the number of function evaluations per step, but
yields a second order method-a winning strategy.

The algorithm for the method is as follows:

The Modified Euler’s method starts with an Euler step, giving a provisional value for y i+1 at
the next timet i+1.

y 'i+1 = y i+ hf (t i , y i) eq(1)

The step actually taken looks like an Euler step, but with ‘f’ replaced by the average of ‘f’ at
the starting point of the step and ‘f’ at the provisional point.

h
y i+1 = y i+ ¿ eq(2)
2

Substituting eq(1) in eq(2):

h
y i+1 = y i+ [f ( t i , y i ) + f (t i+1 , y i+ hf ( t i , y i ) )]
2

which is the Iterative formula for the Modified Euler’s method.

MATLAB CODE:

Task:

MATLAB MANUAL 113

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Illustrate the implementation of Modified Euler’s method in MATLAB to find the solution of
the given first order Differential equation:

dy
=4 e0.8 x −0.5 y
dx

With an Initial condition y (0)=2 using h=0.5 with the range of x=0 to x=1.

Program:
%%
% Modified Eulers method
clc
clear all
close all

f=@(x,y) 4*exp(0.8*x)-0.5*y;

h=0.5; % Step-size
x=0:h:1.5; % range of x variable
n=length(x); % Number of elements in x vector

x(1)=0; % xo
y(1)=2; % yo

fprintf('Modified Euler Method: \n'); % for display text

fprintf(' i x(i) y(i+1) \n'); % for display of headings

for i=1:n-1
y(i+1)=y(i)+h*feval(f,x(i)+h/2,y(i)+0.5*h*feval(f,x(i),y(i)));
fprintf('%3d %10.5f %17.10f\n',i,x(i),y(i+1));
end

fprintf('\n Solution of first order Differential Equation using Modified Euler

method is: %f\n',y(i+1));

OUTPUT:

Modified Euler Method:

i x(i) y(i+1)
1 0.00000 3.7553055163
2 0.50000 6.2051138610
3 1.00000 9.7279236287

Solution of first order Differential Equation using Modified Euler method is: 9.727924

2. RUNGE-KUTTA METHOD OF ORDER-4:

ALGORITHM:

MATLAB MANUAL 114

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Euler’s method is given by

y i +1 = y i + f ( x i , y i ) h
where
x 0=0 , y0= y ( x0 ) ,
h=xi+1−x i

To understand the Runge-Kutta 2nd order method, we need to derive Euler’s method from the
Taylor series.

dy 1 d2 y 2 1 d y
3
3
y i+1 = y i + |x , y ( xi +1−x i ) + 2
|x , y ( xi +1−x i ) + 3
|x , y ( xi +1−x i ) +.. .
dx i i 2 ! dx i i 3 ! dx i i
1 1
= y i +f ( x i , yi ) ( x i+1 −x i ) + f ' (x i , y i )( x i +1−x i ) 2 + f ''( x i , y i ) ( x i+1 −x i )3 + .. .
2! 3!
As you can see the first two terms of the Taylor series

y i +1 = y i + f ( x i , y i ) h

are Euler’s method and hence can be considered to be the Runge-Kutta 1st order method.

The true error in the approximation is given by

f ' ( xi , yi) 2 f ' ' ( xi , yi) 3

Et = h + h +. . .
2! 3!

So what would a 2nd order method formula look like. It would include one more term of the
Taylor series as follows.

1 '
y i+1 = y i +f ( x i , y i ) h+ f x , y h2
2 ! ( i i)
'
However, we already see the difficulty of having to find f ( x , y ) in the above method. What
Runge and Kutta did was write the 2nd order method as

y i +1 = y i + ( a1 k 1 +a2 k 2 ) h

where

k 1 =f ( xi , y i )

k 2 =f ( xi + p 1 h , yi + q11 k 1 h )

This form allows one to take advantage of the 2nd order method without having to calculate
f ' (x, y)
.

The formula described in this chapter was developed by Runge. This formula is same as

Simpson’s 1/3 rule, if f (x , y) were only a function of x . There are other versions of the

MATLAB MANUAL 115

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

4th order method just like there are several versions of the second order methods. The
formula developed by R.K. is:

h
y i+1 = y i + k + 2k 2 +2 k 3 + k 4 )
6( 1
where

k 1 =f ( xi , y i )

1 1
(
k 2 =f x i + h , y i + k 1
2 2 )
1 1
(
k 3 =f x i + h , y i + k 2
2 2 )
k 4 =f ( x i +h , y i + k 3 )

This formula is the same as the Simpson’s 3/8 rule, if f ( x , y ) is only a function of x . In
addition, the fourth order Runge-Kutta method is similar to the Improved Euler’s method in
that multiple estimates of the slope are developed in order to come up with an improved
average slope for the interval.

MATLAB CODE:
Task:

Use the classical 4th order Runge-Kutta method to find the solution of the given first order
differential equation in MATLAB.

dy
=4 e0.8 x −0.5 y
dx

With an initial condition y (0)=2 using h=0.5 with the range of x=0 to x=1

Program:
%%
clc
clear all
close all
h=0.5; % step-size
x=0:h:1;
n=length(x);
x(1)=0; % value of x in initial condition
y(1)=2; % value of y in initial condition
f=@(x,y) (4*exp(0.8*x))-0.5*y;
fprintf('Runge-Kutta method of order-4\n');
fprintf(' i k1 k2 k3 k4 k y(i+1)\n');
for i=1:n
k1=f(x(i),y(i));
k2=f(x(i)+0.5*h,y(i)+0.5*k1);

MATLAB MANUAL 116

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

k3=f(x(i)+0.5*h,y(i)+0.5*k2);
k4=f(x(i)+h,y(i)+k3);
k=(h*(k1+2*(k2+k3)+k4))/6;
y(i+1)=y(i)+k;
fprintf('%3d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f\n',i,k1,k2,k3,k4,k,y(i+1));
end

fprintf('\n k=%f\n',k);
fprintf('\nSolution of First order Ordinary Differential equation:%f\n',y(i+1));

% FOR EXACT SOLUTION

xspan=0:0.5:1; % range for x values
y(1)=2;
fprintf('\n Exact solution is: ');
[x,y]=ode45(f,xspan,y(1)) % Direct solution syntax
plot(x,y)
xlabel('x-axis');
ylabel('y-axis');
title('Exact solution');
grid on
axis tight

OUTPUT:

Runge-Kutta method of order-4

i k1 k2 k3 k4 k y(i+1)
1 3.00000 3.13561 3.10171 3.41644 1.57426 3.57426
2 4.18017 4.45630 4.38727 4.92140 2.23239 5.80665
3 5.99884 6.47009 6.35228 7.20100 3.23705 9.04370

k=3.237049

Solution of First order Ordinary Differential equation: 9.043699

Exact solution is:

x=
0
0.5000
1.0000

y=
2.0000
3.7515
6.1946

MATLAB MANUAL 117

SESSION: 2017-2021 MECHANICAL ENGINEERING DEPARTMENT

Exact solution

5.5

4.5
y-axis
4

3.5

2.5

2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x-axis

MATLAB MANUAL 118