Matlab Programming: Gerald W. Recktenwald Department of Mechanical Engineering Portland State University Gerry@me - Pdx.edu
Matlab Programming: Gerald W. Recktenwald Department of Mechanical Engineering Portland State University Gerry@me - Pdx.edu
Matlab Programming: Gerald W. Recktenwald Department of Mechanical Engineering Portland State University Gerry@me - Pdx.edu
These slides are a supplement to the book Numerical Methods with Matlab: Implementations and Applications, by Gerald W. Recktenwald, c 2001, Prentice-Hall, Upper Saddle River, NJ. These slides are c 2001 Gerald W. Recktenwald. The PDF version of these slides may be downloaded or stored or printed only for noncommercial, educational use. The repackaging or sale of these slides in any form, without written consent of the author, is prohibited. The latest version of this PDF le, along with other supplemental material for the book, can be found at www.prenhall.com/recktenwald.
Version 0.97
Overview
Script m-les Creating Side eects Function m-les Syntax of I/O parameters Text output Primary and secondary functions Flow control Relational operators Conditional execution of blocks Loops Vectorization Using vector operations instead of loops Preallocation of vectors and matrices Logical and array indexing Programming tricks Variable number of I/O parameters Indirect function evaluation Inline function objects Global variables
page 1
Preliminaries
Programs are contained in m-les Plain text les not binary les produced by word processors File must have .m extension m-le must be in the path Matlab maintains its own internal path The path is the list of directories that Matlab will search when looking for an m-le to execute. A program can exist, and be free of errors, but it will not run if Matlab cannot nd it. Manually modify the path with the path, addpath, and rmpath built-in functions, or with addpwd NMM toolbox function . . . or use interactive Path Browser
page 2
Script Files
Not really programs No input/output parameters Script variables are part of workspace Useful for tasks that never change Useful as a tool for documenting homework: Write a function that solves the problem for arbitrary parameters Use a script to run function for specic parameters required by the assignment
Free Advice: Scripts oer no advantage over functions. Functions have many advantages over scripts. Always use functions instead of scripts.
page 3
(1)
Enter statements in le called tanplot.m 1. Choose New. . . from File menu 2. Enter lines listed below Contents of tanplot.m:
theta = linspace(1.6,4.6); tandata = tan(theta); plot(theta,tandata); xlabel(\theta (radians)); ylabel(tan(\theta)); grid on; axis([min(theta) max(theta) -5 5]);
page 4
(2)
1 tan()
-1
-2
-3
-4
-5
2.5
3 (radians)
3.5
4.5
If the plot needs to be changed, edit the tanplot script and rerun it. This saves the eort of typing in the commands. The tanplot script also provides written documentation of how to create the plot.
Example: Put a % character at beginning of the line containing the axis command, then rerun the script
page 5
All variables created in a script le are added to the workplace. This may have undesirable eects because
Variables already existing in the workspace may be overwritten The execution of the script can be aected by the state variables in the workspace.
% %
plot(x,y) % Generate the plot and label it xlabel(x axis, unknown units) ylabel(y axis, unknown units) title(Plot of generic x-y data set)
page 6
The D, x, and y variables are left in the workspace. These generic variable names might be used in another sequence of calculations in the same Matlab session. See Exercise 10 in Chapter 4.
page 7
Occur when a module changes variables other than its input and output parameters Can cause bugs that are hard to track down Cannot always be avoided
Create and change variables in the workspace Give no warning that workspace variables have changed
Because scripts have side eects, it is better to encapsulate any mildly complicated numerical in a function m-le
page 8
Functions are subprograms: Functions use input and output parameters to communicate with other functions and the command window Functions use local variables that exist only while the function is executing. Local variables are distinct from variables of the same name in the workspace or in other functions. Input parameters allow the same calculation procedure (same algorithm) to be applied to dierent data. Thus, function m-les are reusable. Functions can call other functions. Specic tasks can be encapsulated into functions. This modular approach enables development of structured solutions to complex problems.
page 9
outArgs is a comma-separated list of variable names [ ] is optional if there is only one parameter functions with no outArgs are legal
inArgs is a comma-separated list of variable names functions with no inArgs are legal
page 10
twosum.m two inputs, no output threesum.m three inputs, one output addmult.m two inputs, two outputs
page 11
twosum.m
function twosum(x,y) % twosum Add two matrices % and print the result x+y
threesum.m
function s = threesum(x,y,z) % threesum Add three variables % and return the result s = x+y+z;
addmult.m
function [s,p] = addmult(x,y) % addmult Compute sum and product % of two matrices s = x+y; p = x*y;
page 12
Notes:
1. The result of the addition inside twosum is exposed because the x+y expression does not end in a semicolon. (What if it did?) 2. The strange results produced by twosum(one,two) are obtained by adding the numbers associated with the ASCII character codes for each of the letters in one and two. Try double(one) and double(one) + double(two).
page 13
In this example, the x and y variables dened in the workspace are distinct from the x and y variables dened in twosum. The x and y in twosum are local to twosum.
page 14
Note: The last statement produces no output because the assignment expression ends with a semicolon. The value of 24 is stored in b.
page 15
Note: addmult requires two return variables. Calling addmult with no return variables or with one return variable causes undesired behavior.
page 16
Values are communicated through input arguments and output arguments. Variables dened inside a function are local to that function. Local variables are invisible to other functions and to the command environment. The number of return variables should match the number of output variables provided by the function. This can be relaxed by testing for the number of return variables with nargout (See 3.6.1.).
page 17
It is usually desirable to print results to the screen or to a le. On rare occasions it may be helpful to prompt the user for information not already provided by the input parameters to a function.
Inputs to functions:
input function can be used (and abused!). Input parameters to functions are preferred.
Text output from functions:
disp function for simple output fprintf function for formatted output.
page 18
The input function can be used to prompt the user for numeric or string input.
>> x = input(Enter a value for x); >> yourName = input(Enter your name,s);
Prompting for input betrays the Matlab novice. It is a nuisance to competent users, and makes automation of computing tasks impossible.
Free Advice: Avoid using the input function. Rarely is it necessary. All inputs to a function should be provided via the input parameter list. Refer to the demonstration of the inputAbuse function in 3.3.1.
page 19
Output to the command window is achieved with either the disp function or the fprintf function. Output to a le requires the fprintf function. disp fprintf Simple to use. Provides limited control over appearance of output. Slightly more complicated than disp. Provides total control over appearance of output.
page 20
Syntax:
disp(outMatrix)
where outMatrix is either a string matrix or a numeric matrix. Examples: Numeric output
>> disp(5) 5 >> x = 1:3; 1 >> y = 3-x; 1 2 disp(x) 2 3 disp([x; y]) 2 3 1 0
>> disp([x y]) ??? All matrices on a row in the bracketed expression must have the same number of rows.
Note: The last statement shows that the input to disp must be a legal matrix.
page 21
>> t = Earlier versions used LINPACK and EISPACK; >> disp([s; t]) ??? All rows in the bracketed expression must have the same number of columns. >> disp(char(s,t)) MATLAB 6 is built with LAPACK Earlier versions used LINPACK and EISPACK
The disp[s; t] expression causes an error because s has fewer elements than t. The built-in char function constructs a string matrix by putting each input on a separate row and padding the rows with blanks as necessary.
>> S = char(s,t); >> length(s), length(t), ans = 29 ans = 41 ans = 41 length(S(1,:))
page 22
The num2str function is often used to with the disp function to create a labeled output of a numeric value. Syntax:
stringValue = num2str(numericValue)
0 0 1
>> S = num2str(A) S = 1 0 0 0 1 0 0 0 1
page 23
Although A and S appear to contain the same values, they are not equivalent. A is a numeric matrix, and S is a string matrix.
>> clear >> A = eye(3); S = num2str(A); >> A-S ??? Error using ==> Matrix dimensions must agree. >> A-B ans = 0 0 0 >> whos Name A B S ans B = str2num(S);
0 0 0
0 0 0
Bytes 72 72 42 72
page 24
(1)
page 25
(2)
The
disp([x = ,num2str(x)]);
construct works when x is a row vector, but not when x is a column vector or matrix
>> z = y; >> disp([z = ,num2str(z)]) ??? All matrices on a row in the bracketed expression must have the same number of rows.
page 26
(3)
The same eect is obtained by simply entering the name of the variable with no semicolon at the end of the line.
>> z z = 1 2 3 4 (enter z and press return)
page 27
Alternatively, a second parameter can be used to control the precision of the output of num2str
>> disp([pi = ,num2str(pi,2)]) pi = 3.1 >> disp([pi = ,num2str(pi,4)]) pi = 3.142 >> disp([pi = ,num2str(pi,8)]) pi = 3.1415927
page 28
Syntax:
fprintf(outFormat,outVariables) fprintf(fileHandle,outFormat,outVariables)
uses the outFormat string to convert outVariables to strings that are printed. In the rst form (no fileHandle) the output is displayed in the command window. In the second form, the output is written to a le referred to by the fileHandle (more on this later).
Notes to C programmers: 1. The Matlab fprintf function uses single quotes to dene the format string. 2. The fprintf function is vectorized. (See examples below.)
Example:
>> x = 3; >> fprintf(Square root of %g is %8.6f\n,x,sqrt(x)); The square root of 3 is 1.732051
page 29
The outFormat string species how the outVariables are converted and displayed. The outFormat string can contain any text characters. It also must contain a conversion code for each of the outVariables. The following table shows the basic conversion codes.
Code %s %d %f %e %g \n \t Conversion instruction format as a string format with no fractional part (integer format) format as a oating-point value format as a oating-point value in scientic notation format in the most compact form of either %f or %e insert newline in output string insert tab in output string
page 30
In addition to specifying the type of conversion (e.g. %d, %f, %e) one can also specify the width and precision of the result of the conversion. Syntax:
%wd %w.pf %w.pe
where w is the number of characters in the width of the nal result, and p is the number of digits to the right of the decimal point to be displayed. Examples:
Format String %14.5f Meaning use oating point format to convert a numerical value to a string 14 characters wide with 5 digits after the decimal point use scientic notation format to convert numerical value to a string 12 characters wide with 3 digits after the decimal point. The 12 characters for the string include the e+00 or e-00 (or e+000 or e-000 on WindowsTM )
%12.3e
page 31
page 32
The fprintf function is vectorized. This enables printing of vectors and matrices with compact expressions. It can also lead to some undesired results. Examples:
>> x = 1:4; y = sqrt(x); >> fprintf(%9.4f\n,y) 1.0000 1.4142 1.7321 2.0000
The %9.4f format string is reused for each element of y. The recycling of a format string may not always give the intended result.
>> x = 1:4; y = sqrt(x); >> fprintf(y = %9.4f\n,y) y = 1.0000 y = 1.4142 y = 1.7321 y = 2.0000
page 33
Vectorized fprintf cycles through the outVariables by columns. This can also lead to unintended results
>> A = [1 2 3; 4 5 6; 7 8 9] A = 1 2 3 4 5 6 7 8 9 >> fprintf(%8.2f %8.2f %8.2f\n,A) 1.00 4.00 7.00 2.00 5.00 8.00 3.00 6.00 9.00
page 34
Many times a tabular display of results is desired. The boxSizeTable function listed below, shows how the fprintf function creates column labels and formats numeric data into a tidy tabular display. The for loop construct is discussed later in these slides.
function boxSizeTable % boxSizeTable Demonstrate tabular output with fprintf % --- labels and sizes for shiping containers label = char(small,medium,large,jumbo); width = [5; 5; 10; 15]; height = [5; 8; 15; 25]; depth = [15; 15; 20; 35]; vol = width.*height.*depth/10000; % volume in cubic meters fprintf(\nSizes of boxes used by ACME Delivery Service\n\n); fprintf(size width height depth volume\n); fprintf( (cm) (cm) (cm) (m^3)\n); for i=1:length(width) fprintf(%-8s %8d %8d %8d %9.5f\n,... label(i,:),width(i),height(i),depth(i),vol(i)) end
Note: length is a built-in function that returns the number of elements in a vector. width, height, and depth are local variables in the boxSizeTable function.
NMM: Matlab Programming page 35
page 36
File Output with fprintf requires creating a le handle with the fopen function. All aspects of formatting and vectorization discussed for screen output still apply. Example: Writing contents of a vector to a le.
x = ... % content of x fout = fopen(myfile.dat,wt); % open myfile.dat fprintf(fout, k x(k)\n); for k=1:length(x) fprintf(fout,%4d %5.2f\n,k,x(k)); end fclose(fout) % close myfile.dat
page 37
To enable the implementation of computer algorithms, a computer language needs control structures for
page 38
Relational operators are used in comparing two values. Operator < <= > >= ~= Meaning less than less than or equal to greater than greater than or equal to not equal to
The result of applying a relational operator is a logical value, i.e. the result is either true or false. In Matlab any nonzero value, including a non-empty string, is equivalent to true. Only zero is equivalent to false.
Note: The <=, >=, and ~= operators have = as the second character. =<, => and =~ are not valid operators.
page 39
Relational operations can also be performed on matrices of the same shape, e.g.,
>> x = 1:5; >> z = x>y z = 0 0 y = 5:-1:1;
page 40
Logical Operators
Logical operators are used to combine logical expressions (with and or or), or to change a logical value with not Operator & | ~ Meaning and or not
Examples:
>> a = 2; b = 4; >> aIsSmaller = a < b; >> bIsSmaller = b < a; >> bothTrue = aIsSmaller & bIsSmaller bothTrue = 0 >> eitherTrue = aIsSmaller | bIsSmaller eitherTrue = 1 >> ~eitherTrue ans = 0
page 41
Summary:
Relational operators involve comparison of two values. The result of a relational operation is a logical (True/False) value. Logical operators combine (or negate) logical values to produce another logical value. There is always more than one way to express the same comparison
Free Advice:
To get started, focus on simple comparison. Do not be afraid to spread the logic over multiple lines (multiple comparisons) if necessary. Try reading the test out loud.
page 42
Conditional Execution
Conditional Execution or Branching: As the result of a comparison, or another logical (true/false) test, selected blocks of program code are executed or skipped. Conditional execution is implemented with if, if...else, and if...elseif constructs, or with a switch construct. There are three types of if constructs 1. Plain if 2. if...else 3. if...elseif
page 43
if Constructs
Syntax:
if expression block of statements end
Example:
if a < 0 disp(a is negative); end
page 44
if. . . else
Multiple choices are allowed with if. . . else and if. . . elseif constructs
if x < 0 error(x is negative; sqrt(x) is imaginary); else r = sqrt(x); end
page 45
if. . . elseif
Its a good idea to include a default else to catch cases that dont match preceding if and elseif blocks
if x > 0 disp(x elseif x < disp(x else disp(x end
page 46
A switch construct is useful when a test value can take on discrete values that are either integers or strings. Syntax:
switch expression case value1, block of statements case value2, block of statements . . . otherwise, block of statements end
Example:
color = ...; % switch color case red disp(Color case blue disp(Color case green disp(Color otherwise disp(Color end
NMM: Matlab Programming
color is a string
page 47
Repetition or Looping A sequence of calculations is repeated until either 1. All elements in a vector or matrix have been processed or 2. The calculations have produced a result that meets a predetermined termination criterion Looping is achieved with for loops and while loops.
page 48
for loops
for loops are most often used when each element in a vector or matrix is to be processed. Syntax:
for index = expression block of statements end
page 49
Note: In the last example, x is a scalar inside the loop. Each time through the loop, x is set equal to one of the columns of 0:pi/15:pi.
page 50
while loops are most often used when an iteration is repeated until some termination criterion is met. Syntax:
while expression block of statements end
rk =
1 2
rk1 +
x rk1
r = ... % initialize rold = ... while abs(rold-r) > delta rold = r; r = 0.5*(rold + x/rold); end
NMM: Matlab Programming
page 51
It is (almost) always a good idea to put a limit on the number of iterations to be performed by a while loop. An improvement on the preceding loop,
maxit = 25; it = 0; while abs(rold-r) > delta & it<maxit rold = r; r = 0.5*(rold + x/rold); it = it + 1; end
page 52
The break and return statements provide an alternative way to exit from a loop construct. break and return may be applied to for loops or while loops. break is used to escape from an enclosing while or for loop. Execution continues at the end of the enclosing loop construct. return is used to force an exit from a function. This can have the eect of escaping from a loop. Any statements following the loop that are in the function body are skipped.
page 53
page 54
page 55
break is used to escape the current while or for loop. return is used to escape the current function.
function k = demoBreak(n) ... while k<=n if x(k)>0.8 break; end k = k + 1; end function k = demoReturn(n) ... while k<=n if x(k)>0.8 return; end k = k + 1; end
page 56
Vectorization
Vectorization is the use of vector operations (Matlab expressions) to process all elements of a vector or matrix. Properly vectorized expressions are equivalent to looping over the elements of the vectors or matrices being operated upon. A vectorized expression is more compact and results in code that executes faster than a non-vectorized expression. To write vectorized code:
Use vector operations instead of loops, where applicable Pre-allocate memory for vectors and matrices Use vectorized indexing and logical functions
Non-vectorized code is sometimes called scalar code because the operations are performed on scalar elements of a vector or matrix instead of the vector as a whole.
Free Advice: Code that is slow and correct is always better than code that is fast and incorrect. Start with scalar code, then vectorize as needed.
page 57
Scalar Code
for k=1:length(x) y(k) = sin(x(k)) end
Vectorized equivalent
y = sin(x)
page 58
Preallocate Memory
y = ... % some computation to define y s = zeros(size(y)); for j=1:length(y) if y(j)>0 s(j) = sqrt(y(j)); end end
page 59
Thorough vectorization of code requires use of array indexing and logical indexing. Array Indexing: Use a vector or matrix as the subscript of another matrix:
>> x = sqrt(0:4:20) x = 0 2.0000 2.8284 >> i = [1 2 5]; >> y = x(i) y = 0 2
3.4641
4.0000
4.47210
The x(i) expression selects the elements of x having the indices in i. The expression y = x(i) is equivalent to
k = 0; for i = [1 2 5] k = k + 1; y(k) = x(i); end
page 60
Logical Indexing: Use a vector or matrix as the mask to select elements from another matrix:
>> x = sqrt(0:4:20) x = 0 2.0000 2.8284 >> j = find(rem(x,2)==0) j = 1 2 5 >> z = x(j) z = 0 2
3.4641
4.0000
4.47210
The j vector contains the indices in x that correspond to elements in x that are integers.
page 61
Example: Vectorization of Scalar Code We just showed how to pre-allocate memory in the code snippet:
y = ... % some computation to define y s = zeros(size(y)); for j=1:length(y) if y(j)>0 s(j) = sqrt(y(j)); end end
In fact, the loop can be replaced entirely by using logical and array indexing
y = ... s = zeros(size(y)); i = find(y>0); s(y>0) = sqrt(y(y>0)) % % some computation to define y indices such that y(i)>0
If we dont mind redundant computation, the preceding expressions can be further contracted:
y = ... % some computation to define y s = zeros(size(y)); s(y>0) = sqrt(y(y>0))
NMM: Matlab Programming
page 62
Vectorized Code
B(:,1) = A(:,1);
page 63
Vectorized Code
B(1,2:3) = A(2:3,3)
page 64
Deus ex Machina
Variable number of I/O parameters Indirect function evaluation with feval In-line function objects (Matlab version 5.x) Global Variables
page 65
Each function has internal variables, nargin and nargout. Use the value of nargin at the beginning of a function to nd out how many input arguments were supplied. Use the value of nargout at the end of a function to nd out how many input arguments are expected. Usefulness:
Allows a single function to perform multiple related tasks. Allows functions to assume default values for some inputs, thereby simplifying the use of the function for some tasks.
page 66
Consider the built-in plot function Inside the plot function nargin plot(x,y) plot(x,y,s) plot(x,y,s--) plot(x1,y1,s,x2,y2,o) h = plot(x,y) 2 3 3 6 2 nargout 0 0 0 0 1
The values of nargin and nargout are determined when the plot function is invoked. Refer to the demoArgs function in Example 3.13
page 67
Allows routines to be written to process an arbitrary f (x). Separates the reusable algorithm from the problem-specic code.
feval is used extensively for root-nding (Chapter 6), curve-tting (Chapter 9), numerical quadrature (Chapter 11) and numerical solution of initial value problems (Chapter 12).
page 68
page 69
Use of feval
function s = fsum(fun,a,b,n) % FSUM Computes the sum of function values, f(x), at n equally % distributed points in an interval a <= x <= b % % Synopsis: s = fsum(fun,a,b,n) % % Input: fun = (string) name of the function to be evaluated % a,b = endpoints of the interval % n = number of points in the interval x = linspace(a,b,n); y = feval(fun,x); s = sum(y); % % % create points in the interval evaluate function at sample points compute the sum
function y = sincos(x) % SINCOS Evaluates sin(x)*cos(x) for any input x % % Synopsis: y = sincos(x) % % Input: x = angle in radians, or vector of angles in radians % % Output: y = value of product sin(x)*cos(x) for each element in x y = sin(x).*cos(x);
page 70
Use
myFun = inline( x.^2 - log(x) );
Usefulness:
Eliminates need to write separate m-les for functions that evaluate a simple formula. Useful in all situations where feval is used.
NMM: Matlab Programming page 71
Global Variables
localFun.m
(a,b,c) function d = localFun(a,b,c) ... d d = a + b^c
Communication of values via input and output variables and global variables shared by the workspace and function workspace
x = 1; y = 2; global ALPHA; ... ALPHA = 1.2; z = globalFun(x,y) (x,y) (a,b)
globalFun.m
function d = globalFun(a,b) ... global ALPHA ... d = a + b^ALPHA;
Usefulness:
Allows bypassing of input parameters if no other mechanism (such as pass-through parameters) is available. Provides a mechanism for maintaining program state (GUI applications)
page 72