MATLAB for Engineers and Scientists

INDEX
Section-I Section-II Section-III Section-IV Section-V Section-VI Introduction to Matlab Data Plotting in Matlab Writing Scripts & Function Files Programming in Matlab 3D Plots & Data Visualization Case Studies

Section-VII Applications Section-VIII Symbolic Computations

Introduction to MATLAB

Section-I

Introduction to MATLAB

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

1 of 24

Introduction to MATLAB

+ - / * ^ Clear, clc, cla % format diary

Basic algebraic operations Clear variables from the workspace Beginning of a comment line Start saving commands in a file named diary Element by element algebraic operators Create an array of N equally spaced values between x0 and x1 Create an array of values between x0 and x1 spaced by x. Compute Inverse of a square matrix Compute determinant of a square matrix Transpose of a Matrix Compute rank of Matrix Command for opening a File Asking help from Matlab Returns all the commands related to that keyword

.*

linspace (x0, x1, N) x :x : x1 inv( ) det() ' Rank() fopen help lookfor

./

.^

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

2 of 24

Introduction to MATLAB

Introduction to MATLAB
Matlab is an interactive system for doing numerical computations. A numerical analyst called Cleve Moler wrote the first version of Matlab in the 1970s. It has since evolved into a successful commercial software package. Matlab relieves you of a lot of the mundane tasks associated with solving problems numerically. This allows you to spend more time thinking, and encourages you to experiment. Matlab makes use of highly respected algorithms and hence you can be confident about your results. Powerful operations can be performed using just one or two commands. You can build up your own set of functions for a particular application. MATLAB is an acronym for MATrix LABoratory and as the name suggests, was initially developed to provide an environment to deal with matrix based computations. Since then the software has been diversified and now includes numerous collection of special functions. These include Toolboxes, Block sets Virtual world etc. MATLAB also provides an easyto-use interface to display your results graphically on screen as well as on hard copy. Thus MATLAB provides an integrated environment and is an ideal tool to solve problems of practical interest in an interactive fashion.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

3 of 24

Introduction to MATLAB

MATLAB as a Calculation Engine

The basic arithmetic operators are + - * / ^ and these are used in conjunction with brackets: ( ). The symbol ^ is used to get exponents (powers): 2^4=16. Type the commands shown following the MATLAB prompt: >> 2 + 3/4*5 ans = 5.7500 Is this calculation 2 + 3/(4*5) or 2 + (3/4)*5? Matlab works according to the priorities: 1. quantities in brackets, 2. powers 2 + 3^2 2 + 9 = 11, 3. * /, working left to right (3*4/5=12/5), 4. + -, working left to right (3+4-5=7-5), Thus, the earlier calculation was for 2 + (3/4)*5 by priority 3 mentioned above. MATLAB can be used as a calculator and can store the inputs or the outputs into variables. If a semicolon ; is placed at the end of a statement, the value of the variable is calculated but not displayed on the screen. To clear variables from the workspace, use the command clear.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

4 of 24

Introduction to MATLAB

To display the current value of a variable, double-click on its name in command window. To avoid retyping previously used command, use up and down arrows.

Numbers & Formats

Matlab recognizes several different kinds of numbers

The e" notation is used for very large or very small numbers:

All computations in MATLAB are done in double precision, which means about 15 significant figures.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

5 of 24

Introduction to MATLAB The format (how Matlab prints numbers) is controlled by the format command. Type help format for full list. Should you wish to switch back to the default format then format will suffice. The command format compact is also useful in that it suppresses blank lines in the output thus allowing more information to be displayed.

MATLAB Variables Names

Legal variable names: Begin with one of a-z or A-Z Have remaining characters chosen from a-z, A-Z, 0-9, or Have a maximum length of 31 characters Should not be the name of a built-in variable, built-in function, or user-defined function Examples: xxxxxxxxx pipeRadius widgets_per_baubble mySum mysum Note: mySum and mysum are different variables. MATLAB is case sensitive.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

6 of 24

Introduction to MATLAB

Activity Logging in MATLAB

diary diary off Start saving commands in a file named diary diary(filename) Start saving commands in the file filename Stop the diary

The diary command allows saving all command window inputs and outputs (except graphics) in a separate file, which can be opened with most text editors (e.g. Microsoft Word) and printed. Start the diary and save the file as mydiary. Compute the number of seconds in one year. Close the diary. Open the file mydiary using any text editor. In the Command Window:

MATLAB workspace will be saved into a file named diary when no file name is specified.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

7 of 24

Introduction to MATLAB

Vector Operations
Create the vector V = [1 2 3] . Display each of the elements. Multiply the first and second elements. Change the third element to 4 and display the new vector on the screen.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

8 of 24

Introduction to MATLAB

Element by Element Algebraic Operations

.* ./ .^ Element by element algebraic operations Define two vectors
x = [1 2 5 3]

and y = [4 1 0 1] .

Multiply x and y-- element by element. Compute the squares of the elements of x.

A period is needed for element-by-element operations with arrays in front of *, / , and ^ operators. A period is NOT needed in front of + and - .

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

9 of 24

Introduction to MATLAB

Creating Arrays with Equally Spaced Elements

linspace (x0, x1, N) x : x : x 1
Create an array of N equally spaced values between x0 and x1 Create an array of values between x0 and x1 spaced by x.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

10 of 24

Introduction to MATLAB

Evaluating Functions using Vector Operations

Evaluate the function y given below between x = 0 and x = :

x2 y = e sin 2 x + 3
x

To do so, first create a vector x containing equally spaced elements in the domain of the function. Then compute the corresponding vector y. Each element of y will be equal to the value of the function evaluated at the corresponding elements of x.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

11 of 24

Introduction to MATLAB

Matrices & Vectors

MATLAB can store several numerical values in a single variable. For most of the text, well refer to these variables as arrays. MATLAB, however, treats these variables as either matrices or arrays. When treated as arrays, MATLAB addresses these variables like other conventional programming languages, by performing element-by-element operations. When these variables are treated as matrices, MATLAB performs matrix operations with the variables. In MATLAB the default uses matrix algebra rather than array operations. Array operations are defined as element-by-element operations. Array addition and subtraction are the same as matrix addition and subtraction that is elementby-element. The matrix operations are ,/ , ^ ,+ and for matrix multiplication, division, power, addition, and subtraction. These operators perform matrix algebra if the operation is defined for the matrices being used with the operators. The array operators are .*, ./, .^,+ and -. Entering matrices into Matlab is the same as entering a vector, except each row of elements is separated by a semicolon (;) or a return:

In MATLAB, matrices and vectors are created using square bracket[ ]. Rows are separated by semicolon character(;).

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

12 of 24

Introduction to MATLAB

Entering matrices into Matlab

The elements in a row can also be separated by commas (,) rather than spaces

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

13 of 24

Introduction to MATLAB

Matrix Multiplication
Matrices are multiplied by using (*) operator.

A*B is not necessarily equal to B*A because matrix multiplication is not commutative.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

14 of 24

Introduction to MATLAB

Multiplying Matrix by a Scalar.

Creating Row Vector

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

15 of 24

Introduction to MATLAB

Creating Column Vector

Matrix Addition

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

16 of 24

Introduction to MATLAB

Matrix Subtraction

Addition and subtraction can be done between matrices and a scalar number as well.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

17 of 24

Introduction to MATLAB

Inverse of a Matrix
Inverse of a square matrix can be obtained using inv( ) function.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

18 of 24

Introduction to MATLAB

Determinant of a Matrix

Transpose of a Matrix

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

19 of 24

Introduction to MATLAB

Solving System of Linear Equation

5x +7y +3z = 1--------------- (01) -3x + 16y -2.5z = - 4 ----------------- (02) 15x 12y +3.9z = 16.8 --------------- (03)

The equation Ax=b has a solution if the rank of A and rank of the augmented matrix[A b]are both equal to 3.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

20 of 24

Introduction to MATLAB

Rank of A

Rank of [A b]

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

21 of 24

Introduction to MATLAB

Since the rank of A and ([A b]) are the same, we can solve the set of the equations.

X=A-1b

x=1.2921 y=0.2141 z=-1.3206

X=A\b (Left Div.)

x=1.2921 y=0.2141 z=-1.3206

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

22 of 24

Introduction to MATLAB

Opening, Reading and Closing Binary Data Files

A file containing Human Readable Text is also called ASCII text file. It uses a lot of dist space. If we need to create large data files that will only be read by computers, binary files are much more space-efficient. The numbers are written in their binary forms and are unintelligible if we look at them with the help of a text editor. Whether we are reading or writing a file in MATLAB, the very first thing is to open the file to be read or written. This is done with fopen command.

FID=fopen (filename, permissions)

Variable filename is a text string that contains the file name. If we are opening the file for reading it must already exist. If we are opening a file for writing and the file doesnt exist, a new file with the name contained in variable filename will be created. Variable
permission is a code that specifies how a file is opened. The valid codes are:

Permissions Action Open the file for reading. The file must exist in order to open it. An error r r+ w w+ a a+

occurs if the file does not exist. Open the file for reading and writing. The file must exist in order to open it. An error occurs if the file does not exist. Open the file for writing. If the file already exists, delete its contents. If the file does not exist, create a new file. Open the file for reading & writing. If the file already exists, delete its contents. If the file does not exist, create a new file. Open the file for writing. If the file already exists, append new data to the end of the file.. If the file does not exist, create a new file. Open the file for reading & writing. If the file already exists, append new data to the end of the file. If the file does not exist, create a new file.

More than one file can be opened at a time. Variable fid is a code that identifies the file just opened. Each file opened is identified by a unique fid code.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

23 of 24

Introduction to MATLAB

Seeking Help from MATLAB

If you need help with a specific MATLAB command, simply type help command name in the command window. Detailed information will be displayed with everything needed to use this command properly. >>help rank If you dont know the exact name of the command, type lookfor keyword, which returns all the commands related to that keyword. >> lookfor rank CHOLUPDATE Rank 1 update to Cholesky factorization. QRUPDATE Rank 1 update to QR factorization. RANK Matrix rank. SPRANK Structural rank. GFRANK Compute the rank of a matrix over a Galois field. XREGQR qr decomposition with rank chank GANGSTR Zero out 'small' entries subject to structural rank. RANKSUM Wilcoxon rank sum test that two populations are identical. SIGNRANK Wilcoxon signed rank test of equality of medians. TIEDRANK Compute the ranks of a sample, adjusting for ties RANK Matrix rank. RANK Symbolic matrix rank. RANK returns position of component in child list RANKCHECK Perform rank checking on design RANKCHECK Perform rank checking on design RANK Matrix rank. RANKCHECK Perform rank checking on design FRANK Frank matrix.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

24 of 24

Data Plotting in Matlab

Section-II

Data Plotting in Matlab

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

1 of 21

Data Plotting in Matlab

plot title xlabel, ylabel legend hold on / hold off subplot Grid on/off

Create a scatter/line plot in 2-D Add a title to the current figure Add x and y labels to the current figure Add a legend to the current figure Plot multiple curves in the same figure Make several plots in the same window sets scaling for the x- and y-axes on the current plot.

axis([XMIN XMAX YMIN YMAX])

zoom axis square, normal, tight, equal semilogx semilogy loglog Polar Help comet

Axis control commands of a plot. Compute determinant of a square matrix Transpose of a Matrix Compute rank of Matrix Command for opening a File Asking help from Matlab Returns all the commands related to that

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

2 of 21

Data Plotting in Matlab

Plotting in X & Y(2D Plots)

plot title xlabel, ylabel legend hold on, hold off Subplot(m,n,p) Grid on/off
The default is to plot solid lines. A solid white line is produced by plot(x, y, 'w-') The third argument is a string whose 1st character specifies the color (optional) and the second the line style. The options for colors and styles are: Create a scatter/line plot in 2-D Add a title to the current figure Add x and y labels to the current figure Add a legend to the current figure Plot multiple curves in the same figure Breaks the figure window into an m-by-n matrix of small windows, selects the pth window for the current plot

Line Styles & Colors in Matlab Plots

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

3 of 21

Data Plotting in Matlab

The displacement of an object falling under the force of gravity is given by:

y = gt2
where g is the acceleration of gravity [m/s2], y is the distance traveled [m], and t is time [s].The following table gives the values of g for three different planets: Planet Earth Jupiter Mercury i) g 9.80 m/s2 24.79 m/s2 3.72 m/s2

Plot the distance fallen as a function of time near the surface of the Earth between 0 and 10 seconds, with an increment of 1 second. Add a title and x and y labels.

% Define time values t = 0:1:10; % Define acceleration due to gravity g = 9.80; % Compute displacement y y = g * t.^2; % Plot plot(t,y,'b--') % Create labels grid on xlabel('Time [s]') ylabel('Distance fallen [m]') title('Distance fallen as a function of time on Earth')

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

4 of 21

Data Plotting in Matlab

You can define the color and the style of a curve when using the plot command. In the example below, b- - stands for blue dashed line. Other possible options are : for a dotted line, -. for a dash-dotted line, - for a solid line. The colors are designated by the first letter of the desired color (b,r,g,y k = black). If the style or the color is not specified, MATLAB uses its default settings. See help plot for additional information. You can zoom in and out, draw arrows, add titles, labels, legends, or write text directly in the figure using the menu bar.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

5 of 21

Data Plotting in Matlab

Plotting Multiple Functions on Same Figure in one Plot

A call to plot clears the graphics window before plotting the current graph. This is not convenient if we wish to add further graphics to the figure at some later stage. To stop the window being cleared:

>> plot(x,y,'w-'), hold on >> plot(x,z,'gx'), hold off

ii) Plot the displacement curve for both the Earth and Jupiter in the same figure and add a legend. [Code continued from previous one] % Both functions are to be plotted in the same figure hold on % Compute displacement for Jupiter g = 24.79; y = g * t.^2; % Plot displacement for Jupiter plot(t,y,'r-') % Create labels title('Distance fallen as a function of time') xlabel('Time [s]') ylabel('Distance fallen [m]') % Add a legend legend('Earth','Jupiter') % End of the hold command hold off

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

6 of 21

Data Plotting in Matlab

Multiple functions can be plotted on the same figure using the hold on and hold off commands. Several graphs may be drawn on the same figure as in plot(x,y,'w',x,cos(3*pi*x),'g--') A descriptive legend may be included with legend('Sin curve','Cos curve') The legend command will give a list of line-styles, as they appeared in the plot command, followed by a brief description. Matlab fits the legend in a suitable position, so as not to conceal the graphs whenever possible. The legend may be moved manually by dragging it with the mouse.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

7 of 21

Data Plotting in Matlab

Plotting Multiple Functions on Same Figure in Separate Plots

The three curves in the previous example can also be plotted on the same figure in three different plots using the subplot command. iii) Plot the curves for Earth and Jupiter using subplot.

%Compute displacement for Earth t = 0:1:10; g = 9.80; y = g*t.^2; % Plot subplot (2,2,1) plot (t,y,'b--') xlabel ('Time [s]') ylabel('Distance fallen [m]') title ('Distance fallen as a function of time on Earth') grid on % Define displacement for Jupiter g=24.79; y = g*t.^2; % Plot subplot(2,2,2) plot(t,y,'r-') xlabel('Time [s]') ylabel('Distance fallen [m]') title('Distance fallen as a function of time on Jupiter') grid on % Define displacement for Mercury g=3.72; y = g*t.^2; % Plot subplot(2,2,3) plot(t,y,'k-') xlabel('Time [s]') ylabel('Distance fallen [m]') title('Distance fallen as a function of time on Mercury') grid on

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

8 of 21

Data Plotting in Matlab

subplot(221) (or subplot(2,2,1)) specifies that the window should be split into a 2 x 2 array and we select the first sub-window for plotting the data.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

9 of 21

Data Plotting in Matlab

Zooming
We often need to zoom in" on some portion of a plot in order to see more detail. Clicking on the Zoom in" or Zoom out" button on the Figure window is simplest but one can also use the command >> zoom Pointing the mouse to the relevant position on the plot and clicking the left mouse button will zoom in by a factor of two. This may be repeated to any desired level. Clicking the right mouse button will zoom out by a factor of two.
Holding down the left mouse button and dragging the mouse will cause a rectangle to be outlined. Releasing the button causes the contents of the rectangle to fill the window. Zoom off turns off the zoom capability.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

10 of 21

Data Plotting in Matlab

Formatted Text on Plots

It is possible to change to format of text on plots so as to increase or decrease its size and also to typeset simple mathematical expressions (in LATEX form).

The graph produced from the above mentioned code is shown below.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

11 of 21

Data Plotting in Matlab

The salient features of these commands are:

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

12 of 21

Data Plotting in Matlab

Controlling the Plot Axes

Once a plot has been created in the graphics window you may wish to change the range of x and y values shown on the picture. clf N = 100; h = 1/N; x = 0:h:1; y = sin(3*pi*x); plot(x,y) axis([-0.5 1.5 -1.2 1.2]) grid The axis command has four parameters, the first two are the minimum and maximum values of x to use on the axis and the last two are the minimum and maximum values of y. Note the square brackets. The result of these commands is shown below.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

13 of 21

Data Plotting in Matlab

You can experiment with the commands axis equal, axis square, axis normal, axis tight in any order to see the plots. clf N = 100; h = 1/N; x = 0:h:1; y = sin(3*pi*x); subplot(321) plot(x,y) title(default) grid subplot(322) plot(x,y) axis([-0.5 1.5 -1.2 1.2]) title('axis([-0.5 1.5 -1.2 1.2])') grid subplot(323) plot(x,y) axis equal title('axis equal') grid subplot(324) plot(x,y) axis square grid title('axis square') subplot(325) plot(x,y) axis normal grid title('normal') subplot(326) plot(x,y) axis tight grid title('tight')

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

14 of 21

Data Plotting in Matlab

It is easy to graph a curve (f(t), g(t)) in 2D-space.

>> >> >> >>

t = (0:2*pi/100:2*pi)'; plot(cos(t),sin(t)) axis('square') grid

Note the use of the axis('square') command to make the circle appear round instead of elliptic).

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

15 of 21

Data Plotting in Matlab

2D Plot Types
semilogx
The semilogx command works just like the Plot command except that the x-axis is a log scale. >>omega=logspace(1,4,200) >>magnitude=20*log10(abs(1000./(i*omega+1000))) >>semilogx(omega,magnitude) >>axis([10,10000,-20 5]) >>grid on >>ylabel(Gain(dB)) >> xlabel(Frequency(rad/sec)); >>title(Bode Plot)

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

16 of 21

Data Plotting in Matlab

semilogy Plotting y=5x10

3x

x=linspace(0,100,200) y=5.*10.^(3.*x); semilogy(x,y) ylabel(y) xlabel(x) grid

Note that the y-axis is now a log scale that is factors of 10 are equally spaced on the axis.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

17 of 21

Data Plotting in Matlab

loglog Plotting y=5x10

3x

x=linspace(0,100,200) y=5.*10.^(3.*x); loglog(x,y) ylabel(y) xlabel(x) grid

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

18 of 21

Data Plotting in Matlab

polar
A polar plot is used for data with polar coordinates.Lets plot a linear spiral r=2 theta=linspace(0,8*pi,200) r=2*theta; polar(theta,r) grid

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

19 of 21

Data Plotting in Matlab

comet
Comet function displays a plot and shows its trajectory as an animation. theta=linspace(1,20*pi,1000) r=5*log10(theta); x=r.*cos(theta); y=r.*sin(theta); comet(x,y)

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

20 of 21

Data Plotting in Matlab

Plot Annotation Example

>> >> >> >> >> >> >> D=load('pdxTemp.dat');m=D(:,1);T=D(:,2:4); plot(m,T(:,1),ro,m,T(:,2),k+,m,T(:,3),b-); xlabel('Month'); ylabel('Temperature({}^\circF)'); title('Monthly average temperature at PDX'); axis([1 12 20 100]); legend('High','Low','Average',2);

Monthly average temperature at PDX 100 90 80 70 60 50 40 30 20 High Low Average

Temperature(F)

6 7 Month

10

11

12

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

21 of 21

Writing Script & Function Files

Section-III

Writing Script & Function Files

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

1 of 8

Writing Script & Function Files

Scripts & Functions

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

2 of 8

Writing Script & Function Files

Flow Control

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

3 of 8

Writing Script & Function Files

File Input/Output

MATLAB Readable File Formats

DATA FORMATS MAT - MATLAB workspace CSV - Comma separated numbers DAT - Formatted text DLM - Delimited text TAB - Tab separated text SPREADSHEET FORMATS XLS - Excel worksheet WK1 - Lotus 123 worksheet COMMAND load csvread importdata dlmread dlmread xlsread wk1read RETURNS Variables in file. Double array. Double array. Double array. Double array. Double array and cell array. Double array and cell array. Cell array of CDF records Primary or extension table data HDF or HDF-EOS data set MATLAB movie. Truecolor, grayscale or indexed image(s). Truecolor, grayscale or indexed image. Truecolor or indexed image(s). Truecolor or indexed image. Truecolor or grayscale image. Indexed image. Indexed image. Indexed image. Indexed image. Indexed image.

SCIENTIFIC DATA FORMATS CDF - Common Data Format cdfread FITS - Flexible Image Transport System fitsread HDF - Hierarchical Data Format hdfread MOVIE FORMATS AVI - Movie IMAGE FORMATS TIFF - TIFF image PNG - PNG image HDF - HDF image BMP - BMP image JPEG - JPEG image GIF - GIF image PCX - PCX image XWD - XWD image CUR - Cursor image ICO - Icon image
Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

aviread imread imread imread imread imread imread imread imread imread imread 4 of 8

Writing Script & Function Files

RAS - Sun raster image PBM - PBM image PGM - PGM image PPM - PPM image AUDIO FORMATS AU - NeXT/Sun sound SND - NeXT/Sun sound WAV - Microsoft Wave sound imread imread imread imread auread auread wavread Truecolor or indexed. Grayscale image. Grayscale image. Truecolor image. Sound data and sample rate. Sound data and sample rate. Sound data and sample rate.

Between the % sign and the conversion character(e, f or g), you can add a minus sign to specify left adjustment and a field width designator e.g. 6.2, where the first digit specifies the minimum field width and second digit specify the precision(i.e. number of digits to the right of the decimal point.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

5 of 8

Writing Script & Function Files

Example:

Example:

x = 0:.1:1; y = [x; exp(x)]; fid = fopen('exp.txt','w'); fprintf(fid,'%6.2f %12.8f\n',y); fclose(fid);

Open a new script (or .m) file using edit. Write Matlab commands to plot the function as niceplot.m.

y = x sin x , between 0 and 15 using 1000 points. Include labels and a title. Save your file

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

6 of 8

Writing Script & Function Files

In the command window, select the appropriate directory and execute the script by typing niceplot in the command window.

It is a good idea to always start your script with a clear command. This will help eliminate a lot of unexpected errors associated with using previously computed values when you execute the script multiple times To insert comments within a script, start each line with a % symbol. The comments will appear green. If you type help niceplot in the command window, the comments contained at the beginning of the file will be displayed. This is how all help information is stored and displayed for existing MATLAB functions

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

7 of 8

Writing Script & Function Files

Writing Functions in Matlab

When entering a command such as roots, plot, or step into Matlab what you are really doing is running an m-file with inputs and outputs that have been written to accomplish a specific task. These types of m-files are similar to subroutines in programming languages in that they have inputs (parameters which are passed to the m-file), outputs (values which are returned from the m-file), and a body of commands, which can contain local variables. Matlab calls these m-files functions. You can write your own functions using the function command. The new function must be given a filename with a '.m' extension. This file should be saved in the same directory as the Matlab software, or in a directory, which is contained in MATLABs search path. The first line of the file should contain the syntax for this function in the form:

function [output1,output2] = filename(input1,input2,input3)

A function can input or output as many variables as are needed. The next few lines contain the text that will appear when the help filename command is evoked. These lines are optional, but must be entered using % in front of each line in the same way that you include comments in an ordinary m-file. Finally, below the help text, the actual text of the function with all of the commands is included Functions can be rather tricky to write, and practice will be necessary to successfully write one that will achieve the desired goal. Below is a simple example of what the function, add.m, might look like.

function [var3] = add(var1,var2) %Add is a function that adds two numbers var3 = var1+var2;
If you save these three lines in a file called "add.m" in the Matlab directory, then you can use it by typing at the command line:
Obviously, most functions will be more complex than the one demonstrated here. This example just shows what the basic form looks like.
Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

>>y = add(3,8)

8 of 8

Programming in Matlab

Section-IV

Programming in Matlab

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

1 of 11

Programming in Matlab

Inline Function

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

2 of 11

Programming in Matlab

Function of function

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

3 of 11

Programming in Matlab

If elseif else statements

Command : if condition Statement elseif condition statement else statement end

a)

Suppose you are buying apples at a supermarket. If you buy more than 5 apples, you get a 10% discount. Write a program computing the amount owed given the price per apple (assume it to be $1) and the number of apples bought. Suppose now the store offers a super discount of 40% if you buy more than 20 apples. Modify the code in part a) to implement this new condition using the elseif statement. (a) (b)

b)

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

4 of 11

Programming in Matlab

for loops
Command : for variable = initial value : increment : final value statement statement end a) Compute 10!, i.e. the product of all the integers from 1 to 10. b) c) Compute the sum of all even numbers from 0 to 100. Define the following piecewise function between x = 3 and x = 5 using 1000 points and plot it.

y(x ) = x y ( x ) = sin x y ( x ) = ln x
PART_A

x<0 0 x 1 x >1
PART_B

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

5 of 11

Programming in Matlab

PART_C

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

6 of 11

Programming in Matlab

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

7 of 11

Programming in Matlab

While loops
while condition statement end abs (x) returns the absolute value of x sqrt (x) returns the square root of x There are some occasions when we want to repeat section of MATLAB code until some logical condition is satisfied, but we cannot tell in advance how many times we have to go around the loop. This we can do with a whileend construct. What is the greatest value of n that can be used in the sum

Command(s):

12 + 22 + + n2 and get a value less than 100?

CODE RESULT

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

8 of 11

Programming in Matlab

% This function computes the value of n factorial using for loop

function [product]=factorialforloop(n) product=1; for i=1:n product=product*i; end

% This function computes the value of n factorial using if statements

function [product]=factorialifstatement(n) product=1; i=1; while (1) product=product*i; i=i+1; if i>n break end end
% This function computes the value % of n factorial using while loop

function [product]=factorialwhileloop(n) product=1; i=1; while i<=n product=product*i; i=i+1; end

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

9 of 11

Programming in Matlab

Random Numbers, Histogram, Mean and Stem Plot

Command(s) : mean Compute the average of the components of a vector stem Create a stem plot Highway patrol uses radar to measure the speed of the cars on a freeway. The speed limit is 55 mph. The officers pull over all cars whose speed is greater than 80 mph and write a ticket for $150 each. If 8000 vehicles are observed on the highway each day, on the average, how much money does the highway patrol collect each day? Compute the average over one full year and assume A = 10 in the randn command. clear all close all % the average of the daily amount is done over a year for i = 1:365 % create a vector whose components are the speed of each car taking the freeway x = randn(8000,1)*10 + 55; % Initialize the amount collected each day S(i) = 0; % test the speed of each vehicle. If it's greater than 80, add 150 to S for k = 1:length(x) if x(k) > 80 S(i) = S(i) + 150; end end end % Compute the average of vector S containing the amount %collected each day during a year Savg = mean(S) % Plot the amount of money collected each day stem(S) xlabel('days') ylabel('amount')

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

10 of 11

Programming in Matlab

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

11 of 11

3D plots & DATA Visualization

Section-V

3D Plots & Data Visualization

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

1 of 15

3D plots & DATA Visualization

Commands
Creates X and Y arrays for 3-D plots Create a two-dimensional grid of coordinate values Create wireframe parametric surfaces specified by X, Y, and Z [mesh] Draws the mesh surface together with the contour plot beneath the plotted [meshc] surface Creates a shaded surface using Z for the color data as well as surface height [surf ] surfc(...) is the same as surf(...) except that a contour plot is drawn beneath [surfc] the surface. Create a 3D surface on the domain defined by meshgrid [surface] SHADING controls the color shading of SURFACE. SHADING FLAT sets the shading of the current graph to flat. [shading] SHADING INTERP sets the shading to interpolate. SHADING FACETED sets the shading to faceted, which is the default. Adds a Colorbar [colorbar] COLORMAP(MAP) sets the current figure's colormap to MAP. Default is jet. [colormap()] Others are copper, hsv, spring, Autumn, winter, summer etc. CONTOUR(Z) is a contour plot of matrix Z treating the values in Z as heights above a plane. A contour plot is the level curves of Z for some values V. The [contour()] values V are chosen automatically The same as CONTOUR(...) except that the contours are filled. [contourf()] Areas of the data at or above a given level are filled. A three-dimensional analogue of plot(). [meshgrid]
[plot3()]

[sphere()] [cylinder()] [surfnorm()]

sphere command. creates 3 n+1 by n+1 matrices that can then be used in the mesh or surf commands to create unit spheres. surfnorm(X,Y,Z) or surfnorm(Z) plots the surface with the normals emanating from it.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

2 of 15

3D plots & DATA Visualization

In order to create a graph of a surface in 3D-space (or a contour plot of a surface), it is necessary to evaluate the function on a regular rectangular grid. This can be done using the meshgrid command. First, create 1D vectors describing the grids in the x- and y-directions:

>> x = (0:2*pi/20:2*pi)'; >> y = (0:4*pi/40:4*pi)';

Next, ``spread'' these grids into two dimensions using meshgrid:

>> [X,Y] = meshgrid(x,y); >> whos Name Size X 41x21 Y 41x21 x 21x1 y 41x1

Bytes 6888 6888 168 328

Class double double double double

array array array array

Grand total is 1784 elements using 14272 bytes The effect of meshgrid is to create a vector X with the x-grid along each row, and a vector Y with the y-grid along each column. Then, using vectorized functions and/or operators, it is easy to evaluate a function z = f(x,y) of two variables on the rectangular grid:

>> z = cos(X).*cos(2*Y);
Having created the matrix containing the samples of the function, the surface can be graphed using either the mesh or the surf commands.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

3 of 15

3D plots & DATA Visualization

>> mesh(x,y,z)

Plot using the mesh command

>> surf(x,y,z)

Plot using the surf command The difference is that surf shades the surface, while mesh does not.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

4 of 15

3D plots & DATA Visualization

[x,y]=meshgrid(-2:0.1:2,-2:0.1:2); z=sin(x.^2+y.^2); surface(x,y,z) pause view(-45,45) pause shading('interp'); pause grid on; pause colorbar

surf(x,y,z)

view(-45,45)

shading('interp');

grid on; colorbar

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

5 of 15

Example-I:
Plot the function z ( x, y ) =

3D plots & DATA Visualization

sin x 2 + y 2

for x ranging between -3 and 3, and y x 2 + y 2 + 1016 ranging between 4 and 4. Use an increment of 0.1 in both directions. Include labels and a title. Once the figure appears on the screen, rotate it using one of the menu buttons. Use the command shading(interp) to remove the discrete mesh. Add a colorbar to the figure.

The 1016 in the denominator is needed to avoid division by zero when

x= y=0

and to get the correct value in the limit:

lim
x 0 y 0

sin x 2 + y 2 x2 + y 2

) =1

clear; close all; [x,y] = meshgrid(-3:0.1:3,-4:0.1:4); z = sin(x.^2+y.^2)./(x.^2+y.^2+10^-16);

plot3(x,y,z) grid on figure ; surface(x,y,z)

view(-45,45), grid

shading('interp')

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

6 of 15

3D plots & DATA Visualization

shading('interp') colorbar xlabel('x') ylabel('y') zlabel('z') title('Nice 3D plot')

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

7 of 15

Example-II

3D plots & DATA Visualization

Consider the following function : a)

z ( x, y ) = y cos ( x ) sin ( 2 y )

Plot the function z for x&y ranging between 0 and 2, with an increment of 0.1 in both directions using the meshc command. This command draws the mesh surface together with the contour plot beneath the plotted surface. Make a contour plot of the function above using the contour command. This will only display the level curves in a horizontal plane. Include labels and a title. Using the contour command, you can specify the number of levels you would like to be displayed. The default number is 10.

b)

close all ; clear % Define the grid [x,y]=meshgrid(0:0.1:2*pi); % Define the function of two variables z = y.*cos(x).*sin(2*y); % Plot both surface and level curves meshc(x,y,z) grid on xlabel('x') ylabel('y') zlabel('z') title('meshc command') % Draw 20 level curves figure contour(x,y,z,20) grid on xlabel('x') ylabel('y') title('contour command')

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

8 of 15

3D plots & DATA Visualization

3D Sphere
If you want to create a unit sphere, you can use the sphere command. It will create 3 n+1 by n+1 matrices that can then be used in the mesh or surf commands to create unit spheres. In the sphere command n is designated by the number inside the parentheses following the command sphere. Function surf is used to visualize data as shaded surface.

>>[X,Y,Z] = sphere(30); >> mesh(X,Y,Z);

>> surf(X,Y,Z)

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

9 of 15

3D plots & DATA Visualization

3D Cylinder
You can also create a unit cylinder by using the cylinder command. You must state a vector R(defines the radius of the cylinder along the z-axis) and number N(specifying number of points used to define circumference of the cylinder) in the parentheses following cylinder.

>>t = -2*pi:pi/100:2*pi; >>[X,Y,Z] = cylinder(cos(t),20); >>mesh(X,Y,Z)

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

10 of 15

Plotting the Hyperbolic Paraboloid z = y2-x2

mesh(xi,yi,zi)

3D plots & DATA Visualization

x=-1:0.05:1; y=x; [xi,yi]=meshgrid(x,y); zi=yi.^2-xi.^2;

Mesh Surface without the Contour Plot

meshc(xi,yi,zi)

Mesh Surface together with the Contour Plot

Continued..
Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

11 of 15

3D plots & DATA Visualization

x=-1:0.05:1; y=x; [xi,yi]=meshgrid(x,y); zi=yi.^2-xi.^2;

surfc(xi,yi,zi); colormap copper;

shading interp;

view([25,15,20]);

grid off;

title('Hyperbolic Paraboloid z=y2-x2')

MESHC(...) is the same as MESH(...) except that a contour plot is drawn beneath the mesh. Because CONTOUR does not handle irregularly spaced data, this routine only works for surfaces defined on a rectangular grid. The matrices or vectors X and Y define the axis limits only.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

12 of 15

3D plots & DATA Visualization

contour3(zi);

contourf(zi);

[c,h]=contour(zi,'k-');

clabel(c,h);

title('The level curve for z=y2-x2')

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

13 of 15

3D plots & DATA Visualization

Review of 3D plotting commands

The command surfc plots a surface together with the level lines beneath. Unlike the command surfc the command surf plots a surface only without the level curves. The command colormap paints the surface using user-supplied colors(in this case copper).Other color maps are : hsv -hue-saturation-value colormap hot -black-red-yellow-white color map gray -linear gray scale color map bone -gray scale-with tinge of blue color map pink white flag lines color cube vga jet -(default) prism cool autumn spring winter summer etc. The command shading controls the color shading used to paint the surface. Shading interp sets the shading interpolated Shading flat sets the shading of current graph to flat Shading faceted sets the shading to faceted, which is default In order to enhance the graph use contourf instead of contour. countourf plots filled contour lines while contour doesn't. Matrix c and graphic handle h are used as the input parameters for the function clabel, which adds height labels to the current contour plot. Command view is the 3-D graph viewpoint specification. It takes a threedimensional vector, which sets the view angle in Cartesian coordinates.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

14 of 15

3D plots & DATA Visualization

Plotting elliptic Paraboloid with Normals

% initialize memory and figure clear all, clf reset % vector of u-coords on parameter plane u = linspace(-1,1,21); % vector of v-coords on parameter plane v = linspace(-1,1,21); % grid of points on parameter plane [U,V]=meshgrid(u,v); % x-component of vector function for surface,s=f(u,v) SX = U; % y-component of s=f(u,v) SY = V; % z-component of s=f(u,v) SZ = U.^2 + V.^2;

surfnorm(SX,SY,SZ)
xlabel('x-axis'), ylabel('y-axis'), zlabel('z-axis') title('elliptic paraboloid with normals') axis equal

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

15 of 15

Case Studies

Section-VI
CASE STUDY-I CASE STUDY-II CASE STUDY-III CASE STUDY IV CASE STUDY IV -------------------------------------------------------Electric Circuit

Case Studies
Motion of a Vehicle 3D trajectory of a Particle Curve Fitting Random Numbers

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

1 of 13

Case Studies

CASE STUDY-I

Motion of a CAR

Starting from rest, a car (mass = 900 kg) accelerates under the influence of a force of 2000 Newtons until it reaches 13.89 m/s (50 km/h). What is the distance traveled ? How long does it take to travel this distance ?

Newtons law:

dv =F dt

(1)

This equation determines the velocity of the car as a function of time. To obtain the position of the car, we use the following Kinematic relation:

dx =v dt

(2)

Although this problem can be solved analytically, in this exercise we would like to come up with a numerical solution using MATLAB. We first need to discretize the equations, i.e. evaluate the position and velocity at equally spaced discrete points in time, since MATLAB doesnt know how to interpret continuous functions. The solution for v(t ) and x(t ) will then be represented by vectors containing the values of the velocity and the position of the car respectively at these discrete points in time separated by t . The time derivative for the velocity in (1) can be approximated by:

dv v (t + t ) v(t ) dt t
Equation (1) then becomes:

v(t + t ) v(t ) F = t m

v(t + t ) = v(t ) + t

F m

or, if we solve for the velocity at step n+1:

F , where the unknown is m Similarly, equation (2) becomes:

v(n + 1) = v(n ) + t

v(n + 1) .

x(t + t ) x(t ) = v(t ) t

x(t + t ) = x(t ) + t v(t )

and the position at step n+1 will given by:

x(n + 1) = x(n ) + t v(n )

2 of 13

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

Case Studies
Solve the above equations for position and velocity and plot the position as a function of time. Display the distance traveled and the time it takes to reach 50 km/h. % Force applied to the car in Newtons F = 2000; % Mass of the car in kg m = 900; % Time step in seconds deltat = 0.1; % Initialization: at t=0, the car is at rest t(1) = 0; v(1) = 0; x(1) = 0; n = 1; % Iterate until the velocity reaches 50 km/h = 13.89 m/s. %For that, we need a WHILE loop which tests the value of %the velocity at each step.If it's greater than 13.89 m/s, we exit the loop. while v ( n ) < 13.89 % Advance time t(n+1) = t(n) + deltat; % Compute the new velocity v(n+1) = v(n) + F/m*deltat; % Compute the new position x(n+1) = x(n) + v(n)*deltat; % Increment n n = n+1; end % Plot the position as a function of time plot(t,x) grid on title('position of the car as a function of time') xlabel('time [s]') ylabel('position [m]') % Display the position reached by the car (last element of the x vector) x(length(x)) % Time needed to go from 0 to 50 km/h (last element of the t vector) t(length(t))

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

3 of 13

Case Studies

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

4 of 13

Case Studies

CASE STUDY-II

Electric Circuit

A circuit has the following current-voltage characteristic (i.e. a curve giving the output current in Amps as a function of the input voltage V in Volts:
I I 10 V + _
Saturation Linear regime Cut-off regime

1 1 5 V

Given the input voltage V (t ) , write a script that computes the output current I (t ) . Test your code for the following input: V (t ) = 3 + 3 sin (t ) for t 0 Use the time step t = 0.1 . Plot V (t ) and I (t ) on the same set of axes for t ranging between 0 and 2 seconds. clear close % Define the input voltage as a function of t t = 0:0.1:2*pi; V = 3+3*sin(3*t); % Test each component of V and compute the current output for i=1:length(V) if V(i)<1 I(i) = 10; elseif V(i)>1 & V(i)<5 I(i) = -9/4*V(i)+49/4; else I(i) = 1; end end % Plot the results plot(t,V,'-',t,I,'b-.') xlabel('Time') legend('Input Voltage','Output Intensity') title('Input and output of a circuit') grid on

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

5 of 13

Case Studies

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

6 of 13

Case Studies

CASE STUDY-III
Command : plot3 Plot a curve in 3D

3D trajectory of a Particle

The motion of a particle in three dimensions is given by the following parametric equations:
x ( t ) = cos ( t ) y ( t ) = 2sin ( t ) z (t ) = 0

for

0 t 2

a) b) c) vector d) e)

plot the particles trajectory in 3D using the plot3 function

Vx , Vy and the magnitude of the velocity vector compute the acceleration components Ax , Ay and the magnitude of the acceleration
compute the velocity components

r T compute the components of the unit tangent vector = Tx

Ty

Tz

compute the angle between the velocity vector and the acceleration vector clear % time vector dt = 0.1; t = 0:dt:2*pi; % x, y and z vectors x = cos(t); y = 2*sin(t); z = 0*t; a)% plot trajectory in 3D plot3(x,y,z) grid on xlabel('x') ylabel('y') zlabel('z') b)% Computation of the velocity components % Vx = dx/dt Vx = gradient(x,dt); % Vy = dy/dt Vy = gradient(y,dt);

% Speed speed = sqrt(Vx.^2+Vy.^2);

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

7 of 13

Case Studies
c)% Computation of the acceleration components % Ax = dVx/dt Ax = gradient(Vx,dt); % Ay = dVy/dt Ay = gradient(Vy,dt); % Acceleration acceleration = sqrt(Ax.^2+Ay.^2); d)% Computation of the unit tangent vector components % Tx = Vx/|V| Tx = Vx./speed; % Ty = Vy/|V| Ty = Vy./speed; e)% Computation of the angle between the velocity and acceleration vectors % cos(theta) = V.A/(|V||A|) theta = acos((Vx.*Ax+Vy.*Ay)./speed./acceleration);

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

8 of 13

Case Studies

CASE STUDY -IV

Command : polyfit

Curve Fitting

find a polynomial that best fits a set of data

Reynolds number is a non-dimensional parameter frequently used in fluid mechanics:

Re =

VL

where V is the characteristic velocity of the flow is the Kinematic viscosity L is the characteristic length scale of the flow (such as the length or the diameter of a body) The Reynolds number indicates the relative importance of the viscous effects compared to the inertia effects. If the Reynolds number is low, viscous effects are predominant and if its large, viscosity can be ignored. Lets consider the force exerted on a sphere moving in the air. This drag force depends on the viscosity of the fluid, the velocity of the sphere and its dimension. Therefore, it is directly related to the Reynolds number. The drag force can be expressed as follows: 1 FD = CD V 2 r 2 2 where is the density of air r is the radius of the sphere V is the velocity of the sphere CD is drag coefficient (non-dimensional)
In 1851 Stokes had demonstrated that the drag coefficient at low Reynolds numbers was inversely proportional to Re:

CD =

24 Re

(low

Re )

--------

(1)

A group of students have measured this coefficient for a sphere of radius r in a wind tunnel. Here are their results. Confirm that the data satisfies the Stokes law.
Re CD

0.05875 492

0.1003 256.2

0.1585 169.8

0.4786 58.88

3.020 10.86

Taking the logarithm of equation (1), we obtain:

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

9 of 13

Case Studies
CD = 24 Re ln ( CD ) = ln ( 24 ) ln ( Re ) ------- (2)

To check whether the Stokes law is satisfied, we can plot ln ( CD ) versus ln ( Re ) and see whether the data points fall along a straight line. Follow these steps:

Create two vectors containing the logarithm of the experimental data

ln ( CD )

ln ( Re ) and

Find the line that fits best these points. Display its coefficients on the screen and compare them with the theoretical result. In the same figure plot the experimental data, the best-fit line, and the theoretical line. Add a legend.

Hint: The polyfit command finds a polynomial of a specified order that best fits the given data. In this example, we want to see if the experimental data falls along a straight line. We therefore choose 1 as the degree of the polynomial. The polyfit command returns a vector containing the coefficients of the polynomial:

[ n

L 1 0 ] corresponding to P ( X ) = n X n + ... + 1 X 1 + 0

To suppress the line connecting the experimental data points, select the big arrow in the figures menu and double-click on the line. A property editor window will appear. Change the line style to no line and choose circles for the marker style. close all clear all % Experimental data Re = [0.05875 0.1003 0.1585 0.4786 3.020 7.015]; Cd = [492 256.2 169.8 58.88 10.86 5.5623]; % Take the logarithm of these data lnRe = log(Re); lnCd = log(Cd); % Construct the line that fits best data. As we want a % line, the degree of the %polynomial is set to 1. P = polyfit(lnRe,lnCd,1) % Define the fitting line, so it can be plotted lnRefit = linspace(log(0.05875),log(7.015),20); lnCdfit = P(1)*lnRefit + P(2); % Define the theoretical line : ln(Cd) = - ln(Re) + ln(24)
Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

10 of 13

Case Studies
lnRetheo = linspace(log(0.05875),log(7.015),20); lnCdtheo = log(24) - lnRetheo; plot(lnRe,lnCd) hold on plot(lnRefit,lnCdfit,'--') plot(lnRetheo,lnCdtheo,'-.') xlabel('ln(Re)') ylabel('ln(Cd)') title('drag coefficient on a sphere as a function of the Reynolds number') LEGEND('Experimental Data','Best Fit','Theoretical Fit') grid on

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

11 of 13

Case Studies

CASE STUDY-V
Command(s) :
randn hist

Random Numbers

generate a vector whose components are chosen at random plot a histogram

At the assembly line of a big aircraft company, bolts are stored in boxes containing 1000 parts each. The bolts are supposed to have a nominal diameter of 5 millimeters. But due to the uncertainty in the manufacturing process, some turn out to be larger and some turn out to be smaller than nominal. a) Create a random vector whose components are the diameters of each bolt in a box using the randn command. The syntax of the randn command is as follows: randn ( M , N ) * A + B, where: - M and N determine the size of the output vector, i.e. the number of random values generated. M corresponds to the number of rows and N to the number of columns. In our example, we want random diameters for 1000 bolts, so the output will be a vector with M = 1000 and N = 1 - A (called standard deviation) is related to the uncertainty in the diameter of the bolts. Here, A = 0.25 mm. This means that about 67% of the bolts have diameters within 0.25 mm of the nominal diameter. - B is the nominal diameter of the bolts. Here, B = 5 mm b) Plot a histogram showing the number of bolts as a function of their diameter % PART-a) x = randn(1000,1)*0.25 + 5; % PART-b) hist(x,30) title('Bolts contained in a box') xlabel('Diameter in mm') ylabel('Number of bolts') grid on Hint : The hist command takes the vector generated by the randn command as input and counts the number of bolts having a diameter within a certain range (a bin). The number of bins is defined by the second parameter.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

12 of 13

Case Studies

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

13 of 13

Applications

Section-VII

Applications

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

1 of 16

Applications

Polynomials
Polynomials in MATLAB are represented as row arrays. The elements of the row array represent the coefficients of the polynomial. Many computational activities in engineering and science require a polynomial. Suppose we have a function y=x2+2x+1. In MATLAB this polynomial is represented as [1 2 1]. A polynomial of an arbitrary order is represented by

a1 xn + a2xn-1 + a3xn-2 +. an x + an+1

In MATLAB, this polynomial is represented as of the row array is n+1.

[a1 a2 a3 +.+an an+1 ]. The length

Evaluating Polynomials

How do we can evaluate a polynomial for an arbitrary value of x (say x=5). MATLAB provides the polyval function for evaluating polynomials.

A polynomial can be evaluated at all the desired individual points at the same time by creating a row array for x. Polyval returns an array of results.

In this case Polyval command returns an array of results.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

2 of 16

Applications

Plotting Polynomials

After knowing how to represent polynomials in MATLAB, we can easily plot them.

>>poly=[1 0 1] %y=(x+1)(x-1)=x2-1 >>x=linspace(-2,2,100); >>y=polyval(poly,x); >>plot(x,y) >>grid

The MATLAB function conv multiplies two polynomials and returns the resultant polynomial as a row array. Suppose, we want to multiply the following three simple 1st order polynomials to get the resultant 3rd order polynomial.

Multiplying Polynomials

(x-1), (x+1), (x-2)

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

3 of 16

Applications

Roots of Polynomials

Consider the polynomial y=f(x)=x3-2x2-x+2. MATLAB provides a function called roots, Which finds the solution x to the equation y=f(x)=0 where f(x) is a polynomial.

Complex Roots of Polynomials

The roots function can also find the complex roots of a polynomial. Consider solving y=f(x)=x2+x+1 .

This result says that the solution to y=f(x)=x2+x+1 is

&

x= -0.5000 + 0.8660i x= -0.5000 - 0.8660i

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

4 of 16

Applications

Polynomial Curve fitting

S.No. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 Input (x) 0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 Output (y) 0.05 0.1 0.08 1.15 1.13 .95 2.1 2.2 2.05 3.1 3.05 2.95 4.2 4.2 4.05 5.6 4.9 5.1 6.1 6.03 6.12

In many engineering applications we need to find a polynomial to represent the given data set. We can use this polynomial in modelling the system under consideration. The MATLAB function polyfit, fits an nth order polynomial to a given data set.Consider the data set given on left. We will use polyfit to fit a curve to these data. The equation of the fitted curve will be verified graphically.

The equation representing the data is y=1.0107x+0.0731.To verify the equation with the given data set enter the following commands.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

5 of 16

Applications

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

6 of 16

Applications

Roots of Equation

Consider the problem of solving x +3x+5= sin(x). We are required to know the value of x 3 that makes this equation true. The equation can be re stated as x +3x+5 -sin (x)=0. The MATLAB function fzero finds the roots of an equation.. There appears to be a root close to x=0.(How?)

At x=-1.28273076924844, value of function y is such a small number that we can call it zero..

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

7 of 16

Applications

Numerical Integration of a function

Numerical Integration is most useful for finding integrals of functions that cant be integrated mathematically.

Example-I: Integrate x2+2x+1 from 0 to 1.

function y=NUMINT(x) y=x.^2+2*x+1-log(x) Or function y=NUMINT(x) %y=x.^2+2*x+1-log(x) y=x.*(x+2)+1

>> Integ_y=quad8('NUMINT',0,1)

y= Columns 1 through 2 1.00000000000000 1.05793318891372 Columns 3 through 4 1.19192179527508 1.39021391475144 Columns 5 through 6 1.62917960675006 1.90939445477740 Columns 7 through 8 2.25000000000000 2.61854405376385 Columns 9 through 10 2.97082039324994 3.31577394178878 Columns 11 through 12 3.64141153805826 3.88658043600016 Column 13 4.00000000000000 y= Columns 1 through 2 1.19192179527508 1.62917960675006 Columns 3 through 4 2.25000000000000 2.97082039324994 Column 5 3.64141153805826 Integ_y = 2.33333333333333

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

8 of 16

Applications

Example-II: Integrate e
function y=NUMINT(x) y=exp(-0.5*x.*x);

(-x^2/2)

from -5 to 5.

Trapz

Computes the integral of a function over a given domain using the trapezoidal method Compute the integral of the function f ( x ) = x + x + 1 for 0 trapezoidal method. Compare the result to the exact analytical solution. clear close
2

x 1

using the

x = linspace(0,1,100)

y = x .^2+x+1; I = trapz(x,y)

Analytical result : I = 11/6 = 1.8333

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

9 of 16

Applications

Solving Ordinary Differential Equations

In Matlab, there are a set of ODE solvers, namely ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb.General information regarding Matlab ODE solvers are given below: ode45 is based on an explicit Runge-Kutta (4,5) formula, the Dormand-Prince pair. It is a one-step solver in computing y(t n ), it needs only the solution at the immediately preceding time point, y(t n1 ). In general, ode45 is the best function to apply as a first try for most problems. ode23 is an implementation of an explicit Runge-Kutta (2,3) pair of Bogacki and Shampine. It may be more efficient than ode45 at crude tolerances and in the presence of moderate stiffness. Like ode45, ode23 is a one-step solver. ode113 is a variable order Adams-Bashforth-Moulton PECE solver. It may be more efficient than ode45 at stringent tolerances and when the ODE file function is particularly expensive to evaluate. Ode113 is a multistep solver it normally needs the solutions at several preceding time points to compute the current solution. The above algorithms are intended to solve non-stiff systems. If they appear to be unduly slow, try using one of the stiff solvers below. ode15s is a variable order solver based on the numerical differentiation formulas, NDFs. Optionally, it uses the backward differentiation formulas, BDFs (also known as Gears method) that are usually less efficient. Like ode113, ode15s is a multistep solver. If you suspect that a problem is stiff or if ode45 has failed or was very inefficient, try ode15s. ode23s is based on a modified Rosenbrock formula of order 2. Because it is a one-step solver, it may be more efficient than ode15s at crude tolerances. It can solve some kinds of stiff problems for which ode15s is not effective. ode23t is an implementation of the trapezoidal rule using a free interpolant. Use this solver if the problem is only moderately stiff and you need a solution without numerical damping. ode23tb is an implementation of TR-BDF2, an implicit Runge-Kutta formula with a first stage that is a trapezoidal rule step and a second stage that is a backward differentiation formula of order two. By construction, the same iteration matrix is used in evaluating both stages. Like ode23s, this solver may be more efficient than ode15s at crude tolerances.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

10 of 16

Applications
Stiffness is a special problem in the solution of ODEs. A stiff system is one which involves rapidly changing components together with slowly changing ones. The rapidly changing component is actually transient which dies out quickly, after which the system becomes dominated by the slowly varying components. Both individual and systems of ODEs can be stiff. Although the transient phenomena exist for only a short part of the integration interval, they can dictate the time step for the entire solution. An example of single stiff ODE is given below:

dy/dt = -1000y + 3000 2000e-t

for y(0)=0 the analytical solution can be developed as

y = 3 0.998e-1000t 2.002e-t
3 3.053 .10

3 2.953 .10 0 10 20 30 40 50 60 70 80 90 100 0 t 100

The above figure shows that the solution is initially dominated by the fast exponential term(e-1000t ). After a very short period t<0.005, this transient dies out and the solution becomes dictated by the slow exponential (002e-t).
Syntax:
Where

[T,Y] = ode45 ('Filename',tspan,y0)

Name of the file containing the ODE(s) to be solved A vector specifying the interval of integration [t0 tfinal]. To obtain solutions at specific times (all increasing or all decreasing), use tspan = [t0,t1, ..., tfinal]. A vector of initial conditions.

Filename tspan

y0

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

11 of 16

Applications

RKCOMP.m
tspan=[0 4] y0=0.5 [x,num_y]=ode23('g1',tspan,y0); an1_y=x.^3-7.5; subplot(2,2,1) plot(x,num_y,x,an1_y,'o') title('Solution using ode23 ') xlabel('x'); ylabel('dy/dx=x.^3-7.5'); grid on; [x,num_y]=ode45('g1',tspan,y0); an1_y=x.^3-7.5; subplot(2,2,2) plot(x,num_y,x,an1_y,'o') title('Solution using ode45'); xlabel('x'), ylabel('dy/dx=x.^3-7.5'); grid on; [x,num_y]=ode23t('g1',tspan,y0); an1_y=x.^3-7.5;
subplot(2,2,3)

plot(x,num_y,x,an1_y,'o') title('Solution solution using ode23t') xlabel('x'); ylabel('dy/dx=x.^3-7.5'); grid on; [x,num_y]=ode15s('g1',tspan,y0); an1_y=x.^3-7.5; subplot(2,2,4) plot(x,num_y,x,an1_y,'o') title('Solution using ode15s') xlabel('x') ylabel('dy/dx=x.^3-7.5') grid

g1.m function dy = g1(x,y) dy = 3*x^2;

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

12 of 16

Applications
Solution using ode23 80 dy/dx=x. -7.5 dy/dx=x. -7.5 60 40 20 0 -20 0 2 3 4 x Solution solution using ode23t 80 dy/dx=x. -7.5 60 40 20 0 -20 0 1 2 x 3 4 0 1 2 x 3 4 1 80 60 40 20 0 -20 0 2 3 x Solution using ode15s 1 4 Solution using ode45

80 dy/dx=x. -7.5 60 40 20 0 -20

Exercises:

Write Matlab programs to solve the following differential equations using 4th order Runge-Kutta method for the range and increment x shown. Also compare the accuracy of the solution by comparing the computed value with the values computed from the corresponding true solutions.

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

13 of 16

Applications

System of ordinary diff. Equations

An example of a non-stiff system is the system of equations describing the motion of a rigid body without external forces:

To simulate this system, create a function M-file rigid.m (you can use any file name) containing the equations:

Then use ode45 to calculate y(t) for t=0 to 12. The initial conditions are given in the third vector.

function dy = rigid(t,y) dy = zeros(3,1); % a column vector dy(1) = y(2) * y(3); dy(2) = y(1) * y(3); dy(3) = 0.51 * y(1) * y(2);

[t,y] = ode45('rigid',[0 12],[0 1 1]);

Now type [t,y] at command prompt and press enter. Matlab shall respond as under. The first column correspond to t and 2nd 3rd and 4th column represent the corresponding values of y1 , y2 , y3.
ans = 0 0 1.0000 1.0000 0.0001 0.0001 1.0000 1.0000 0.0001 0.0001 1.0000 1.0000 0.0002 0.0002 1.0000 1.0000 . . . .. . . . Type plot(t,y) at MATLAB command prompt. MATLAB shall respond as under.
1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 2 4 6 8 10 12

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

14 of 16

Applications

Finding Local Minima of a Function

Command(s) : fminsearch
a) find a local minimum of a function

For functions of a single variable: find the minimum of the function

f ( x ) = 1 xe2 x ,
first analytically and then numerically using the fminsearch command. Hint : Unlike most MATLAB operations, you dont need to define a vector x before using fminsearch. Simply type in the function as shown below. fminsearch looks for the local minimum which is closest to the initial guess you specify as the third parameter.

fminsearch ( ' f ( x ) ', x0 )

Analytically:

f ( x ) = e2 x ( 2 x 1) = 0
Numerically:

x=

1 2

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

15 of 16

Applications
b) For functions of several variables:
0 < x < 0 < y <

find the minimum of the function f ( x, y ) = cos ( xy ) sin ( x ) over the domain using the fminsearch command.

Follow these steps: plot the function over the given domain to approximately locate the minimum. use these approximate values of the coordinates as your initial guess Hint: The syntax for functions of several variables is slightly different. Define X = [ x y ] . Then in MATLAB, x and y will be represented by the two components: X (1) = x and X ( 2 ) = y , and the command becomes:

fminsearch ' f ( X (1) , X ( 2 ) ) ', X 0 (1) , X 0 (1)

% Plot the surface to visualize the minimum [ x , y ] = meshgrid ( 0:0.1:pi , 0:0.1:pi); f = cos(x.*y).*sin(x); surface(x,y,f);

% Transform x and y in X(1) and X(2) fminsearch( ' cos( x(1) * x(2)) * sin(x(1)) ',[1.5 2])

% The minimum is around the point [ 1.5 2 ]

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

16 of 16

Symbolic Computations

Section-VIII

Symbolic Computations

Prepared by: Engr. Tarique Ahmad Reviewed by Prof. Dr. Asad-ur-Rehman

1 of 10

MATLAB has the capability to manipulate expressions symbolically. There are tools to perform algebraic operations, differentiate and integrate functions, solve systems of equations, and solve ordinary differential equations. MATLAB adapted these tools from the software program Maple developed at the University of Waterloo, Canada.

Creating symbolic expressions

Variables can be declared symbolic variables with the command for example we can declare variables x; y; z; a; b; c, etc. as symbolic variable using MATLAB syms command

syms x y z a b c
This command is a short cut for the more elaborate command sym('x', 'y','z', 'a', 'b', ' c' ) , or even more deliberately, x = sym('x'),y = sym('y'),etc .We can then define expressions using these variables and these expressions can be manipulated symbolically. For example a matrix A can be defined by

> A = [ a b 1; 0 1 c; x 0 0 ]

Since A is a symbolic expression, we can calculate its determinant in terms of the variables a; b; c; x with the usual MATLAB command.

>> d = det(A)

S 3ymbolic Functions
A function f(x) can be defined in terms of a symbolic expression.Suppose we want to express f(x)=ax2+bx+c+2cos(x) as a symbolic expression.

>> syms x y z a b c >> f = a*x^2 + b*x +c + 2*cos(x)

Notice that we havent used array operators in defining the function. Actually, the array operators .^, .*, ./ cannot be used in symbolic expressions because

symbolic expressions are not applied directly to vectors and matrices.

Evaluating symbolic expressions

Next, how do we specify the values of the parameters in the expression, and how do we evaluate the symbolically defined function at a point? This is done using the substitution command subs. The syntax is subs(F, Old, New) Where, Old is a symbolic variable, a string representing variable name, or a string (quoted) expression.

New is a symbolic or numeric variable or expression. F is the symbolic expression to be evaluated

For example, if we wish to evaluate the function f at x = 2, we cannot execute f(2)to get f at x=2, as we did in case of using MATLAB calculation engine. To evaluate the symbolic function f defined above at x = 2, (leaving in the parameters a; b; c,) we will use MATLAB subs command from the Symbolic Toolbox:

>> subs(f,x,2)

The result is still a symbolic expression. Now specify the values of the parameters a, b, c as a=2; b=3; c=9, and use subs command to evaluate f in terms of x only.

>> g = subs(f, [a b c], [2 -3 9])

We still have a symbolic expression depending on variable x. Now let us try to evaluate this function at a particular point, say x = 1.5. For this, we make another substitution subs(g,x,-1.5) with the answer 9.1415.

>> subs(f,x,-1.5)

The result is a real number now. Off course, many variations are possible.

Symbolic Summations
MATLAB symbolic engine is capable of doing symbolic summations, even those with infinitely many terms. The symsum command is used for symbolic summation. Consider the symbolic sum
k =m

f (k ) ,

where f(k) is a function of an to calculate this using MATLAB

integer k. The syntax Symbolic Toolbox is: >>syms k >>symsum(f,k,m,n).

Summation expression often involves factorial (!). Since (!) is also a numeric command, we need to tell MATLAB

to use the symbolic form by including the term involving in the factorial. For example the expression (1+2*k)!, is translated as, >>sym('(1+2*k)!') Computing Infinite sums are particularly interesting in MATLAB symbolic toolbox. The following solutions for the three summation expression illustrate the necessary concepts of symbolic summation using MATLAB Symbolic Toolbox.

1. Evaluate

k =0

x 2 k +1 (1) (2k + 2)!

k

>> clear, syms k x; >>symsum((-1)^k*x^ (2*k+1)/ sym('(2*k+1)!'),k,0,Inf)

2. Evaluate

k =1

k2

>> clear, syms k; >> symsum((1/k^2),k,1,Inf)

xk 3. Evaluate k = 0 k!

>>clear, syms k x; >>symsum(x^k/sym('(k!)'),k,0,I nf)

GRAPHING EQUATIONS AND FUNCTIONS

MATLAB symbolic engine has commands for plotting symbolic equations and functions. Most of these graphing commands start with ez pronounced as easy. These commands include ezplot, ezsurf and ezmesh for 2D and 3D plotting. 2D D Plots The command ezplot is used primarily to graph functions which are defined symbolically. If f is defined by a symbolic expression and we wish to graph it on the interval [1; 5], we can do it with the one command, ezplot(f, [1,5]). For example

>> syms x >> f = cos(x)^2*exp(x) >> ezplot(f, [1,5]) >>grid on

Prepared by Tarique Ahmad

As another example consider the equation for a circle of radius 2, i.e. x2+y2 = 4. Step-1: Move everything to the left-hand side of the equation, i.e. x2 + y2 4 = 0 Step-2: Clear the workspace and define x and y as symbolic. clear; syms x y; Step-3: Define a variable equal to the left-hand ide of the equation, inside single quotation marks. LHS = 'x^2+y^2 - 4'; Step-4: Use ezplot. ezplot(LHS,[-3,3,-2.5,2.5])

>> clear; >>syms x y; >>LHS = 'x^2+y^2 - 4'; >>ezplot(LHS) >>axis([-3,3,-2.5,2.5]) >> grid on

The [-3,3,-2,2] tells MATLAB that the x scale is to go from -3 to +3, and the y scale is to go from -2.5 to +2.5. Note that we have assigned x2 + y2 - 4 to the character variable LHS (see the Workspace window or type >> whos). You can edit the figure in the Figure window. First save it as circle.fig using File/Save As so you can retrieve it later from the Command window, without having to recreate it. Here are some things you can do in the Figure window:

Click on the arrow icon and then double-click on any

portion of the figure youd like to edit. In this way

you can change the color and thickness of the curve, add points, change the axes, edit the title. Edit / Copy Figure to place it on the Windows Clipboard. Then you can paste it into a Word document, for example. While in the Figure window, press the Alt and Print Screen keys at the same time, to place the entire figure window on the clipboard. This can also be pasted into a Word document. 3D D Plots Three-dimensional plots can be created using the ezmesh and ezsurf commands. For these commands an explicit form must be used, i.e. z = f(x, y). It is the f(x, y) that is entered into the command. The code below illustrates some possibilities.

>>clear >>syms V x y >>V=sin(x)*sin(y); >>figure(2); >>ezsurf(V); >>colormap summer

In this example the function V ( x, y ) = sin( x) cos( x) is plotted. The 1st line clear variables and functions from memory. The 2nd line defines the variables V, x, y as symbolic variables. The 3rd line creates a blank figure window titled Figure 2

The 4th line is a symbolic implementation of MATLAB surf command. The MATLAB Surf command plots the colored parametric surface of the function V. Surf(V) uses x = 1:n and y = 1:m. The height, V, defined over a geometrically rectangular grid of x and y. The 5th line selects and implements a colormap scheme for the plotted figure. In this case the selected colormap scheme is summer. MATLAB provides many other colormap schemes like jet, bone, colorcube, cool, copper, flag, grey, hot, hsv, lines, pink, prism, spring, white and winter.

>>clear >>syms V x y >>V=sin(x)*sin(y) >>figure(3); >>ezmesh(V); >>colormap summer

In this example the function V ( x, y ) = sin( x) cos( x) is plotted once again. However this time we used ezmesh instead of ezsurf command. Ezmesh is a symbolic implementation of MATLAB mesh command. It gives 3D-mesh surfaces.

Note the difference ezmesh commands.

in

results

produced

by

ezsurf

and