CS计算机代考程序代写 matlab Module 1: Programming in MATLAB

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=Y] Is X greater than or equal to Y? ‚ [X<=Y] Is X less than or equal to Y? ‚ [X==Y] Is X equal to Y? ‚ [X~=Y] Is X not equal to Y? The symbols used between X and Y are called the relational operators. 5/148 Logical Variables and Logical Operators ‚ A relational statement evaluates to either True(1) or False(0); these are called logical variables or boolean variables. ‚ As arithmetic operators (+,-,*,/) put together two numbers and produce other numbers, logical operators combine two logical variables to produce other logical variables. ‚ Logical Operators: and, or, not, xor 6/148 Logical Operator: && (AND) Let A and B be two logical variables. The && operation is completely dened by the following truth table: A B A && B FFF FTF TFF TTT NotethatA && BistrueifandonlyifbothAandBaretrue. 7/148 Logical Operator: || (OR) Let A and B be two logical variables. The || operation is completely dened by the following truth table: A B A||B FFF FTT TFT TTT NotethatA || BisfalseifandonlyifbothAandBarefalse. 8/148 Logical Operator: xor (exclusive or) This is a special variant of the || operator. A B xor(A,B) FFF FTT TFT TTF Note that xorg(A,B) is true if only one of A or B is true. 9/148 Logical Operator: ~ (NOT) This is a negation operator. A ~A FT TF 10/148 Combination of Logical Operations LetAandBbelogicalvariables.Then~(A && B)and~A || ~Bare equivalent: A B A&&B ~(A&&B) FF FT TF TT A B ~A ~B ~A||~B FF FT TF TT 11/148 Quadratics Revisited Consider a monic quadratic function qpxq “ x2 ` bx ` c on a close interval rL, Rs. ‚ Critical point: xc “ ́b{2 ‚ If xc P pL, Rq, qpxq attains the (global) minimum at xc; otherwise, the minimum occurs at one of the endpoints x “ L or x “ R. Question Write a program which determines whether the critical point of qpxq falls on the interval. 12/148 Initialization b = input(’Enter b: ’); c = input(’Enter c: ’); L = input(’Enter L: ’); R = input(’Enter R (L= R’)
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 = 1 :
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 modied within the loop body.
‚ However, it is a recommended practice that the value of the loop variable is never modied 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 =1:
‚ 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 dierent ways:
szA = size(A)
[m, n] = size(A)
Q. How are they dierent?
‚ 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 specic 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 (specied 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 (dierent) powers:
x .^ y
‚ Toraiseanumbertomultiplepowers(specied 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 modied 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 dierent 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 dierence 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 dierent 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 dened 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 dening f1pxq
107 / 148

Examples
Expressions in function denitions 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: Dierent Ways of Dening a Function
The function
can be dened in two dierent 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 dene 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 aected if elementwise operations were not used in the denition 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 modications 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:
‚ Dene f pxq “ x3 ` x as an anonymous function.
‚ Find f1 and f2 and dene 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 specic task;
‚ has specic input and output arguments;
‚ is encapsulated to be independent of the rest of the program.
‚ In MATLAB, functions are dened 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.
Benets 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 specic 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: Specication ‚ The comments written below function header statement function ... = myfun( ... ) % MYFUN: this awesome function calculates ... .... end canbeeasilyaccessedfromthecommandlineusinghelp myfun. ‚ Use this feature to write the function specication such as its purpose, usage, and syntax. ‚ Write a clear, complete, yet concise specication. 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 modied 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); Conrm that the function does not aect 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 eect, 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 specic 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

Dening/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