Module 1: Programming in MATLAB
1/148
Relational and Logical Operations
2/148
FPRINTF: Alternate Displaying Function
Combine literal text with numeric data. ‚ Number of digits to display
fprintf(’There are %d days in a year.\n’, 365) ‚ Complex number
z = exp(1i*pi/4);
fprintf(’%f+%fi\n’, real(z), imag(z));
3/148
FPRINTF: Formatting Operator
%[field width][precision][conversion character]
e.g. %12.5f.
‚ %: marks the beginning of a formatting operator
‚ [field width]: maximum number of characters to print; optional
‚ [precision] number of digits to the right of the decimal point; optional ‚ [conversion character]
%d integer
%f xed-point notation
%e exponential notation
%g the more compact of %f or %e %s string array
%x hexadecimal
4/148
Relational Operators
How are two numbers X and Y related? ‚ [X>Y] Is X greater than Y?
‚ [X
end
14/148
Main Fragment – another way
if xc <= L || xc >= R
disp(’Either xc <= L or xc >= R’)
else
fprintf(’Interior critical point at x_c = %5.2f\n’, xc)
end
15/148
Main Fragment – yet another way
if ~(xc <= L || xc >= R)
fprintf(’Interior critical point at x_c = %5.2f\n’, xc)
else
disp(’Either xc <= L or xc >= R’)
end
16/148
The simplest if statement?
So far, we have seen
‚ if-else statement
‚ if-elseif-else statement
The simplest if statement is of the form
if [condition]
[statements to run]
end
17/148
Input Errors
If a user mistakenly provides L that is larger than R, x it silently by swapping L and R.
if L > R
tmp = L;
L = R;
R = tmp;
end
I will show you how to send an error message and halt a program later.
18/148
Exercise 1: Simple Minimization Problem
Question
Write a program which xmin P rL, Rs at which qpxq is minimized and the minimum value qpxminq.
‚ This can be done with if-elseif-else
19/148
Exercise 2: Leap Year
Question
Write a script which determines whether a given year is a leap year or not. A year is a leap year if
‚ it is a multiple of 4;
‚ it is not a multiple of 100; ‚ it is a multiple of 400.
Useful: mod function.
20/148
Pseudocode
if [YEAR] is not divisible by 4
it is a common year
elseif [YEAR] is not divisible by 100
it is a leap year
elseif [YEAR] is not divisible by 400
it is a common year
else
it is a leap year
end
21/148
Exercise 3: Angle Finder
Question
Let x and y be given, not both zero. Determine the angle ◊ P p ́fi, fis between the positive x-axis and the line segment connecting the origin to px, yq.
Four quandrants:
‚ 1st or 4th (x °“ 0): ◊ “ tan ́1py{xq
‚ 2nd(x†0,y°“0):◊“tan ́1py{xq`fi ‚ 3rd(x†0,y†0):◊“tan ́1py{xq ́fi
Useful: atan (inverse tangent function)
22/148
Extended Inverse Tangent
if x > 0
theta = atan(y/x)
elseif y >= 0
theta = atan(y/x) + pi
else
theta = atan(y/x) – pi
end
‚ MATLAB provides a function that exactly does this: atan2(x,y).
‚ Further Exploration: What would you do if you are asked to nd the angle ◊ P r0, 2fiq, with atan alone or with atan2?
23/148
FOR-Loops
24/148
Approximating fi
Supposethecirclex2 `y2 “n2,nPN,isdrawnon
graph paper.
‚ The area of the circle can be approximated by
counting the number uncut grids, Nin. fin2 « Nin,
and so
‚ Using symmetry, may only count the grids in the rst quadrant and modify the formula accordingly:
fi « 4Nin,1 , n2
where Nin,1 is the number of inscribed grids in the rst quadrant.
fi « Nin . n2
25/148
Approximating fi
Problem Statement
Write a script that inputs an integer n and displays the approximation of fi by fl “4Nin,1,
n n2 along with the (absolute) error |fln ́ fi|.
Note. The approximation gets enhanced and approaches the true value of fi as n Ñ 8.
26/148
Strategy: Iterate
The key to this problem is to count the number of uncut grids in the rst quadrant programmatically.
row 10
row 9 row 8 row 7 row 6 row 5 row 4 row 3 row 2 row 1
Set Nin,1 “ 0.
Count the number of uncut grids in row 1. Add that to Nin,1.
Count the number of uncut grids in row 2. Add that to Nin,1.
.
Set fl10 “ 4Nin,1{102.
Count the number of uncut grids in row 10. Add that to Nin,1.
27/148
MATLAB Way
The repeated counting can be delegated to MATLAB using for-loop. The procedure outlined above turns into
Assume n is initialized and set Nin,1 to zero. for k = 1:n
Count the number of uncut grids in row k. Add that to Nin,1.
end
Set fl10 “ 4Nin,1{102.
28/148
Counting Uncut Tiles
The problem is reduced to counting the number of uncutgridsineachrow.
row 10
row 9
row 8
row 7
row 6
row 5
row 4
row 3
row 2
row 1
‚
The x-coordinate of the intersection of the top edge of the kth row and the circle
x2 ` y2 “ n2 is a
x“ n2 ́k2.
px,yq For n “ 10 and k “ 5:
x“?n2 ́k2
“ ?102 ́ 52 “ 8.6602 . . .
‚ Thenumberofuncutgridsinthekthrowisthe largest inteager less than or equal to this value,
i.e.,
t n2 ́ k2u. ( oor function)
‚ MATLAB provides floor.
29/148
Main Fragment Using FOR-Loop
N1 = 0;
for k = 1:n
m = floor(sqrt(n^2 – k^2));
N1 = N1 + m;
end
rho_n = 4*N1/n^2;
Exercise. Complete the program.
30/148
Exercise 1: Overestimation
Question
Note that fln is always less than fi. If N1 denotes the total number of grids, both cut and uncut, within the quarter disk, then μn “ 4N1{n2 is always larger than fi. Modify the previous (complete) script so that it prints fln, μn, and
μn ́fln.
‚ ceil, an analogue of floor, is useful.
31/148
Notes on FOR-Loop
‚ The construct is used when a code fragment needs to be repeatedly run. The number of repetition is known in advance.
for
end
‚ Examples:
for k = 1:3
fprintf(’k = %d\n’, k)
end
nIter = 100;
for k = 1:nIter
end
fprintf(’k = %d\n’, k)
32/148
Caveats
Run the following script and observe the displayed result.
for k = 1:3
disp(k)
k = 17;
disp(k)
end
‚ The loop header k = 1:3 guarantees that k takes on the values 1, 2, and 3, one at a time even if k is modi ed within the loop body.
‚ However, it is a recommended practice that the value of the loop variable is never modi ed in the loop body.
33/148
Simulation Using rand
rand is a built-in function which generate a (uniform) “random” number between 0 and 1. Try:
for k = 1:10
x = rand();
fprintf(’%10.6f\n’, x);
end
Let’s use this function to solve:
Question
A stick with length 1 is split into two parts at a random breakpoint. On average, how long is the shorter piece?
34/148
Program Development – Single Instance
Consider breaking one stick.
‚ Random breakage can be simulated with rand; denote by x P p0, 1q.
‚ The length of the shorter piece can be determined using if-construct; denote by s P p0, 1{2q.
x = rand();
if x <= 0.5
s = x;
else
s = 1-x
% x: the location of breakage % if x§0.5
% shorter part has length x
% otherwise
% shorter part has length 1 ́x
end
35/148
Program Development – Multiple Instances
‚ Repeat the previous multiple times using a for-loop. Pseudocode: if 1000 breaks are to be simulated:
nBreaks = 1000;
for k = 1:nBreaks end
‚ But how are calculating the average length of the shorter pieces?
36/148
Calculating Average Using Loop
Recall how the total number of uncut Similarly, we can compute an average grids were calculated using iterations. by:
Assume n is initialized and set Nin,1 to zero.
for k = 1:n
Assume n is initialized and set s to zero.
for k = 1:n
Count the number of uncut grids in row k. Add that to Nin,1.
end
The value of N,1 is the total numbers of uncut grids.
end
Set savg “ s{n.
Simulate a break and nd the length of the shorter piece. Add that to s.
37/148
Complete Solution
nBreaks = 1000;
s = 0;
for k = 1:nBreaks
x = rand();
if x <= 0.5
s = s + x;
else
end end
s_avg = s/nBreaks;
s = s + (1-x);
38/148
Exercise 2: Game of 3-Stick
Game: 3-Stick
Pick three sticks each having a random length between 0 and 1. You win if you can form a triangle using the sticks; otherwise, you lose.
Question
Estimate the probability of winning a game of 3-Stick by simulating one million games and counting the number of wins.
Win
Lose
39/148
WHILE-Loops
40/148
Understanding Loops
Question 1
How many lines of output are produced by the following script?
for k = 100:200
disp(k)
end
A 99 B 100 C 101 D 200
41/148
Understanding Loops
Question 2
How many lines of output are produced by the following script?
for k = 100:200
if mod(k,2) == 0
disp(k) end
end
A 50
B 51
C 100
D 101
42/148
FOR-Loop: Tips
‚ Basic loop header:
for
‚ To adjust starting value:
for
‚ To adjust step size:
for
43/148
Examples
‚ Toiterateover1,3,5,...,9: for k = 1:2:9
or
for k = 1:2:10
‚ To iterate over 10, 9, 8, ..., 1: for k = 10:-1:1
[stepsize=2]
[ negative step size ]
44/148
Need for Another Loop
‚ For-loops are useful when the number of repetitions is known in advance.
"Simulate the tossing of a fair coin 100 times and print the number of Heads." ‚ It is not very suitable in other situations such as
"Simulate the tossing of a fair coin until the gap between the number of Heads and that of Tails reaches 10."
We need another loop construct that terminates as soon as |NH ́ NT| “ 10.
45/148
WHILE-Loop Basics
WHILE-loop is used when a code fragment needs to be executed repeatedly while a certain condition is true.
while
end
‚ The number of repetitions is not known in advance.
‚ The continuation criterion is a boolean expression, which is evaluated at
the start of the loop.
‚ If it is true, the loop body is executed. Then the boolean expression is evaluated again.
‚ If it is false, the ow of control is passed to the end of the loop.
46/148
Simple WHILE-Loop Examples
k = 1; n = 10;
while k <= n
k = 1;
while 2^k < 5000
k = k+1;
end
fprintf(’k = %d\n’, k)
end
fprintf(’k = %d\n’, k)
k = k+1;
47/148
FOR-Loop to WHILE-Loop
A for-loop can be written as a while-loop. For example,
FOR
s = 0;
for k = 1:4
s = s + k;
fprintf(’%2d %2d\n’, k, s)
end
WHILE
k = 0; s = 0;
while k < 4
k = k + 1; s = s + k;
fprintf(’%2d %2d\n’, k, s)
end
‚ Note that k needed to be initialized before the while-loop.
‚ The variable k needed to be updated inside the while-loop body.
48/148
Up/Down Sequence
Question
Pick a random integer between 1 and 1,000,000. Call the number n and repeat the following process:
‚ If n is even, replace n by n{2.
‚ Ifnisodd,replacenby3n`1.
Does it ever take more than 1000 updates to reach 1?
‚ To generate a random integer between 1 and k, use randi, e.g.,
randi(k)
‚ To test whether a number n is even or odd, use mod, e.g., mod(n, 2) == 0
49/148
Attempt Using FOR-Loop
for step = 1:1000
if mod(n,2) == 0
n = n/2;
else
n = 3*n + 1;
end
fprintf(’ %4d %7d\n’, step, n)
end
‚ Note that once n becomes 1, the central process yields the following
pattern: 1,4,2,1,4,2,1,...
‚ This program continues to run even after n becomes 1.
50/148
Solution Using WHILE-Loop
step = 0;
while n > 1
if mod(n,2) == 0
n = n/2;
else
n = 3*n + 1;
end
step = step + 1;
fprintf(’ %4d %7d\n’, step, n)
end
‚ This shuts down when n becomes 1!
51/148
Exercise: Gap of 10
Question
Simulate the tossing of a fair coin until the gap between the number of Heads and that of Tails reaches 10.
52/148
Summary
‚ For-loop is a programming construct to execute statements repeatedly.
for
end
‚ While-loop is another construct to repeatedly execute statements. Repetition is controlled by the termination criterion.
while
end
53/148
Arrays in MATLAB
54/148
Introduction to Arrays
Vectors and matrices are often collectively called arrays. Notation
‚ Rm (or Cm): the set of all real (or complex) column vectors with m elements.
‚ Rmˆn (or Cmˆn): the set of all real (or complex) m ˆ n matrices.
‚ IfvPRm withv“pv1,v2,...,vmqT,thenfor1§i§m,vi PRiscalled
the ith element or the ith index of v.
‚ IfAPRmˆn withA“pai,jq,thenfor1§i§mand1§j §n,ai,j PRis the element in the ith row and jth column of A.
55/148
Creating Arrays
‚ A row vector is created by x = [1 3 5 7];
x = [1,3,5,7];
‚ A column vector is created by y = [6; 1; 4];
y = [6 1 4].’;
‚ A matrix is formed by
A = [3 1 2 3;
1 5 6 5;
4 9 5 8];
x
1
3
5
7
3
1
2
3
1
5
6
5
4
9
5
8
6
1
4
yA
The MATLAB expression x.’ means xT while x’ means xH “ px ̊qT.
56/148
Shape of Arrays
‚ To nd the number of elements of a vector: length(x)
length(y)
‚ To nd the number of rows/columns of an
array:
size(A,1) % # of rows
size(A,2) % # of cols
size(A) % both
‚ To nd the total number of elements of an array:
numel(A)
4
1
3
5
7
1
1
6
1
4
secondary
3
size(A,2)
4
3
1
2
3
1
5
6
5
4
9
5
8
3
57/148
size(A,1)
primary
Shape of Arrays (Notes)
‚ For a matrix A, length(A) yields the larger of the two dimensions. ‚ The result of size(A) can be stored in two di erent ways:
szA = size(A)
[m, n] = size(A)
Q. How are they di erent?
‚ All of the following generate empty arrays.
[]
[1:0]
[1:0].’
Q. What are their sizes? What are their numel values?
58/148
Getting/Setting Elements of Arrays
‚ To access the ith element of a vector: x(2)
y(3)
‚ To access the pi, jq-element of a matrix: A(2,4)
‚ To assign values to a speci c element: x(2) = 2
A(2,4) = 0
‚ Indices start at 1 in MATLAB, not at 0!
1st
2nd 3rd 4th
1
3
5
7
6
1
4
secondary
1st 2nd 3rd
1st
1st 2nd 3rd
2nd 3rd 4th
3
1
2
3
1
5
6
5
4
9
5
8
59/148
primary
Linear Indexation and Straightening of Matrix
‚ MATLAB uses column-major layout by default, meaning that the elements of the columns are contiguous in memory.
‚ Consequently, one can get/set an element of a matrix using a single index.
A(8)
‚ An array can be put into a column vector using
A(:)
A A(:)
3123 1565 4 9 5 8
3 1st
1 2nd
4 3rd
1 4th
5 5th
9 6th
2 7th
6 8th
5 9th 3 10th
5 11th 8 12th
60/148
Two Kinds of Transpose
‚ The transpose of an array: AT A.’
‚ The conjugate transpose of an array: AH “A ̊ “AT
A’
‚ IfAPRmˆn,AH “AT.So,ifAisareal array, A.’ and A’ are equivalent.
A
3
1
2
3
1
5
6
5
4
9
5
8
3
1
4
1
5
9
2
6
5
3
5
8
A.'
61/148
Standard Arithmetic Operation
Standard arithmetic operations seen in linear algebra are executed using the familiar symbols.
‚ LetA,BPRmˆn andcPR.
‚ A ̆B:elementwiseaddition/subtraction (A ̆B) ‚ A ̆c:shiftingallelementsofAby ̆c (A ̆c)
‚ LetAPRmˆp,BPRpˆn,andcPR.
‚ A*B:themˆnmatrixobtainedbythelinearalgebraicmultiplication (AB) ‚ c*A: scalar multiple of A (cA)
‚ LetAPRmˆm andnPN.
‚
A^n: the n-th power of A; the same as A*A* ̈ ̈ ̈ *A (n times) (An)
62/148
Standard Arithmetic Operation – Inner Products
Let x, y P Rm be column vectors. The inner product of x and y is calculated by T ÿm
x y“x1y1 `x2y2 ` ̈ ̈ ̈`xmym “ xjyj PR. j“1
In MATLAB, simply type x’*y.
*
1x4
4x1 1x1
63/148
Standard Arithmetic Operation – Outer Products
Let x P Rm, y P Rn be column vectors. The outer product of x and y is calculated by
»— x 1 y 1 x 1 y 2 ̈ ̈ ̈ x 1 y n fi
xyT “ —x2y1 x2y2 ̈ ̈ ̈ x2yn P Rmˆn.
In MATLAB, simply type x*y’. *
3x1
1 x4 3 x 4
—–. . ... . fl xmy1 xmy2 ̈ ̈ ̈ xmyn
3
5
5
8
64/148
Elementwise Multiplication (.*)
‚ To multiply entries of two arrays of same size, element by element:
x .* y
x y
x .* y
65/148
Elementwise Division (./)
‚ To divide entries of an array by corresponding entries of another same-sized array:
x ./ y
‚ To divide a number by multiple numbers (speci ed by entries of an array):
s ./ y
‚ To divide all entries of an array by a common number:
x/s
x y
x ./ y
s
x
s ./ x
x s
x/ s
66/148
Elementwise Exponentiation (.^)
‚ To raise all entries of an array to (di erent) powers:
x .^ y
‚ Toraiseanumbertomultiplepowers(speci ed by entries of an array):
s .^ x
‚ To raise all entries of an array to a common power:
x .^ s
x y
x .^ y
s
x
s .^ x
x s
x .^ s
67/148
Mathematical Functions
‚ Built-in mathematical functions accept array inputs and return arrays of function evaluation, e.g.,
sqrt(A)
sin(A)
mod(A)
...
input
output
FUNCTION
68/148
Colon Operator
Supposea < b.
‚ To create an arithmetic progression from a to b (increment by 1):
a:b
Theresultisarowvector[a, a+1, a+2, ..., a+m],where
m =tb-au.
‚ To create an arithmetic progression from a to b with steps of size d > 0:
a:d:b
Theresultisarowvector[a, a+d, a+2*d, ..., a+m*d],where m =t(b-a)/du.
69/148
LINSPACE and LOGSPACE
‚ To create a row vector of n numbers evenly spaced between a and b: linspace(a, b, n)
Theresultis[a, a+d, a+2*d, ..., b],where d = (b-a)/(n-1).
‚ To create a row vector of n numbers that are logarithmically evenly spaced between 10a and 10b:
logspace(a, b, n)
The result is r10a, 10a`d, 10a`2d, . . . , 10bs, where
d = (b-a)/(n-1).
70/148
ZEROS, ONES, and EYE
‚ To create an pm ˆ nq zero matrix: zeros(m, n)
‚ To create an pm ˆ nq matrix all whose entries are one:
ones(m, n)
‚ To create the pm ˆ mq identity matrix: eye(m)
zeros(1,4)
ones(3,1)
0
0
0
0
1
1
1
eye(3)
1
0
0
0
1
0
0
0
1
71/148
Random Arrays
Each of the following generates an pm ˆ nq array of random numbers: ‚ rand(m,n): uniform random numbers in p0, 1q
‚ randi(k,m,n): uniform random integers in r1, ks
‚ randn(m,n): Gaussian random numbers with mean 0 and standard deviation 1
72/148
Random Arrays (Application)
To generate an pm ˆ nq array of
‚ uniform random numbers in pa, bq:
a + (b - a)*rand(m, n)
‚ uniform random integers in rk1, k2s:
randi([k1, k2], m, n)
‚ Gaussian random numbers with mean μ and standard deviation ‡: mu + sig*randn(m, n)
73/148
Concatenation
If two arrays A and B have comparable sizes, we can concatenate them.
‚ horizontally by [A B] AB
‚ vertically by [A; B]
A [A; B]
B
[A B]
74/148
RESHAPE and REPMAT
‚ reshape(A, m, n) reshapes the array A into an m ˆ n matrix whose elements are taken columnwise from A.
‚ repmat(A, m, n) replicates the array A, m times vertically and n times horizontally.
repmat(A,2,3)
123
A
13 242424
24 131313
2
A
reshape(A,3,2)
1
4
2
5
3
6
131313
1
3
5
2
4
6
1
242424
75/148
FLIP
‚ Type help flip on the Command Window and learn about flip function.
‚ Do the same with its two variants, flipud and fliplr
76/148
Creating Diagonal Matrices
‚ To create a diagonal matrix
»—v1 0 0 ̈ ̈ ̈ 0fi —0 v2 0 ̈ ̈ ̈ 0 —0 0 v3 ̈ ̈ ̈ 0 :
—– . . . ... . fl 0 0 0 ̈ ̈ ̈ vn
diag(v) Note.
‚ diag(v,k) puts the elements of v on the k-th super-diagonal. ‚ diag(v,-k) puts the elements of v on the k-th sub-diagonal.
v diag(v)
1
8
3
6
2
1
0
0
0
0
0
8
0
0
0
0
0
3
0
0
0
0
0
6
0
0
0
0
0
2
77/148
Extracting Diagonal Elements
Use diag(A,k) to extract the k-th diagonal of A. 1
‚ k > 0 for super-diagonals: 012
12345
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
diag(A, 1) diag(A, 2)
‚ k < 0 for sub-diagonals:
A
diag(A)
-2 12345 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
1 2 3
7 8 9 13 14 15 19 20
25
1 6 11
7 12 17 13 18 23 19 24
25
0 -1
A
diag(A)
diag(A,-1) diag(A,-2)
1 diag(A) short for diag(A,0) .
78/148
...
...
Using Vectors as Indices
To get/set multiple elements of an array at once, use vector indices.
‚ To grab 3rd, 4th, and 5th elements of x:
x(3:5) % or x([3 4 5])
‚ To grab 3rd to 8th elements of x:
x(3:8)
x(3:end)
‚ To grab 3rd to 7th elements of x:
x(3:7)
x(3:end-1)
x(3:5)
1
3
5
7
9
3
5
7
x(3:8) or x(3:end)
1
3
5
7
9
3
5
7
x(3:7) or x(3:end-1)
1
3
5
7
9
3
5
7
79/148
Using Vectors as Indices – Example
‚ To extract 2nd, 3rd, and 4th columns of the 2nd row of A:
A(2,2:4) % or A(2,[2 3 4]) ‚ To extract the entire 2nd row of A:
A(2,1:5)
A(2,1:end)
A(2,:)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
80/148
Using Vectors as Indices – Example
‚ To extract 2nd through 5th elements of the 4th column of A:
A([2 3 4 5],4)
A(2:5,4)
A(2:end,4)
‚ To extract the entire 4th column of A:
A(1:5,4)
A(1:end,4)
A(:,4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
81/148
Using Vectors as Indices – Example
‚ To grab the interior block of A: A(2:4,2:4)
A(2:end-1,2:end-1)
‚ To extract every other elements on every
other rows as shown:
A(1:2:5,1:2:5)
A(1:2:end,1:2:end)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
82/148
More on Arrays
83/148
Arithmetic Progressions
Question
Create the following periodic arithmetic progressions using ONE MATLAB
statement.
p1,2,3,4,0,1,2,3,4,0,1,2,3,4,0q.
m = 5;
n = 15;
mod([1:n], m)
84/148
Exercise: Arithmetic Progressions
Question
Create each of the following row vectors using ONE MATLAB statement. ‚ v “ p1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0q
‚ w “ p1,2,3,4,1,2,3,4,1,2,3,4q
85/148
Geometric and Other Progressions
Question
Create each of the following column vectors using ONE MATLAB statement. ‚ v “ p1,2,4,8,...,1024qT
‚ w “ p1,4,9,16,...,100qT Using the colon operator:
v = ( 2.^[0:10] )’
w = ( [1:10].^2 )’
Using the linspace function:
v = ( 2.^linspace(0, 10, 11) )’
w = ( linspace(1, 10, 10).^2 )’
86/148
Function Evaluation
Recallthatmathematicalfunctionssuchassin, sind, log, expaccept array inputs and return arrays of function evaluation.
Question
Create each of the the following row vectors using ONE MATLAB statement. ‚ u “ p1!,2!,3!,...,n!q
‚ v “ psin0 ̋,sin30 ̋,sin60 ̋,...,sin180 ̋q ‚ w “ pe1,e4,e9,...,e64q
v = sind(0:30:180)
w = exp([1:8].^2)
87/148
Matrices with Patterns
Question
Generate each of the following matrices using ONE MATLAB statement.
»1 2 3 4fi
—5 6 7 8 A “ – 9 10 11 12fl ,
13 14 15 16
A = reshape(1:16, 4, 4)’
B = ((1:10)’).^(1:10)
»—1 12 13
—2 22 23 B “ —– . . .
10 102 103
̈ ̈ ̈ 110fi ̈ ̈ ̈ 210
. . . . fl . ̈ ̈ ̈ 1010
88/148
Matrices with Patterns
Question
Suppose n is already stored in MATLAB. Generate each of the following matrices using ONE MATLAB statement. All the elements not shown are 0’s.
»2
— 4 C“—
6
fi »—cos1 ́3 fi — cos2 ́3 ,D“— cos3 ́3 .
— –
.. — ... ...
. 2 n fl —–
c o s p n ́ 1 q ́ 3 fl cos n
C = diag(2:2:2*n)
D = diag(cos(1:n)) - 3*diag(ones(n-1,1), 1)
89/148
Data Manipulation Functions
There are a number of MATLAB functions with spreadsheet functionalities that are suitable for data manipulation.
‚ max and min
‚ sum and prod
‚ cumsum and cumprod (cumulative sum and product) ‚ diff
‚ mean, std, and var (simple statistics)
‚ sort
90/148
Example 1: Finding the Maximum Value of a Vector
Question
Write a program to nd the maximum value of a vector.
‚ With loops:
% input: x
% output: m % DON’T USE max FOR THE VARIABLE NAME
m = x(1); % CODE ABORTS IF THE VECTOR IS EMPTY
for r = 2:length(x)
if m < x(r)
m = x(r);
end end
‚ Vectorized code: m = max(x)
91/148
Example 1 (cont’)
Question
Now modify the previous program to nd both the maximum value of a vector and the corresponding index.
‚ With loops:
% input: x
% output: m, index_m
m = x(1);
index_m = 1;
for r = 2:length(x)
if m < x(r)
m = x(r);
index_m = r;
end end
‚ Vectorized code:
[m, index_m] = max(x)
92/148
Example 2: Summing Elements in a Vector
Question
Sum all elements in a vector.
‚ With loops:
% input: x
% output: s % DON’T USE sum FOR THE VARIABLE NAME
s = 0; % s begins before the first iteration
for el = 1:length(x)
s = s + x(el);
end
‚ Vectorized code: s = sum(x)
93/148
FIND Function
Basic Usage of FIND
Let v be an array of numbers (can be a vector or a matrix). Then
find(
returns the (linear) indices of v satisfying
‚ Some examples of
‚ v > korv >= k ‚ v < korv <= k ‚ v == korv ~= k
‚ To combine more than two conditions:
‚ &(and) ‚ |(or)
94/148
Example 3: Comparing Elements in Vectors
Question
Compare two real vectors of the same length, say x and y, elementwise and determine how many elements of the former are larger than the latter.
‚ With loops:
% input: x, y
% output: nr_gt
nr_gt = 0;
for k = 1:length(x)
if x(k) > y(k)
nr_gt = nr_gt + 1;
end end
‚ Vectorized code:
nr_gt = length(find(x > y))
95/148
CPU Time
cputime reads total CPU time used by MATLAB from the time it was started. ‚ Single measurement:
ct = cputime; % total cputime as of now
t = cputime - ct; ‚ Average CPU time:
ct = cputime; % total cputime as of now
for i = 1:nr_reps
end
t_avg = (cputime - ct)/nr_reps;
96/148
Elapsed Time
At the execution of tic, MATLAB records the internal time (in seconds); at toc command, MATLAB displays the elapsed time.
‚ Single measurement:
tic % starts a stopwatch timer
toc % reads the elapsed time from tic ‚ Average elapse time:
tic % starts a stopwatch timer
for i = 1:nr_reps
end
t_avg = toc/nr_reps;
97/148
What Do You Think It Does?
Below is a modi ed version of an example code from MATLAB’s Help documentation for tic. What do you think it’s doing?
REPS = 1000; minTime = Inf; nSum = 10;
tic;
for i = 1:REPS
tStart = tic;
s = 0;
for j = 1:nsum
s = s + besselj(j,REPS);
end
tElpased = toc(tStart);
minTime = min(tElapsed, minTime);
end
t_avg = toc/REPS;
98/148
Example 4: Timing Elementwise Operations
Question
Generate a 107 ˆ 1 random vector and measure the internal time and CPU time when computing elementwise squares.
n = 1e7;
x = rand(n, 1);
t = cputime;
x1 = x.^2;
time1 = cputime - t;
tic
x2 = x.^2;
time2 = toc();
disp([time1, time2])
99/148
Pythagorean Triples
Question
Given n P N, nd all triples pa, b, cqN3 , with a, b § n, satisfying a2 ` b2 “ c2.
Notation.
‚ N:thesetofallnaturalnumbers,1,2,3,.... ‚ Nr1,ns “ t1,2,...,nu.
‚ N3 “ tpa,b,cq | a,b,c P Nu.
100 / 148
Pythagorean Triples – Solution Using Loops
% input: n
% output: M
iM = 0;
M = [];
for a = 1:n
for b = 1:n
c = sqrt(a^2 + b^2);
if mod(c, 1) == 0
iM = iM + 1;
M(iM, :) = [a, b, c];
end end
end
101 / 148
Pythagorean Triples – Solution Without Loops
% input: x
% output: M
A = repmat([1:n], n, 1);
B = repmat([1:n]’, 1, n);
C = sqrt(A.^2 + B.^2);
M = [A(:), B(:), C(:)];
lM = ( mod(M(:, 3), 1) ~= 0 );
M(lM, :) = [];
102 / 148
Birthday Problem
Question
In a group of n randomly chosen people, what is the probability that everyone has a di erent birthday?
1 Find this probability by hand.
2 Let n “ 30. Write a script that generates a group of n people randomly
and determines if there are any matches.
3 Modify the script above to run a number of simulations and numerically calculate the sought-after probability. Try 1000, 10000, and 100000 simulations. Compare the result with the analytical calculation done in 1.
103 / 148
Birthday Problem (Hints)
‚ For simplicity, ignore leap years.
‚ Create a random (column) vector whose elements represent birthdays of
individuals (denoted by integers between 1 and 365).
‚ Line up the birthdays in order and take the di erence of successive pairs. What does the resulting vector tell you?
‚ For 3, to run simulation multiple times, consider creating a random matrix whose rows represent birthdays of individuals and the columns correspond to di erent simulations.
104 / 148
Graphics in MATLAB
105 / 148
Anonymous Functions
Mathematical functions such as
f1pxq “ cos x sin`cosptan xq ̆, f2p◊q “ pcos 3◊ ` 2 cos 2◊q2,
f3px,yq“ sinpx`yq, 1`x2 `y2
can be de ned in MATLAB using anonymous functions:
f1 = @(x) cos(x).*sin(cos(tan(x)));
f2 = @(th) ( cos(3*th) + 2*cos(2*th) ).^2;
f3 = @(x, y) sin(x + y)./(1 + x.^2 + y.^2);
106 / 148
Anonymous Functions – Syntax
Take a closer look at one of them.
f1 = @(x) cos(x).*sin(cos(tan(x)));
‚ f1 : the function name or the function handle
‚ @ : marks the beginning of an anonymous function
‚ (x) : denotes the function (input) argument
‚ cos(x).*sin(cos(tan(x))) : MATLAB expression de ning f1pxq
107 / 148
Examples
Expressions in function de nitions can get very complicated. For example,
h1 = @(x) [2*x, sin(x)];
h2 = @(x) [2*x, sin(x); 5*x, cos(x); 10*x, tan(x)];
r = @(a,b,m,n) a + (b-a)*rand(m,n);
108 / 148
Exercise: Di erent Ways of De ning a Function
The function
can be de ned in two di erent ways:
f4p◊; c1, c2, k1, k2q “ pc1 cos k1◊ ` c2 cos k2◊q2
f4s = @(th,c1,c2,k1,k2) c1*cos(k1*th) + c2*cos(k2*th)
f4v = @(th,c,k) c(1)*cos(k(1)*th) + c(2)*cos(k(2)*th)
Question
Use f4s and f4v to de ne yet another anonymous functions for ‚ gp◊q“3cosp2◊q ́2cosp3◊q
‚ hp◊q “ 3 cosp◊{7q ` cosp◊q
109 / 148
Exercise: Understanding Anonymous Functions
Type in the following statements in MATLAB:
f1 = @(x) cos(x).*sin(cos(tan(x)));
f2 = @(th) ( cos(3*th) + 2*cos(2*th) ).^2;
x1 = 5; y1 = f1(x1)
x2 = [5:-2:1]; y2 = f1(x2)
TH = diag(0:pi/2:2*pi); R = f2(TH)
Question
1 What are the types of the input and output variables?
‚ x1andy1 ‚ x2andy2 ‚ THandR
2 Which of the three outputs will be a ected if elementwise operations were not used in the de nition of f1 and f2?
110 / 148
The PLOT Function
To draw a curve in MATLAB:
‚ Construct a pair of n-vectors x and y corresponding to the set of data points tpxi,yiq | i “ 1,2,...,nu which are to appear on the curve in that order.
‚ Then type plot(x, y) . For example:
x = linspace(0, 2*pi, 101);
y = sin(x);
plot(x, y)
or
% or simply plot(x, sin(x))
f = @(x) 1 + sin(2*x) + cos(4*x); % anonymous function
x = linspace(0, 2*pi, 101);
plot(x, f(x))
111 / 148
Example: Wiggly Curve
1 2 3 4
First, run the following script.
f1 = @(x) cos(x).*sin(cos(tan(x)));
x = 2*pi*[0:.0001:1]; % or x = linspace(0, 2*pi, 10001);
plot(x, f1(x))
shg
Play Around!
Observe what happens after applying the following modi cations one by one.
‚ Change line 3 into plot(x, f1(x), ’r’).
‚ Change line 3 into plot(x, f1(x), ’r--’).
‚ After line 3, add axis equal, axis tight.
‚ Then add text(4.6, -0.3, ’very wiggly’).
‚ Then add
xlabel(’x axis’), ylabel(’y axis’), title(’A wiggly curve’).
112 / 148
Note: Line Properties
‚ To specify line properties such as colors, markers, and styles:
plot(x, y, ’--’)
plot(x, y, ’g:’)
plot(x, y, ’.’, ’MarkerSize’, 3)
plot(x, y, ’b-’, ’LineWidth’, 5)
% dashed line
% dotted line in green
% adjust marker size
% adjust line width
Line Styles
- solid
: dotted -. dashdot -- dashed
Colors Markers
b blue
g green
r red
c cyan
m magenta y yellow
k black
w white
. o x +
* s
d
point circle x-mark plus
star square diamond
113 / 148
Note: Labels and Saving
‚ To label the axes and the entire plot, add the following after plot statement:
xlabel(’x axis’)
ylabel(’y axis’)
title(’my awesome graph’)
‚ Save gures using print function. Multiple formats are supported. print -dpdf ’wiggly’
% or print(’-dpdf’, ’wiggly’)
print -djpeg ’wiggly’
% or print(’-djpeg’, ’wiggly’)
print -deps ’wiggly’
% or print(’-deps’, ’wiggly’)
[pdf]
[jpeg]
[eps]
114 / 148
Note: Drawing Multiple Figures
‚ To plot multiple curves:
plot(x1, y1, x2, y2, x3, y3, ...)
‚ To create a legend, add
legend(’first graph’, ’second graph’, ’third graph’, ...)
115 / 148
Note: Miscellaneous Commands
‚ shg: (show graph) to bring Figure Window to the front ‚ figure: to open a new blank gure window
‚ clf: (clear gure) to clear previously drawn gures
‚ axis equal: to put axes in equal scaling
‚ axis tight: to remove margins around graphs
‚ axis image: same as axis equal and axis tight ‚ grid on: to put light gray grid lines
116 / 148
Exercise
Question
Do the following:
‚ De ne f pxq “ x3 ` x as an anonymous function.
‚ Find f1 and f2 and de ne them as anonymous functions. ‚ Plot all three functions in one gure in the interval r ́1, 1s. ‚ Include labels and title in your plot.
‚ Add legend to the graph.
‚ Save the graph as a pdf le.
117 / 148
Multiple Figures – Stacking
To draw multiple curves in one plot window as in Figure 1: ‚ One liner:
plot(x1, y1, x2, y2, x3, y3)
‚ Or, add curves one at a time using hold
command.
plot(x1, y1)
hold on
plot(x2, y2)
plot(x3, y3)
‚ hold on: holds the current plot for further overlaying
‚ hold off: turns the hold o
Figure 1: Multiple curves in one plot window
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
118 / 148
Multiple Figures – Subplots
To plot multiple curves separately and arrange them as in Figure 2:
subplot(1,3,1)
plot(x1, y1)
subplot(1,3,2)
plot(x2, y2)
subplot(1,3,3)
plot(x3, y3)
subplot(m,n,p):
‚ m, n: determine grid dimensions
‚ p: determines grid is to be used
111
0.9 0.9 0.9
0.8 0.8 0.8
0.7 0.7 0.7
0.6 0.6 0.6
0.5 0.5 0.5
0.4 0.4 0.4
0.3 0.3 0.3
0.2 0.2 0.2
0.1 0.1 0.1
000
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Figure 2: Multiple plots in 1 ˆ 3 grids
119 / 148
Exercise: Multiple Figures
Do It Yourself
Generate Figures 1 and 2.
‚ Common: Generating sample points
x = linspace(0, 1, 101);
y1 = x.^2; y2 = x.^4; y3 = x.^6;
‚ Figure 1:
hold off
plot(x, y1, ’*’)
hold on
plot(x, y2, ’g:o’)
plot(x, y3, ’r-s’)
‚ Figure 2:
subplot(1, 3, 1)
plot(x, y1)
subplot(1, 3, 2)
plot(x, y2, ’g’)
subplot(1, 3, 3)
plot(x, y3, ’r’)
120 / 148
The POLAR Function
To draw the polar curve r “ fp◊q, for ◊ P ra, bs:
‚ Grab n sample points
tp◊i,riq | ri “fp◊iq, 1§i§nu on the curve and form vectors th and r.
‚ Then type polar(th, r) .
‚ For example, to plot
r “ f2p◊q “ pcos 3◊ ` 2 cos 2◊q2, for ◊ P r0, 2fis:
th = linspace(0, 2*pi, 361);
f2 = @(th) (cos(3*th) + 2*cos(2*th)).^2;
polar(th, f2(th));
120
10
8
6
60
90
150
30
4
2
180
0
210
330
240
300
270
121 / 148
Exercise: Drawing Polar Curves
Question
1 Draw the graph of two-petal leaf given by
r“fp◊q“1`sinp2◊q, ◊Pr0,2fis.
2 Draw the graphs of
r“fp◊ ́fi{4q, r“fp◊ ́fi{2q, r“fp◊ ́3fi{4q
on the same plotting window.
3 Does your gure make sense?
122 / 148
Curves and the PLOT3 Function
Curves in R3 are plotted in an analogous fashion.
‚ Grab n sample points
tpxi,yi,ziq | i “ 1,2,...,nu on the curve and form vectors x, y, and z.
‚ Then type plot3(x, y, z) .
‚ For example, to plot the helix given by the
parametrized equation
rptq “ x10 cosptq, 10 sinptq, ty, for t P r0, 10fis:
t = linspace(0, 10*pi, 1000);
plot3(10*cos(t), 10*sin(t), t);
123 / 148
Exercise: Corkscrew
Question
Modify the code to generate a corkscrew by putting the helix outside of an upside down cone.
Hint: Use rptq “ xt cosptq, t sinptq, ty.
124 / 148
Surfaces and the SURF Function
To plot the surface of z “ f px, yq on R “ ra, bs ˆ rc, ds:
‚ Collect samples points on the intervals ra, bs and rc, ds and form vectors x and y.
‚ Based on x and y, generate grid pointstpxi,yjq | i“ 1,2,...,m,j “ 1,2,...,nu on the domain R and separate coordinates into matrices X and Y using meshgrid.
‚ Type surf(X, Y, f(X,Y))
r ́3, 3s ˆ r ́2, 2s
Figure3:Graphofz “ 29px2 ́y2qon
125 / 148
Note: How MESHGRID Work
>> x = [1 2 3 4]; y = [5 6 7]; >> [X, Y] = meshgrid(x,y)
X=
1234 1234 1234
7 6 5
Y= 1234 5555
6666 7777
6
6
7
7
y567 x1234
INPUT
MESHGRID
replicate
5555 66 7 7
1234 1234 4
OUTPUT
Y
5555
X 6666
7777 1234
1234 1234
1
2
3
126 / 148
replicate
Example: Saddle
Question
Plot the saddle parametrized by
z “ x2 ́ y2 c a2 b2
for your choice of a, b, and c.
x = linspace(-3, 3, 13);
y = linspace(-2, 2, 9);
[X, Y] = meshgrid(x, y);
a = 1.5; b = 1.5; c = .5;
g2 = @(x,y) c*( x.^2 /a^2 - y.^2 /b^2);
surf(X, Y, g2(X,Y))
axis equal, box on
Figure 3 was generated using this code.
127 / 148
Example: Oblate Spheroid
The gure for Problem 5 of Homework 1 was generated by the following code.2
a = 1; b = 1.35; c = 1;
nr_th = 41; nr_ph = 31;
x = @(th, ph) a*cos(th).*sin(ph);
y = @(th, ph) b*sin(th).*sin(ph);
z = @(th, ph) c*cos(ph);
th = linspace(0, 2*pi, nr_th);
ph = linspace(0, pi, nr_ph);
[T, P] = meshgrid(th, ph);
surf(x(T,P), y(T,P), z(T,P))
colormap(winter)
axis equal, axis off, box off
Figure 4: Oblate spheroid.
128 / 148
Function M-File
129 / 148
Basics
‚ A function is a piece of code which ‚ performs a speci c task;
‚ has speci c input and output arguments;
‚ is encapsulated to be independent of the rest of the program.
‚ In MATLAB, functions are de ned in .m les just as scripts. ‚ The name of the le and that of the function must coincide.
130 / 148
Why Functions?
Two Important Principles
‚ A piece of code should only occur once in a program. If it tries to appear more times, put it into its own function.
‚ Each function should do precisely one task well. If it has to carry out a number of tasks, separate them into their own functions.
Bene ts of writing functions:
‚ elevated reasoning by hiding details ‚ facilitates top-down design
‚ ease of software management
131 / 148
Writing a Function
A function m- le must be written in a speci c manner. ‚ When there is no input/output argument
function myfun()
....
end % <-- optional
‚ Where there are multiple input and output arguments
function [out1, out2, out3] = myfun(in1, in2, in3, in4)
....
end % <-- optional
132 / 148
Calling a Function
‚ If the function m- le myfun.m is saved in your current working directory3, you can use it as you would use any other built-in functions:
% when no input/output argument is required
myfun
or
% multiple inputs/outputs
[out1, out2, out3] = myfun(in1, in2, in3, in4)
‚ When not all output arguments are needed:
out1 = myfun(in1, in2, in3, in4) % only 1st output
[~, ~, out3] = myfun(in1, in2, in3, in4) % only 3rd output Note that tilde ( ~ ) is used as a placeholder.
3The path function gives more exibility in this regard.
133 / 148
Practical Matters: Speci cation
‚ The comments written below function header statement
function ... = myfun( ... )
% MYFUN: this awesome function calculates ...
....
end
canbeeasilyaccessedfromthecommandlineusinghelp myfun.
‚ Use this feature to write the function speci cation such as its purpose,
usage, and syntax.
‚ Write a clear, complete, yet concise speci cation.
134 / 148
Properties of Function M-Files
‚ The variable names in a function are completely isolated from the calling code. Variables, other than outputs, used inside a function are unknown outside of the function; they are local variables.
‚ The input arguments can be modi ed within the body of the function; no change is made on any variables in the calling code. (“ pass by value ” as opposed to “ pass by reference ”)
‚ Unlike script m- les, you CANNOT execute a function which has input argumentsusingtheRUN icon.
135 / 148
Example: Understanding Local Variables
‚ The function on the right nds the maximum value of a vector and an index of the maximalelement.
‚ Run the following on the Command Window. Pay attention to m.
>> m = 33;
>> x = [1 9 2 8 3 7];
>> [M, iM] = mymax(x)
>> disp(m)
function [m, el] = mymax(x)
if isempty(x)
el = [];
m = [];
return
end
el = 1;
m = x(el);
for i = 2:length(x)
if m < x(i)
el = i;
m = x(el);
end end
end
136 / 148
Example: Understanding Pass-By-Value
Consider the function
fpxq“sin|x| e ́|x|2 `lnp1`|x|q 1`|x| 1`2|x|
written as a MATLAB function:
function y = funky(x)
x = abs(x); % x is redefined to be abs(x)
y = sin(x).*exp(-x)./(1 + x.^2)...
+ log(1 + x)./(1 + 2*x);
Con rm that the function does not a ect x in the calling routine:
>> x = [-3:3]’;
>> y = funky(x);
>> disp([x, y])
end
137 / 148
Quiz
Question
Script.
x = 2;
x = myfun(x);
y = 2*x
Function.
function y = myfun(x)
x = 2*x;
y = 2*x;
end
What is the output when the script on the left is run?
1 y = 2 2 y = 4 3 y = 8 4 y = 16
138 / 148
Description of Problem
Problem
Given m and ◊ (in degrees), draw m equilateral triangles so that for any two successive triangles,
‚ the vertices of the smaller triangle are lying on the sides of the larger;
‚ the angle between the two triangles is ◊ (degrees)
Figure 5: A spiral triangle with m “ 21 and ◊ “ 4.5 ̋.
139 / 148
Generation and Transformation of Polygons
Let pxj,yjq, j “ 1,2,...,n be the coordinates of vertices of an n-gon. ‚ To plot the polygon:
x = [x1 x2 ... xn x1]; % note: x1 and y1 are repeated at
y = [y1 y2 ... yn y1]; % end to enclose the polygon.
plot(x, y)
‚ To plot the polygon obtained by scaling the original by a factor of s: plot(s*x, s*y)
‚ To plot the polygon obtained by rotating the original by an angle theta (in degrees) about the origin:
V = [x; y];
R = [cosd(theta) -sind(theta);
sind(theta) cosd(theta)];
Vnew = R*V;
plot(Vnew(1,:), Vnew(2,:))
% all vertices
% 2-D rotation matrix
% rotated vertices
140 / 148
Special Case: Inscribed Regular Polygons
‚ The vertices of the regular n-gon inscribed in the unit circle can be found easily using trigonometry, e.g.,
pxj,yjq“ˆcos2fipj ́1q,sin2fipj ́1q ̇, j“1,...,n. nn
‚ Thus we can plot it by
theta = linspace(0, 360, n+1);
V = [cosd(theta); % 2-by-(n+1) matrix whose cols are
sind(theta)]; % coordinates of the vertices
plot(V(1,:), V(2,:)), axis equal
‚ To stretch and rotate the polygon:
Vnew = s*R*V;
plot(Vnew(1,:), Vnew(2,:)), axis equal
141 / 148
Scale to Spiral
To create the desired spiraling e ect, the scaling factor must be calculated carefully.
‚ Useful:
‚ Compute the scaling factor s:
30 ̋ s
1
◊
sin– “ sin— “ sin“ abc
142 / 148
Put All Together
m = 21; d_angle = 4.5;
th = linspace(0, 360, 4) + 90;
V = [cosd(th);
sind(th)];
C = colormap(hsv(m));
s = sind(150 - abs(d_angle))/sind(30);
R = [cosd(d_angle) -sind(d_angle);
sind(d_angle) cosd(d_angle)];
hold off
for i = 1:m
if i > 1
V = s*R*V;
end
plot(V(1,:), V(2,:), ’Color’, C(i,:))
hold on
end
set(gcf, ’Color’, ’w’)
axis equal, axis off
143 / 148
Exercise: Script to Function
Question
Modify the script so that it can create a spiral n-gon consisting of m regular n-gons successively rotated by ◊ degrees. Then turn the script into a function m- le spiralgon.m.
function V = spiralgon(n, m, d_angle)
% SPRIALGON plots spiraling regular n-gons
% input: %
%
%
n = the number of vertices
m = the number of regular n-gons
d_angle = the degree angle between successive n-gons
(can be positive or negative)
% output: V = the vertices of the outermost n-gon
% Fill in the rest.
....
144 / 148
Local Functions
There are “general purpose” functions which may be used by a number of other functions. In this case, we put it into a separate le. However, many functions are quite speci c and will only be used in one program. In such a case, we keep them in whatever le calls them.
‚ The rst function is called the primary function and any function which follows the primary function is called a local function or a subfunction.
function
end
function
....
end
‚ The primary function interact with local functions through the input and output arguments.
‚ Only the primary function is accessible from outside the le.
145 / 148
De ning/Evaluating Mathematical Functions in MATLAB
‚ using an anonymous function ‚ using a local function
‚ using a primary funciton
‚ using a nested function
‚ putting it in situ, i.e., write it out completely in place – not recommended See time_anon.m.
146 / 148
Passing Function to Another Function
‚ A function can be used as an input argument to another function.
‚ The function must be stored as a function-handle variable in order to be
passed as an argument.
‚ This is done by prepending “ @ ” to the function names. For instance, try
>> class(mymax) % WRONG
>> class(@mymax) % CORRECT
147 / 148
When a Function Is Called
Below is what happens internally when a function is called:
‚ MATLAB creates a new, local workspace for the function.
‚ MATLAB evaluates each of the input arguments in the calling statements (if there are any). These values are then assigned to the input arguments of the function. These are now local variables in the new workspace.
‚ MATLAB executes the body of the function. If there are output arguments, they are assigned values before the execution is complete.
‚ MATLAB deletes the local workspace, saving only the values in the output arguments. These values are then assigned to the output arguments in the calling statement.
‚ The code continues following the calling statement.
148 / 148