Programming with
MATLAB
Cheng-Liang Chen
PSE
LABORATORY
Department of Chemical Engineering National TAIWAN University
Chen CL
1
Relational Operators
Six relational operators Rela Re lati tion onal al oper operat ato or
Mean Me anin ing g
< <= > >= ==
Less than. Less than or equal to. Greater Great er than. Greater Great er than or equal to. Equal to.
∼=
Not equal to.
Note: single = is assignment or replacement
Chen CL
1
Relational Operators
Six relational operators Rela Re lati tion onal al oper operat ato or
Mean Me anin ing g
< <= > >= ==
Less than. Less than or equal to. Greater Great er than. Greater Great er than or equal to. Equal to.
∼=
Not equal to.
Note: single = is assignment or replacement
Chen CL
Result of comparison: (can be used as a variable) 1 if the comparison is true 0 if the comparison is false x y z w
2
= = = =
2; 5; x < y,... x > y
z = 1 w = 0
Compare arrays x y z p q
= = = = =
[ 6, 3, 9]; [14, 2, 9]; (x < y) (x y ),... (x ~= y),... (x > 8)
z = 1 p = 1 q = 0
0
0
z = x( x
1
0
z = 6
0
1
Chen CL
3
Logical Operators and Functions
Three logical operators Operator
Name
∼
NOT
Definition
∼ A returns an array the same dimension as
A; the new array has ones where A is zero and zeros where A is nonzero. &
AND A&B returns an array the same dimension as A and B; the new array has ones where both A and B have nonzero elements and zeros where either A or B is zero.
|
OR A B returns an array the same dimension as A and B; the new array has ones where at least one element in A or B is nonzero and zeros where A and B are both zero.
|
Chen CL
4
Element-by-element operations (except NOT) x y z p q r %
= [ 6, 3, 9]; = [14, 2, 9]; = ~x,... = x > y = x <= y =~(x > y) = (x <= y)
z = 0&3 w = 2&3 z = 0 w = 1
z =
0
0
0
0
1
0
1
0
1
1
0
1
p = q = r =
z = [5,-3,0,0]&[2,4,0,5] z = 1
1
0
0
Chen CL
5
Order of Precedence for Operator Types First Second Third Fourth Fifth
Parentheses ; evaluated starting with the innermost pair. Arithmetic operators and NOT ( ); evaluated from left to right. Relational operators; evaluated from left to right. Logical AND . Logical OR.
z = 5 < 6 & 1 % (5<6)&1 z = 1 x y a z % w % z
= [ 6, 3, 9]; = [14, 2, 9]; = [ 4, 3, 12]; = x > y & a ,... = (x>y) & a = x > y & x > a = (x>y) & (x>a) = 0 1 0 w = 0 0 0
∼
~(4 > 5) is equivalent to 5 >= 4 5 < x < 10 is equivalent to (5
= = = =
[5,-3,0,0] | [2,4,0,5] 3<5|4==7 %(3<5)|(4==7) 1|0&0 % 1|(0&0) 1&0|0 % (1&0)|0
z = 1 w = 1 p = 1 q = 0
1
0
1
Chen CL
6
z = ~3==7|4==6 %((~3)==7)|(4==6) z = 0 x y z x y z
x = [1, 1, 0, 0]’; = xor([3,0,6],[5,0,0]) y = [1, 0, 1, 0]’; = [3,0,6]&[5,0,0] Truth_table = [x, y, ~x, x|y, x&y, = [3,0,6]|[5,0,0] Truth_table = 1 1 0 1 1 = 1 0 0 1 0 0 0 1 0 1 1 1 0 = 0 0 1 0 0 1 0 0 = 1 0 1
x
y
true
true false true
true
false false true false true
false true
∼x
xy
|
x&y xor(x, y) true false
true
true false true
false false true
false false false
xor(x,y)] 0 1 1 0
Chen CL
7
Logical Operators and find Function Logical function
Definition
any(X)
Returns a scalar, which is 1 if any of the elements in the vector x is nonzero and 0 otherwise. Returns a row vector having the same number of column as A and containing ones and zeros, depending on whether or not the corresponding column of the matrix A contains any nonzero elements. Returns a scalar, which is 1 if all the elements in the vector x are nonzero and 0 otherwise. Returns a row vector having the same number of column as the matrix A and containing ones and zeros, depending on whether or not the corresponding column of A has all nonzero elements. Computes an array containing the indices of the nonzero elements of the array A.
any(A)
all(X) all(A)
find(A)
Chen CL
8
Logical function
Definition
[u,v,w] = find(A)
Computes the arrays u and v containing the row and column indices of the nonzero elements of the array A and computes the array w containing the values of the nonzero elements. The array w may be omitted. Returns an array of the same dimension as A with ones where the elements of A are finite and zeros elsewhere. Returns an array of the same dimension as A with ones where A has ’NaN’ and zeros elsewhere. (’NaN’ stands for ”Not a Number,” which means an undefined result.) Returns an array of the same dimension as A, with ones where A has ’inf ’ and zeros elsewhere. Returns a 1 if A is an empty matrix and 0 otherwise. Returns a 1 if A has no elements with imaginary parts and 0 otherwise. Returns an array the same dimension as A and B; the new array has ones where either A or B is nonzero, but not both, and zeros where A and B are either both nonzero or both zero.
finite(A) isnan(A)
isinf(A) isempty(A) isreal(A) xor(A,B)
Chen CL
9
Logical Operators and find Function Example L_nonzero = x = [5,-3,0,0,8]; 1 2 y = [2, 4,0,5,7]; L_nonzero = find(x&y) ,... values = 2 4 values = y(x&y) ,... how_many = length(values) how_many = 3
5 7
Chen CL
10
Test Your Understanding
T4.2-1 If x = [5, 3, 18, 4] and y = [ 9, 13, 7, 4], what will be the result of the following operations ? Use MATLAB to check your answer.
−
1. 2. 3. 4.
z z z z
= y>x = x&y =xy = xor(x, y)
∼
|
−
Chen CL
11
Example: Height and Speed of A Projectile The height and speed of a projectile (such as a thrown ball) launched with a speed of v0 at an angle A to the horizontal are given by h(t) = v 0tsinA 0.5gt 2, v(t) = v02 2v0gtsinA + g 2t2, where g is the acceleration due to gravity.
−
−
The projectile will strike the ground when h(t) = 0, which gives the time to hit, thit = 2(v0/g) sin(A). Suppose that A = 40o, v0 = 20 meters/second, and g = 9.81 meters/second2. Use the MATLAB relational and logical operators to find the times when the height is no less than 6 meters and the speed is simultaneously no greater than 16 meters/second. In addition, discuss another approach to obtaining a solution.
Solution: The key to solving this problem with relational and logical operators is to use the
find command to determine the times at which the logical expression (h>=6)&(v<=16) is true
Chen CL
12
Example: Height and Speed of A Projectile % Set the values for initial speed, gravity, and angle. v0 = 20; g = 9.81; A=40*pi/180; % Compute the time to hit. t_hit = 2*v0*sin(A)/g; % Comput the arrays containing time, height, and speed. t = [0 : t_hit/100 : t_hit]; h = v0*t*sin(A) - 0.5*g*t.^2; v = sqrt(v0^2 - 2*v0*g*sin(A)*t + g^2*t.^2); h6=6*t.^0; v16=16*t.^0; subplot(3,1,1) plot(t,h,t,h6,’Linewidth’,2) set(gca,’LineWidth’,2,’FontSize’,12) xlabel(’\bf Time’,’FontSize’,12) ylabel(’\bf Height’,’FontSize’,12) subplot(3,1,2) plot(t,v,t,v16,’Linewidth’,2) set(gca,’LineWidth’,2,’FontSize’,12) xlabel(’\bf Time’,’FontSize’,12) ylabel(’\bf Velocity’,’FontSize’,12)
Chen CL
13
Example: Height and Speed of A Projectile % Determine when the height is no less than 6, % and the speed is no greater than 16. y = h>=6 & v<=16; u = find( y ); % u = find( h>=6 & v<=16 ); subplot(3,1,3) plot(t,y,’o’,’MarkerSize’,4) set(gca,’LineWidth’,2,’FontSize’,12) xlabel(’\bf Time’,’FontSize’,12) ylabel(’\bf T (1) or F (0)’,’FontSize’,12) % Compute the corresponding times. t_1 = (u(1)-1)*(t_hit/100) t_2 = u(length(u)-1)*(t_hit/100)
Chen CL
14
Example: Height and Speed of A Projectile
t_1 = 0.8649 t_2 = 1.7560
Chen CL
15
Test Your Understanding T4.2-2 Consider the problem given in the previous Example. Use relational and logical operators to find the times for which either the projectile’s height is less than 4 meters or the speed is greater than 17 meters per second. Plot h(t) and v(t) to confirm your answer.
Chen CL
16
Conditional Statements % Decision making (I) If I get a raise, I will buy a new car . (period) % Decision making (II) If I get at least a $100 per week raise, I will buy a new car; else, I will put the raise into savings . (period) % Decision making (III) If I get at least a $100 per week raise, I will buy a new car; else, if the raise is greater than $50, I will buy a new stereo; otherwise, I will put the raise into savings
Chen CL
17
The if Statement
if logical expression statements end
if x >= 0 y = sqrt(x) end
if x >= 0, y = sqrt(x), end
Chen CL
18
Nested if Statements
if logical expression 1 if logical expression 2 statements end end
⇓ if (logical expression 1) & (logical expression 2) statements end
Chen CL
z = w =
19
√ x + √ y log x − 3log y
z = 0;w = 0; if (x >= 0) if (y >= 0) z = sqrt(x) + sqrt(y) w = log(x) - 3*log(y) end %%% Nested if statements end if logical expression 1 z = 0;w = 0; statement group 1 if (x >= 0)&(y >= 0) if logical expression 2 z = sqrt(x) + sqrt(y) statement group 2 w = log(x) - 3*log(y) end end end
Chen CL
20
The else Statement
if logical expression statement group 1 else statement group 2 end
√ x y= e −1 x
for x
≥0
for x < 0
if x >= 0 y = sqrt(x) else y = exp(x) - 1 end
Chen CL
21
Note: For logical expression of an array, the test returns a value of true only if ALL the elements of the logical expression are nonzero ! x = [4, -9, 25]; if x < 0 disp(’All elements of’) disp(’x are negative.’) else y = sqrt(x) end
x = [4, -9, 25]; if x >= 0 y = sqrt(x) else disp(’Some of the elements’) disp(’ of x are negative.’) end
y = 2
Some of the elements of x are negative.
0+3.000i
5
Chen CL
22
The elseif Statement if logical expression 1 statement group 1 elseif logical expression 2 statement group 2 else statement group 3 end
Chen CL
23
ln x y= √ x
if x if
ln√ x y= x e −1 x
≥ 0≤
if x >= 5 y = log(x) else if x >= 0 5 y = sqrt(x) x < 5 end end % Note : no actio n
if x > 10 if 0
≤ x ≤ 10
if x < 0
if x >= 5 y = log(x) elseif x >= 0 y = sqrt(x) end for x=-2 ?
if x > 10 y = log(x ) elseif x >= 0 y = sqrt(x (x) ) else y = exp(x) - 1 end
Chen CL
24
if x > 10 y = lo g(x) if y >= 3 z = 4*y elseif y >= 2.5 z = 2*y else z = 0 end else y = 5*x z = 7*x end
Chen CL
25
Test Your Understanding
T4.3-1 Enter the script file shown above, whose flowchart is show in Figure 4.3-6. Run the file for the following values of x. x . Check the program results by hand: x = 11 11,, 25 25,, 2, 13 13..
T4.3-2 Giv iveen a nu num mbe berr x and and th thee qu quad adra rant nt q (q = 1, 2, 3, 4) 4),, write a program to compute sin−1(x) in degrees, taking into account the quadrant. The program should display an error message if x > 1.
| | |
Chen CL
26
Strings and Conditional Statements name = ’Leslie Student’ full_name = [name(1:6),’ C. ’,name(7:14)] name = Leslie Student full_name = number = ’123’ Leslie C. Student number = full_name(8) = ’F’ 123 full_name = size(name) Leslie F. Student ans = findstr(full_name,’e’) 1 14 ans = name(8) 2 6 15 ans = % letter ’e’ occurs in 2nd, S % 6th, and 15th columns first_name = name(1:6) first_name =
Chen CL
Input Prompts and Output Messages One of the most important applications for strings is to create Input Prompts and Output Messages x = input(’prompt’,’string’) response = input(’Do you want to continue? Y/N [Y]: ’,’s’); if ( isempty(response) | (response == ’Y’) | (response == ’y’) ) response = ’Y’ else response = ’N’ end
27
Chen CL
28
for Loops
for loop variable = m:s:n statements end
for k = 5:10:35 x = k^2 end
for x = 0:2:10, y = sqrt(x), end
Chen CL
29
Nested Loops function A = specmat(n) specmat(5) A = ones(n); ans = for r = 1:n 1 for c = 1:n 1 if (r>1)&(c>1) 1 s = A(r-1,c) + A(r,c-1); 1 if s<20 1 A(r,c) = s; else A(r,c) = max(A(r-1,c),A(r,c-1)); end end end end
1 2 3 4 5
1 3 6 10 15
1 4 10 10 15
1 5 15 15 15
Chen CL
30
Test Your Understanding T4.4-1 Write a program to produce the following matrix:
4 10 A= 16
8 14 20
12 18 24
22 26 30
Chen CL
31
Rules for Using for Loops m : s : n
The step value s may be negative. For example, k = 10 : 2 : 4 produces k = 10, 8, 6, 4.
−
If s is omitted, the step value defaults to one.
If s is positive, the loop will not be executed if m > n.
If s is negative, the loop will not be executed if m is less than n.
If m equals n, the loop will be executed only once.
If the step value s is not an integer, round-off errors can cause the loop to execute a different number of passes than intended.
Chen CL
32
break and continue Commands for k = 1:10 x = 50 - k^2; if x < 0 break end y = sqrt(x) end % The program execution % jumps to % ------HERE-----% if the break command % is executed
x = [10,1000,-10,100]; y = NaN*x; for k = 1:length(x) if x(k) < 0 continue end y(k) = log10(x(k)); end y y = 1
3
NaN
2
Chen CL
33
Vector as Loop Variables A = [1,2,3,4,5,6]; for v = A disp(v) end 1 2 3 4 5 6
A = [1,2,3,4,5,6]; n = 3; for k = 1:n v = A(:,k) end v = 1 v = 2 v = 3
Chen CL
34
Implied Loops x = [0:5:100]; y = cos(x);
y = find(x>0)
for k = 1:21 x = (k-1)*5; y(k) = cos(x); end
j=0; for i=1:length(x) if x(i)>0 j = j + 1; y(j)=i; end end
Chen CL
35
Test Your Understanding T4.4-1 Write a for loop that is equivalent to the command sum(A), where A is a matrix.
Chen CL
36
Example: Data Sorting A vector x has been obtained from measurements. Suppose we want to consider any data value in the range 0.1 < x < 0.1 as being erroneous. We want to remove all such elements and replace them with zeros at the end of the array. Develop two ways of doing this. An example is given in the following table.
−
x(1) x(2) x(3) x(4) x(5) x(6) x(7)
Solution:
Before
After
1.92 0.05 2.43 0.02 0.09 0.85 0.06
1.92 2.43 0.85 0.00 0.00 0.00 0.00
− − −
−
Chen CL
37
%%%%% Method 1: use a for loop %%%%% x = [1.92,0.05,-2.43,-0.02,0.09,0.85,-0.06]; y = [];z = []; for k = 1:length(x) if abs(x(k)) >= 0.1 y = [y,x(k)]; else z = [z,x(k)]; end end xnew = [y, zeros(size(z))] xnew = 1.9200
-2.4300
0.8500
0
0
0
0
%%%%% Method 2: use find function %%%%% x = [1.92,0.05,-2.43,-0.02,0.09,0.85,-0.06]; y = x(find(abs(x) >= 0.1)); z = zeros(size(find(abs(x)<0.1))); xnew = [y,z]
Chen CL
38
EX: Flight of An Instrumented Rocket All rockets lose weight as burn fuel; thus the mass of the system is variable. The following equations describe the speed v and height h of a rocket launched vertically, neglecting air resistance. They can be derived from Newton’s law [Beer and Johnston 1997].
0 v(t) = u ln mm gt 0 −qt u h(t) = (m0 qt) l n (m0 q
−
− qt) gt m u +u ln (m + 1) t − − ln(m ) 2 q −
2
0
0
0
where m0 is the rocket’s initial mass, q is the rate as which the rocket burns fuel mass (assumed constant), u is the exhaust velocity of the burned fuel relative to the rocket (also assumed constant), and g is the acceleration due to gravity. Let b be the burn time , after which all the fuel is consumed. Thus the rocket’s mass without fuel is me = m 0 qb. At the burn time , the rocket height is
−
h(b) =
u (me) l n (me) + u ln (m0 + 1) b q
−
gb 2 2
− mq u ln(m ) 0
0
Chen CL
39
For t > b the rocket engine no longer produces thrust, and the speed and height are given by v(t) = v(b)
− g(t − b) (for t > b) g(t − b) h(b) + v(b)(t − b) − 2 2
h(t) =
( )
∗
The time t p to reach the peak height is found by setting v(t) = 0. The result is t p = b + v(b)/g. Substituting this expression into the expression (*) for h(t) gives the following expression for the peak height: h p = h(b) + v 2(b)/(2g). The time at which the rocket hits the ground is thit
= t + 2h /g. p
p
Suppose the rocket is carrying instruments to study the upper atmosphere, and we need to determine the amount of time spent above 50, 000 feet as a function of the burn time b (and thus as a function of the fuel mass qb). Assume that we are given the following values: me = 100 slugs, q = 1 slug per second, u = 8000 feet per second, and g = 32.2 feet per second2. If the rocket’s maximum fuel load is 100 slugs, the maximum value of b is 100/q = 100 seconds. Write a MATLAB program to solve this problem.
Chen CL
40
Solution Pseudocode for developing the program appears in the following. A for loop is a logical choice to solve this problem because we know the burn time b and thit, the time it takes to hit the ground. A MATLAB program to solve this problem appears in the following. It has two nested for loops. The inner loop is over time and evaluates the equations of motion at times spaced 1/10 of a second apart. This loop calculates the duration above 50, 000 feet for a specific value of the burn time b. We can obtain more accuracy by using a smaller value of the time increment dt. The outer loop varies the burn time in integer values from b = 1 to b = 100. The final result is the vector of durations for the various burn time . Th following Figure gives the resulting plot.
Chen CL
41
Enter data. Increment burn time from 0 to 100. For each burn-time value: Compute m0, vb, hb, h p. If h p hdesired, Compute t p, thit.
≥
Increment time from 0 to thit. Compute height as a function of time, using the appropriate equation, depending on whether burnout has occurred. Compute the duration above desired height. End of the time loop. If h p < hdesired, set duration equal to zero. End of the burn-time loop. Plot the result.
Chen CL file rocket1.m: % Script % Computes flight duration as a function of burn time. % Basic data values. m_e = 100; q = 1; u = 8000; g = 32.2; dt = 0.1; h_desired = 50000; for b = 1:100 % Loop over burn time. burn_time(b) = b; % % Following lines implement the formulas in text. m_0 = m_e + q*b; v_b = u*log(m_0/m_e) - g*b; % h_b = ((u*m_e)/q)*log(m_e/(m_e+q*b))+u*b-0.5*g*b^2; h_b = ((u*m_e)/q)*log(m_e) + u*(log(m_0)+1)*b ... - 0.5*g*b^2 - ((u*m_0)/q)*log(m_0); h_p = h_b + v_b^2/(2*g); if h_p >= h_desired; % % Calculate only if peak height > desire height. t_p = b + v_b/g; %Compute peak time. t_hit = t_p + sqrt(2*h_p/g); %Compute time to hit. % for p = 0:t_hit/dt:t_hit for p = 0:t_hit/dt % % Use a loop to compute the height vector. k = p+1; t = p*dt; time(k) = t; % k = p + 1; t = p*dt; time(k) = t; if t <= b % % Burnout has not yet occurred. h(k) = (u/q)*(m_0 - q*t)*log(m_0 - q*t) ... + u*(log(m_0) + 1)*t - 0.5*g*t^2 ... (m_0*u/q)*log(m_0); else % % Burnout has occurred. h(k) = h_b + v_b*(t - b) - 0.5*g*(t - b)^2; end
42
%
% Compute the duration. duration(b) = ... length(find(h>=h_desired))*dt; else % % Rocket did not reach desired heig duration(b) = 0; end end % Plot the results. plot(burn_time,duration) ,... xlabel(’Burn Time (sec)’),... ylabel(’ Duration (sec)’),... title(’Duration Above 50,000 Feet’) Duration Above 50,000 Feet 180
160
140
120 ) c e 100 s ( n o i t a r u D
80
60
40
20
0
0
10
20
30
40
50 60 Burn Time (sec)
70
80
90
100
Chen CL
43
while Loops
The while loop is used when the looping process terminates because a specified condition is satisfied, and the number of passes in not known in advance. while logical expression statements end
A practical application of while loop is when we want the loop to continue as long as a certain statement in true
Chen CL
x = 5; while x < 25 disp(x) x = 2*x-1; end 5 9 17
44
x=1; while x ~= 5 disp(x) x = x+1; end 1 2 3 4
Chen CL
45
Example: Time to Reach a Specified Height Consider the variable-mass rocket treated in previous Example. Write a MATLAB program to determine how long it takes for the rocket to reach 40, 000 feet if the burn time is 50 seconds.
Solution: The pseudocode appears in the following Table. Because we do not know the time required, a while loop is convenient to use. The program performs the task and is a modification of the previous program. Note that the new program allows for the possibility that the rocket might not reach 40, 000 feet. It is important to write your program to handle all such foreseeable circumstances. The answer given by the program is 53 seconds.
Chen CL
46
Enter data. Compute m0, vb, hb, h p. If h p hdesired, Use a while loop to increment time and compute height until desired height is reached.
≥
Compute height as a function of time, using the appropriate equation, depending on whether burnout has occurred. End of the results. Display the results. If h p < hdesired, rocket cannot reach desired height.
Chen CL % Script file rocket2.m % Compute time to reach desired height. % Set the data values. h_desired = 40000; m_e = 100; q = 1; u = 8000; g = 32.2; dt = 0.1; b = 50; % Compute values at burnout, peak time, and height. m_0 = m_e + q*b; v_b = u*log(m_0/m_e) - g*b; h_b = ((u*m_e)/q)*log(m_e/(m_e+q*b))+u*b- .5*g*b^2; t_p = b + v_b/g; h_p = h_b + v_b^2/(2*g); % If h_p > h_desired, compute time to reached h_desired. if h_p > h_desired h = 0; k = 0; while h < h_desired % Compute h until h = h_desired. t = k*dt; k = k + 1; if t <= b h = (u/q)*(m_0 - q*t)*log(m_0 - q*t) ... % Burnout has not yet occurred. + u*(log(m_0) + 1)*t - 0.5*g*t^2 ... - (m_0*u/q)*log(m_0); else h = h_b - 0.5*g*(t - b)^2 + v_b*(t - b); % Burnout has occurred. end end % Display the results. disp(’The time to reach the desired height is:’) disp(t) else disp(’Rocket cannot achieve the desired height.’) end The time to reach the desired height is: 53
47
Chen CL
48
Test Your Understanding T4.4-3 Rewrite the following code using a while loop to avoid using the break command. for k = 1:10 x = 50 - k^2; if x < 0 break end y = sqrt(x) end
T4.4-4 Find to two decimal places the largest value of x before the error in the series approximation ex 1 + x + x2/2 + x3/6 exceeds 1 percent (answer: x = 0 83)
≈
Chen CL
49
The switch Structure The switch structure provides an alternative to using the if , elseif , and else commands. Anything programmed using switch can also be programmed using if structures. However, the switch can provide more readable codes. switch input expression (scalar or string) case value1 statement group 1 case value2 statement group 2 ... ... otherwise statement group n end
switch angle case 45 disp(’Northeast’) case 135 disp(’Southeast’) case 225 disp(’Southwest’) case 315 disp(’Northwest’) otherwise disp(’Direction Unknow end
Chen CL
Use a string variable for input expression to result in very readable programs
50
The switch statement can handle multiple conditions in a single case statement by enclosing the case value in a cell array
t = [0:100]; x = exp(-t).*sin(t); switch angle response = ... case {0,360} input(’Type min, max, or sum.’,’s’) disp(’North’) response = lower(’response’); case {-180,180} switch response disp(’South’) case min case {-270,90} minimum = min(x) disp(’East’) case max case {-90,270} maximum = max(x) disp(’West’) case sum otherwise total = sum(x) disp(’Direction Unknown’) otherwise end disp(’You have not entered a proper choice.’) end
Chen CL
51
Test Your Understanding T4.5-1 Write a program using the switch structure to input one angle, whose value may be 45, 45, 135 degrees, and display the quadrant (1, 2, 3, or 4) containing the angle.
−
Chen CL
52
Debugging MATLAB Programs function y = fun1(x) avg = sum(x)/length(x); y = fun2(avg,x); function above = fun2(x,avg) above = length(find(x>avg)); above = fun1([1,2,3,4,10]) above = 3
Chen CL
53
Debugging A Loop function z = invest(x,r) z = 0; y = 1+0.01*r; for k = 1:length(y) z = z*y+x(k); end >>total = invest([1000,1500,2000],10) total = 1000
Chen CL
54
Example: A College Enrollment Model (I) As a example of how simulation can be used for operations research, consider the following college enrollment model. A certain college wants to analyze the effect of admissions and freshman retention rate on the college’s enrollment so that it can predict the future need for instructors and other resources. Assume that the college has estimates of the percentages of students repeating a grade or leaving school before graduating. Develop a matrix equation on which base a simulation model that can help in this analysis.
Solution: Suppose that the current freshman enrollment is 500 students and the college decides to admit 1000 freshman per year from now on. The college estimates that 10 percent of the freshman class will repeat the year. The number of freshman in the following year will be 0.1(500) + 1000 = 1050, then it will be 0.1(1050) + 1000 = 1105, and so on. Let x1(k) be the number of freshman in year k, where k = 1, 2, 3, 4, 5, 6, . . .. Then in year k + 1, the number of freshman
Chen CL
55
is given by x1(k + 1) = 10 percent of previous freshman class repeating freshman year +1000 new freshman = 0.1x1(k) + 1000 Because we know the number of freshman in the first year of our analysis (which is 500), we can solve this equation step to predict the number of freshman in the future. Let x2(k) be the number of sophomores in year k. Suppose that 15 percent of the freshman do not return and that 10 percent repeat freshman year. Thus 75 percent of the freshman class returns as sophomore. Suppose also 5 percent of the sophomores repeat the sophomore year and that 200 sophomores each year transfer from other schools. Then in year k + 1, the number of sophomores is given by x2(k + 1) = 0.75x1(k) + 0.05x2(k) + 200 To solve this equation we need to solve the ”freshman” equation at the same time, which is easy to do with MATLAB. Before we solve these equations, let us develop the rest of the model. Let x3(k) and x4(k) be the number of juniors and seniors in year k. Suppose that 5 percent of the sophomores and junior leave school and that 5 percent of the sophomores, juniors, and seniors repeat the grade. Thus 90 percent of the
Chen CL
56
sophomores and juniors return and advance in grade. The models for juniors and seniors are x3(k + 1) = 0.9x2(k) + 0.05x3(k) x4(k + 1) = 0.9x3(k) + 0.05x4(k) These four equations can be written in the following matrix form:
x (k + 1) 0.1 x (k + 1) 0.75 = x (k + 1) 0 1 2 3
x4(k + 1)
0
0 0.05 0.9 0
0 0
0 0
0.05 0 0.9 0.05
x (k) 1000 x (k) 200 + x (k) 0 1 2 3
x4(k)
0
In the following example, we will see how to use MATLAB to solve such equations.
Chen CL
57
Test Your Understanding T4.7-1 Suppose that 70 percent of the freshman, instead of 75 percent, return for the sophomore year. How does the previous equation change ?
Chen CL
58
Example: A College Enrollment Model (II) To study the effects of admissions and transfer policies, generalize the enrollment model in previous example to allow for varying admissions and transfers. Solution: Let a(k) be the number of new freshman admitted in the spring of year k for the following year k + 1 and let d(k) be the number of transfers into the following year’s sophomore class. Then the model becomes
x1(k + 1) = c11x1(k) + a(k) x2(k + 1) = c21x1(k) + c22x2(k) + d(k) x3(k + 1) = c32x2(k) + c33x3(k) x4(k + 1) = c43x3(k) + c44x4(k)
where we have written the coefficients c21, c22, and so on in symbolic, rather than numerical, form so that we can change their values if desired. This model can be represented graphically by a state transition diagram, like the one shown in the following Figure.
Chen CL
59
Such diagrams are widely used to represent time-dependent and probabilistic processes. The arrows indicate how the model’s calculations are updated for each values of x1(k), x2(k), x3(k), and x4(k); that is, by the vector x(k), which is called the state vector . The elements of the state vector are the state variables depend on both the previous values and the input a(k) and d(k). The four equations can be written in the following matrix form:
1) c xx (k + (k + 1) = c x (k + 1) 0 1
11
2
21
0 c22
0
c32 c33 0 0 c43 c44
3
x4(k + 1)
0 0
0 0
a(k) xx (k) (k) + d(k) x (k) 0 1 2 3
x4(k)
or more compactly as x(k + 1) = C x(k) + b(k)
0
Chen CL
where
60
x (k + 1) x (k + 1) x(k) = x (k + 1)
a(k) d(k) b(k) = 0
1 2 3
x4(k + 1)
and
cc C = 0
11 21
0
0 c22
0 0
0 0
c32 c33 0 0 0 c43 c44 Suppose that the initial total enrollment of 1480 consists of 500 freshman, 400 sophomores, 300 juniors, and 280 seniors. The college wants to study, over a 10-year period, the effects of increasing admissions by 100 each year and transfers by 50 each year until the total enrollment reaches 4000; then admissions and transfers will be held constant. Thus the admissions and transfers for the next 10 years are given by a(k) = 900 + 100k d(k) = 150 + 50k for k = 1, 2, 3, . . . until the college’s total enrollment reaches 4000; then admissions and transfers are held constant at the previous year’s levels. We cannot
Chen CL
61
determine when this event will occur without doing a simulation. The following Table gives the pseudocode for solving this problem.
Enter the coefficient matrix C and the initial enrollment vector x. Enter the initial admissions and transfers, a(1) and d(1). Set the first column of the enrollment matrix E equal to x. Loop over years 2 to 10. If the total enrollment is 4000, increase admissions by 100 and transfers
≤
by 50 each year. If the total enrollment is > 4000, hold admissions and transfers constant. Update the vector x, using x = Cx + b. Update the enrollment matrix E by adding another column composed of x. End of the loop over years 2 to 10. Plot the results. Because we know the length of the study ( 10 years), a for loop is a natural choice. We use an if statement to determine when to switch from the increasing admissions and transfer schedule to the constant schedule. A MATLAB script file
Chen CL
62
to predict the enrollment for the next 10 years appears in the following. The Figure shows the resulting plot. Note that after year 4 there are more sophomores than freshmen. The reason is that increasing transfer rate eventually overcomes the effect of the increasing admission rate. % Script file enroll1.m Computes college enrollment. % Model’s coefficients. C = [.1, 0, 0, 0; .75, .05, 0, 0; 0, .9, .05, 0; 0, 0, .9, .05]; % Initial enrollment vector. x = [500; 400; 300; 280]; % Initial admissions and transfers. a(1) = 1000; d(1) = 200; % E is the 4 x 10 enrollment matrix. E(:,1) = x; % Loop over years 2 to 10. for k = 2:10 % The following describes the admissions % and transfer policies if sum(x) <= 4000 % Increase admissions and transfers. a(k) = 900+100*k; d(k) = 150+50*k; else % Hold admissions and transfers constant. a(k) = a(k - 1); d(k) = d(k - 1); end % Update enrollment matrix.
Chen CL
63
b = [a(k);d(k);0;0]; x = C*x+b; E(:,k) = x; end % Plot the results. (Demo in the class) plot(E’), hold,... plot(E(1,:),’o’), plot(E(2,:),’+’), ... plot(E(3,:),’*’), plot(E(4,:),’x’), ... xlabel(’Year’), ylabel(’Number of Students’),... gtext(’Frosh’),gtext(’Soph’),gtext(’Jr’),gtext(’Sr’),... title(’Enrollments as a Function of Time’)
In actual practice this program would be run many times to analyze the effects of different admissions and transfer policies and to examine what happens if different values are used for the coefficients in the matrix C (indicating different dropout and repeat rates).
Chen CL
64
Test Your Understanding T4.7-1 In the previous program, line 16 and 17 compute the values of a(k) and d(k). These lines are repeated here: a(k) = 900 + 100 k
∗
d(k) = 150 + 50 k;
∗
Why does the program contain the line a(1) = 1000; d(1) = 200 ?
