程式流程控制 方煒 台大生機系
迴圈指令 MATLAB 提供兩種迴圈指令: for 迴圈 while 迴圈
for 迴圈 for 迴圈的使用語法如下: 運算式 end 另一種 for 迴圈的使用語法如下: for 變數 = 向量,
程式流程控制- for loop 下列 for 迴圈會產生一個長度為 6 的調和數列(Harmonic Sequence): 範例 1 : forLoop01.m x = zeros(1,6); % 變數 x 是一個 1×6 大小的零矩陣 for i = 1:6 x(i) = 1/i; end x % 顯示 x x = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 上例中,矩陣 x 最初是一個 1×6 大小的零矩陣,在 for 迴圈中,變數 i 的值依次是 1 到 6,因此矩陣 x 的第 i 個元素的值依次被設為 1/i。 可用分數形式來顯示此數列: >> format rat >> disp(x) 1 1/2 1/3 1/4 1/5 1/6
For k=m:s:n Statements End
練習 Summing a series with a for loop % script file: summing.m s=0; % set a variable to 0 so that 1/n^2 can be repeatedly added to it N=10000; % set the upper limit of the sum for n=1:N % add 1/n^2 to s each time, then put the answer back into s s=s+1/n^2; end fprintf(’ Sum = %g \n’,s) % print the answer % calculate the sum of the squares of the reciprocals of the % integers from 1 to 10,000 n=1:10000; sum(1./n.^2)
練習 Products with a for loop % script file: product.m P=1; % set the first term in the product N=20; % set the upper limit of the product for n=2:N P=P*n; end fprintf(’ N! = %g \n’,P) % print the answer factorial(20) gamma(21)
練習 Recursion relations with for loops %script file: recursion.m a(1)=1; % put the first element into the array N=19; % the first one is loaded, so let’s load 19 more for n=1:N a(n+1)=(2*n-1)/(2*n+1)*a(n); End format short disp(a) % display the resulting array of values format rat disp(a)
程式流程控制- for loop for 迴圈可以是多層或巢狀式(Nested)的,在下例中即產生一個 6×6 的Hilbert 矩陣 h,其中位於第 i 列、第 j 行的元素為 : 範例 2 : forLoop02.m h = zeros(6); % 變數 x 是一個 6×6 大小的零矩陣 for i = 1:6 for j = 1:6 h(i,j) = 1/(i+j-1); end format rat h h = 1 1/2 1/3 1/4 1/5 1/6 1/2 1/3 1/4 1/5 1/6 1/7 1/3 1/4 1/5 1/6 1/7 1/8 1/4 1/5 1/6 1/7 1/8 1/9 1/5 1/6 1/7 1/8 1/9 1/10 1/6 1/7 1/8 1/9 1/10 1/11
程式流程控制- for loop 在下例中,for 迴圈列出先前產生的 Hilbert 矩陣的每一直行的平方和: 範例 3 : forLoop03.m format short % 回到預設形式來顯式所有數值 for i = h disp(norm(i)^2); % 印出每一行的平方和 End 1.4914 0.5118 0.2774 0.1787 0.1262 0.0944 在上例中,由於 h 是一個矩陣,因此每一次 i 的值就是矩陣 h 的一直行的內容。 for i=1:6,disp(sum(h(i,:).^2)), end
程式流程控制- for loop 若要跳出 for 迴圈,可用 break 指令。例如,若要找出最小的 n 值,滿足 n! > 10100,可輸入如下: 範例 4 : break01.m for i = 1:1000 if prod(1:i) > 1e100 fprintf('%g! = %e > 1e100\n', i, prod(1:i)); break; % 跳出 for 迴圈 end 70! = 1.197857e+100 > 1e100
程式流程控制- for loop 在一個迴圈內若要直接跳至到此迴圈下一回合的執行,可使用 continue 指令。 範例 5 : continue01.m x = [1 -2 3 -4 5]; posTotal = 0; for i = 1:length(x) if x(i)<0, continue; end % 若 x(i) 小於零,跳到此迴圈的下一回合 posTotal=posTotal+x(i); end posTotal % 顯示 posTotal 的值 posTotal = 9
while 迴圈 while 迴圈使用語法如下: while 條件式 運算式; end
程式流程控制– while loop 先前產生調和數列的例子,亦可用 while 迴圈改寫如下: 範例 6 : while01.m x = zeros(1,6); i = 1; while i<=6 x(i) = 1/i; i = i+1; end x % 顯示 x x = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667
程式流程控制– while loop 若要用 while 指令找出最小的 n 值,使得 n! > 10100 ,可輸入如下: 範例 7 : while02.m n = 1; while prod(1:n) < 1e100 n = n+1; end fprintf('%g! = %e > 1e100\n', n, prod(1:n)); 70! = 1.197857e+100 > 1e100 與前述的 for 迴圈相同,在任何時刻若要跳出 while 迴圈,亦可使用 break 指令;若要跳到下一回合的 while 迴圈,也可以使用 continue 指令。 無論是 for 或 while 迴圈,均會降低 MATLAB 的執行速度,因此盡量使用向量化的運算(Vectorized Operations)而盡量少用迴圈。 break 指令若用在多重迴圈中,每次只跳出包含break指令的最內部迴圈。
練習 x = 5;k = 0; while x < 25 k = k + 1; y(k) = 3*x; x = 2*x-1; end
練習 total = 0;k = 0; while total < 1e+4 k = k + 1; total = 5*k^2 - 2*k + total; end disp(’The number of terms is:’) disp(k) disp(’The sum is:’) disp(total)
條件/分支指令 二種條件指令(Branching Command): if-then-else 條件指令 switch - case - otherwise 條件指令
if-then-else條件指令 if - then - else,其使用語法: if 條件式 運算式一; else 運算式二; end 當條件式成立時,執行運算式一,否則,就執行運算式二。 若不使用運算式二,則可省略 else 和運算式二。
程式流程控制- if-then-else 在數值運算的過程中,若變數值為 NaN(即 Not A Number )時,我們要立刻印出警告訊息,可輸入如下例: 範例8 : if01.m x = 0/0; if isnan(x) disp('Warning: NaN detected!'); end Warning: Divide by zero. … Warning: NaN detected! 在上例中,第一個警告訊息是 MATLAB 自動產生的,第二個警告訊息則是我們的程式碼產生的,其中 isnan(x) 可用於判斷 x 是否為 NaN,若是,則傳回 1(真),否則即傳回 0(偽)。
程式流程控制- if-then-else 可根據向量 y 的元素值為奇數或偶數,來顯示不同的訊息: 範例9 : if02.m for i = 1:length(y) if rem(y(i), 2)==0 fprintf('y(%g) = %g is even.\n', i, y(i)); else fprintf('y(%g) = %g is odd.\n', i, y(i)); end y(1) = 0 is even. y(2) = 3 is odd. y(3) = 4 is even. y(4) = 1 is odd. y(5) = 6 is even. 上述的 if - then - else 為雙向條件,亦即程式只會執行「運算式一」或「運算式 二」,不會有第三種可能。
練習 %script file: logic.m clear; a=1;b=3; % If the number a is positive set c to 1; if a is 0 or negative set c to 0 if a>0 c=1 else c=0 end % if either a or b is non-negative, add them to obtain c; % otherwise multiply a and b to obtain c if a>=0 | b>=0 % either non-negative c=a+b c=a*b % otherwise multiply them to obtain c
練習 x = [4,-9,25]; if x < 0 disp(’Some elements of x are negative.’) else y = sqrt(x) end 結果會是甚麼?
練習 x = [4,-9,25]; if x >= 0 y = sqrt(x) else disp(’Some elements of x are negative.’) end 結果會是甚麼?
elseif if logical expression 1 statement group 1 elseif logical expression 2 statement group 2 else statement group 3 end
練習 if x > 10 y = log(x) elseif x >= 0 y = sqrt(x) else y = exp(x) - 1 end
練習
程式流程控制- if-then-else 若要進行更多向的條件,只需一再重覆 elseif 即可。例如,欲判斷 y 向量之元素是屬於 3n、3n+1、 或 3n+2,可輸入如下: 範例10 : if03.m y = [3 4 5 9 2]; for i = 1:length(y) if rem(y(i),3)==0 fprintf('y(%g)=%g is 3n.\n', i, y(i)); elseif rem(y(i), 3)==1 fprintf('y(%g)=%g is 3n+1.\n', i , y(i)); else fprintf('y(%g)=%g is 3n+2.\n', i , y(i)); end y(1)=3 is 3n. y(2)=4 is 3n+1. y(3)=5 is 3n+2. y(4)=9 is 3n. y(5)=2 is 3n+2.
練習 陣列 A中不小於0的元素取平方根,其他元素則加上50,求新的陣列 A =[0 -1 4;9 -14 25;-34 49 64]; for m=1:size(A,1) for n = 1:size(A,2) if A(m,n)>=0 B(m,n)=sqrt(A(m,n)); else B(m,n)=A(m,n)+50; end A =[0 -1 4;9 -14 25;-34 49 64]; C=(A>=0); A(C)=sqrt(A(C); A(~C)=A(~C)+50;
Implied loop 內涵迴圈 for k=1:21 x=(k-1)*5 x=[0:5:100]; y(k)=cos(x); 等效於 end x=[0:5:100]; y=cos(X) 等效於 x=[-10,-5,0,5,10]; m=0; for k=1:length(x) If x(k)>0 m=m+1; y(m)=k; end x=[-10,-5,0,5,10]; y=find(x>0) 等效於
switch-case-otherwise 指令 MATLAB 在第五版開始支援 switch-case-otherwise 的多向條件指令,其使用語法如下: switch expression case value(1) statement(1) case value(2) statement(2) case value(n-1) statement(n-1) otherwise statement(n) end 在上述語法中,expression 為一數值或字串,當其值和 value(k) 相等時,MATLAB 即執行 statement(k) 並跳出 switch 指令。若 expression 不等於 value(k),k=1, 2, …, n-1,則 MATLAB 會執行 statement(n) 並跳出 switch 指令。
程式流程控制 – switch case 欲根據月份來判斷其季別,可輸入如下: 範例11 : switch01.m for month = 1:12 switch month case {3,4,5} season = 'Spring'; case {6,7,8} season = 'Summer'; case {9,10,11} season = 'Autumn'; case {12,1,2} season = 'Winter'; end fprintf('Month %d ===> %s.\n', month, season); Month 1 ===> Winter. . . . . Month 12 ===> Winter.
程式流程控制 - switch case 如果 expression 是字串,那麼若要在 case 之後比對多個字串,就必需使用字串的異值陣列(Cell Array of Strings): 範例12 : switch02.m month = {'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep'}; for i = 1:length(month) switch month{i} case {'Mar','Apr','May'} season = 'Spring'; case {'Jun','Jul','Aug'} season = 'Summer'; case {'Sep','Oct','Nov'} season = 'Autumn'; case {'Dec','Jan','Feb'} season = 'Winter'; end fprintf('%s is %s.\n', month{i}, season);
程式流程控制- switch case 上述範例output如下: Jan is Winter. Feb is Winter. Mar is Spring. Apr is Spring. May is Spring. Jun is Summer. Jul is Summer. Aug is Summer. Sep is Autumn. MATLAB 的 switch 指令和 C 語言的 switch 指令略有差別:在 C 語言的 switch 敘述內,每個 case 敘述需加上 break 以跳出該 switch 敘述,而在 MATLAB 則不必多此一舉。 一般而言,switch–case–otherwise 的執行效率優於 if–then–else 。
補充 - Using fzero function f=fz(x) % evaluate the function fz(x) whose % roots are being sought f=exp(-x)-x; %**************************************** % Here is the matlab code that uses fz.m to find % a zero of f(x)=0 near the guess x=.7 % Note that the @ sign is used to tell Matlab that % the name of an M-file is being passed into fzero %*************************************** x=fzero(@fz,0.7)
補充:inline, Example 1 This example creates a simple inline function to square a number. >>g = inline('t^2') g = Inline function: g(t) = t^2 >>g(3) ans = 9 You can convert the result to a string using the char function. >>char(g) t^2
補充:inline, Example 2 f = inline('3*sin(2*x.^2)') argnames(f) Inline function: f(x) = 3*sin(2*x.^2) argnames(f) ans = 'x' formula(f) 3*sin(2*x.^2)
補充:inline, Example 3 >>f = inline('sin(alpha*x)') Inline function: f(alpha,x) = sin(alpha*x) If inline does not return the desired function variables or if the function variables are in the wrong order, you can specify the desired variables explicitly with the inline argument list. >>g = inline('sin(alpha*x)','x','alpha') g = g(x,alpha) = sin(alpha*x)
補充 – inline , Example 4 func=inline(’exp(-x)-x’,’x’); x=0:.01:1; f=func(x); plot(x,f,’r-’)
切線法求根
切線法 (Secant method) 求根 % set chk, the error, to 1 so it won’t trigger % script file: secant01.m clear;close all; %************************************ % Define the function as an in line function func=inline(’exp(-x)-x’,’x’); % First plot the function x=0:.01:2; f=func(x); plot(x,f,’r-’,x,0*x,’b-’) % From the plot the solution is near x=.6 % Secant method to solve the exp(-x)-x = 0 % Use an initial guess of x1=0.6 x1=0.6; % find f(x1) f1=func(x1); % find a nearby second guess x2=0.99*x1; % set chk, the error, to 1 so it won’t trigger % the while before the loop starts chk=1; % start the loop while chk>1e-8 % find f(x2) f2=func(x2); % find the new x from the straight line approximation and print it xnew = x2 - f2*(x2-x1)/(f2-f1) % find chk the error by seeing how closely f(x)=0 is approximated chk=abs(f2); % load the old x2 and f2 into x1 and f1; then put the new x into x2 x1=x2;f1=f2;x2=xnew; end
切線法 (Secant method) 求根 script file: secant02.m func=inline('exp(x)-3*x','x'); x=0.6159 x=1.5121
function function showNearestPoint1(action) global h0 h1 h2 h3 h4 x1 y1 if nargin==0 action='b'; end x1=linspace(0,4*pi,30); y1=sin(x1); plot(x1,y1,'r-o') axis equal %
switch(action) case 'b' set(gcf,'windowbuttonmotionfcn','showNearestPoint1 a') case 'a' currPt=get(gca,'CurrentPoint'); x2=currPt(1,1); y2=currPt(1,2); Minx=sqrt((x2-x1).^2+(y2-y1).^2); n=find(Minx==min(Minx)); line(x1(n),y1(n),'marker','o') location=[num2str(x1(n)) ',' num2str(y1(n))]; %text(x1(n),y1(n),loction);%顯示在點上 text(x2+.2,y2,location);%顯示在游標旁 end