1/26/07

第四章 機構之自由度

在參考書第三章中討論到機構之自由度。如前所言,在二維運動中任何元件均應有三個自由度,即XY方向的移動與Z方向之轉動,亦即有三個自由度。在三維的空間則更多,亦即有XYZ各方向有移動而XYZ軸方向各有轉動,故有六個自由度。實際上,無論二維或三維,任何物件若運動太過自由,實際上也構不成機構,正如一個團體或國家放任每個人都能完全自由的話,就無法成為一個團體或國家一樣。

有些物件在運動過程中因而必須有約束力,這些約束力約會減少物件的自由度,或稱為拘束度。不過一個元件除了當固定桿外,其自由度不能為零,否則就無動產生運動。機構並不僅由一元件組成,通常由許多元件互相一運動結連繫而成。這一連繫就會使元件之運動受到約束,其自由度也會減少。機構之自由度減到為零時,表示這個機構變成靜物,不再有運動現象。這也是屬於一般靜力的範圍了。

一個機構,例如四連桿,由四個連桿元件及四個運動結組成。照前面所言,在二維的機構中,每一連桿應有三個自由度,四個連桿應有十二個自由度。但事實上可能這麼多嗎?讓我們來算它一算。由於成為機構之條件是要有一桿為固定,故要減去三個自由度;四個運動結對連桿而言目的在拘束連桿之運動,若每一個運動結若均為旋轉結,則每接一個結僅能傳遞一個自由度,亦即需減去二個拘束度,因此四個結必須減去八個拘束度。剛剛說,四連桿總共十二個自由度,減去三個再減去八個,僅剩一個自由度。所以四連桿僅能傳遞一個輸入變數。

上述之四連桿有一個自由度,為與一般元件的自由度判別,機構的自由度統稱為可動度,因為它是指向一個機構到底可以不可以動起來。若可動度為零,則屬靜態結構;若為負值則表示為靜不定結構,都是不可動的機構。

計算一個機構的可動度,通常可用古魯伯公式,下節將詳加介紹。


4.1古魯伯公式



古魯伯公式主要是根據下列計算公式演繹而來,可以計算一個機構之可動度M。可動度大於或等於一,才能成為可運動的機構,小於一時則變成結構。就計算的過程而言,一個機構之可動度應為其總自由度扣除總拘束度,其計算法則又稱為古魯伯(Grueber)公式,其規則如下:

  M = 3*(N-1)-(3J-Σf)
   = 3*(N-J-1)+Σf


式中,

  M = 可動度,即系統之自由度
  N = 連桿之總數
  J = 運動結之總數
   = 第i運動結之連結度

程式設計


程式gruebler是針對上式公式設計的一個函數。其主要功能在於先找出結的型式與數目。在本程式中係用利用jointype來達成。jointype是採用陣列的方式,依序安排節的位置與數量,其順序依序為:
  • 1 R型結
  • 2 滑動結
  • 3 複式結
  • 4 球結
  • 5 圓筒結
  • 6 平面結
  • 7 滾筒結
  • 8 凸輪
  • 9 螺旋結 
  • 10 球滑結
  • 11 點觸結

例如:jointype=[3 1],即表示有三個R型結,一個滑動結。另一個輸入參數為nlink,即為連桿之總數。jointype的安排上最大的缺點是若中間沒有該結時,必須留空,或在該位置填以零值。例如:一個R型結及三個球結,則jointype=[1 0 0 4]。

不過比較特別是,在jointype中若包括第四項以上時,這些結均為三維空間的結。故程式會以此作為判斷gruebler公式是計算二維或三維空間之桿結問題。因此,若一個三維空間的結,不含第四項球結以後的各項時,則至少在第四項必須以0表示,例如:jointype=[3 1 0 0],表示一組三維空間之連桿組,但具有三個R型結及一個滑動結。

在程式中,為求得對應於jointype中對應結之各結自由度,則利用code作為對應值,故每次必須選擇對應結之自由度。

程式內容


程式開始時,先決定jointype陣列之長度,除確定其維度外,同時作為計算各結自由度之基礎。利用njoint確定總結數,ff累計Σf之值,最後得到答案。


function [df]=gruebler(nlink,jointype)
%
% [df]=gruebler(nlink,jointype)
% nlink:no. of total links
% jointype:row matrix for number of joints for each type,
% the order of elements is:
% 1 R-joint 2 slider 3 compound joint(sliding & rolling)
% 4 ball 5 cylinder 6 planar 7 cylinder rolling
% 8 cam 9 helix 10 ball & 11 point contact
% Example: df=gruebler(4,[4])
% Author:D.S.Fon Bime,NTU. Date:Jan. 30, 2007
code=[1 1 2 3 2 3 1 2 1 3 5];
n=length(jointype);
dim=3;if n>3, dim=6;end;
ff=0;njoint=0;
for i=1:n,
njoint=njoint+jointype(i);
ff=ff+jointype(i)*code(i);
end;
df=dim*(nlink-njoint-1)+ff;
%-------------------

應用方法


其呼叫方式為:

[df]=gruebler(nlink,jointype)


其中之參數內容如下:

nlink:總桿數。
Jointype:行矩陣,代表屬於表3.1中所列各種型結之數目,其中各元素之內容為:
1 旋轉結(R) 2 滑動結(P) 3 複式結(滑動與轉動) (S) 4 球結(B) 5 圓筒結(C)
6 平面結(PL) 7 圓柱滾動(CR) 8 凸輪接觸(CP) 9 螺旋結(H)
10 球結(滾動與滑動)(BR) 11 點接觸結(SP)。

若以各元素代號表示即為

  Jointtype=[R P S B C PL CR CP H BR SP]

其中各代號之位置為該結之結數。中間沒有該結時,以零填之。例如:

     jointype=[2 0 1 0 1]

表示有兩個旋轉結,一個複式結,一個圓筒結。
例如:jointype=[2 0 1 0 1],表示有兩個旋轉結,一個複式結,一個圓筒結。

範例一


決定一組四連桿之可動度:N=4, R=4

>> df=gruebler(4,[4])
df = 1


範例二


三度空間連桿組具五連桿但有三個旋轉結(R型結)及三個球結。

>> df=gruebler(5,[3 0 0 3])
df = 0


範例三


但是同樣的四連桿若,若為三維,則完全沒有自由度,例如:

>> df=gruebler(4,[4 0 0 0])
df = -2

顯然,例一之平面四連桿僅是三維中之一特例。

gruebler函數之撰寫只是簡單的matlab程式之應用,但其輸入參數中使用的jointype比較複雜,你能想出更好的處理方法嗎?

4.2 機構之元件

機構之元件


機構之元件除前節所述之桿形(linkshape)外,其餘如滑塊、接地結、球結、圓筒結等都是常見的型態。本節就各種型狀再進一步論述。

接地結(ground)


接地結通常可以提供滑塊、圓筒結等在運動時與地接觸之表示法。通常接地結除繪一橫線表示地面外,為表示其不動的狀態,在其橫線下,用許多斜線表示接地。此接地結會將其相對應之值繪出,然後輸出其相對之座標。若僅想取其座標,不事先繪圖,則可以設dd值為負值。此處dd值代表由接地之原點往前與往後延伸之距離。

此函數ground執行時,必須有旋轉點(x0,y0)之資料,以及旋轉之角度theta。此外,尚有參數dd作為以旋轉點往兩側延伸之距離。其使用例將配合其他元件使用。


function [x,y]=ground(x0,y0,dd,theta)
% The ground line
% Inputs:
% x0,y0:the pivot point
% dd:distance on both sides
% theta:inclined angle, deg
% Example: ground(0,0,5,30)
d=abs(dd);d8=d/8;
xm=[-d/2:d/8:d/2];
xu=[xm;-d/16+xm;xm];xu=xu(:);
ym=zeros(size(xm));
yu=[ym;ym-d/8;ym];yu=yu(:);
x1=[-d;xu;d]+x0;y1=[0;yu;0]+y0;
th=-theta*pi/180;
x=x1*cos(th)-y1*sin(th);
y=x1*sin(th)+y1*cos(th);
if dd<0,
line(x,y,'linewidth',3);
axis equal;
end

此程式將配合後面之範例應用。

底座(anchor)


底座anchor函數為配合其他連桿結端點必須固定時,所需顯示之固定座。其外圍之圓弧半徑為rr,中心點(x0,y0)為桿端之連結位置,也是旋轉之迴轉中心。

此程式若直接接地,則會呼叫上述之ground函數,為使其僅傳回座標,參數rr值可設為正值,若要繪圖則設為負值。但rr之正負值並不影響圓弧之半徑。底座需要旋轉時,可以利用theta調整,其單位以度數表示。

function [x,y]=anchor(rr,x0,y0,theta)
% The anchor of a joint
% Inputs:r:outer radius of anchor
% x0,y0:center point
% theta:inclined angle, ang
% Example: anchor(5,0,0,30)
r=abs(rr);
th=0:pi/50:2*pi;r15=r*1.5;
x1=r*cos(th);y1=r*sin(th);
xx=x1/5;yy=y1/5;
th=-theta*pi/180;
[xm,ym]=ground(0,-r15,-r15,0);
x1=[0 xx x1 r15 xm' -r15 x1(50) 0]+x0;
y1=[0 yy y1 -r15 ym' -r15 y1(50) 0]+y0;
x=x1*cos(th)-y1*sin(th);
y=x1*sin(th)+y1*cos(th);
if rr<0,line(x,y);axis equal;end

此程式使用例與後面之範例合併。

滑塊(sldblk)


滑塊是方塊的物件,必須在一個平滑的面上移動。因此必須配合前述之接地結ground使用。但為使程式活用,本函數sldblk並沒有直接呼叫ground,因此必須由外部合併呼叫。其中LL,W分別為方塊之長與寬度,LL同時作為是否立即繪圖之指標,若使用負值,則會同時繪圖,否則僅傳回其外型之座標值。

方塊之中心(x0,y0)為參考點,可由連桿之端點在此連接。若方塊必須旋轉在某一方向,則方塊將以此點為旋轉中心,其旋轉角度為theta。


function [x,y]=sldblk(LL,W,x0,y0,theta)
% circular motion
% Inputs:LL,W:length and width of slider
% x0,y0:center point
% theta:inclined angle, degrees
% Example: sldblk(3,2,0,0,30)
L=abs(LL);th=0:pi/50:2*pi;
xx=W/5*cos(th);yy=W/5*sin(th);
th=-theta*pi/180;
%The coordinates of the blok
x1=[0 xx L/2 L/2 -L/2 -L/2 0 L/2 L/2]+x0;
y1=[0 yy 0 W/2 W/2 -W/2 -W/2 -W/2 0]+y0;
%Rotating the block by an angle th in radians
x=x1*cos(th)-y1*sin(th);
y=x1*sin(th)+y1*cos(th);
if LL<0,line(x,y);axis equal;end

執行例



% demo1.m
clf;
anchor(-0.6,0,0,0);%設底座,rr=-0.6表示要繪出圖形
linkxy([0 0],[2,2],-1);%繪第一桿
linkxy([2 2],[7,-1],-1);%繪第二桿
[x,y]=sldblk(-1.5,1,7,-1,30);%繪滑塊
ground(x(end-2),y(end-2),-2,30);%繪接地結



% demo2.m
clf;
anchor(-0.7,0,0,270);%繪底座一
linkxy([0 0],[-1,3],-1);%繪桿一
linkxy([-1 3],[-4 3],-1);%繪桿二
linkxy([-4 3],[-6 8],-1);%繪桿三
anchor(-0.7,-6,8,180);%繪底座二




滾動結(rolls)


滾動為一圓在另一平面上滾動。圓心則與另一桿相連,而圓周則在一平面上滾動。此平面可能為傾斜面,亦可能為平面,由前面的接地結函數ground決定。但是為符合滾動的原則,圓與地應有一起始點,以表示與接地結有一定的關係,為此可以在圓上畫一半徑,作為起點標示之用。

若地為水平狀態,圓之半徑與地垂直,成為最基本但也最簡單的狀態。若滾動接觸的地面傾斜某一角度時,此時必須連將圓與其一個半徑迴轉同一角度。程式之圓筒半徑為r,其中心點為(x0,y0),同時地面之傾斜角度為theta。

function [x,y]=rolls(r,x0,y0,theta)
% circular motion
% Inputs:r:radius of circle
% x0,y0:center of the circle
% theta:inclined angle, ang
% Example: rolls(5,0,0,30)
th=-pi/2:pi/50:1.5*pi;
xx=r*cos(th);yy=r*sin(th);
xi=xx/5;yi=yy/5;
th=theta*pi/180;
%x1,y1 are coordinates of circle and radius
x1=[0 xx xi 0];y1=[0 yy yi 0];
x=x1*cos(th)-y1*sin(th)+x0;
y=x1*sin(th)+y1*cos(th)+y0;
line(x,y);axis equal

執行例



% demo3.m
clf;
anchor(-0.7,-6,8,180);%上方底座
linkxy([-4 3],[-6 8],-1);%第一連桿
linkxy([0 0],[-4,3],-1);%第二連桿
[x,y]=rolls(2,0,0,30);%滾輪接地成30度
ground(x(2),y(2),-2,30);%接地結,與地成30度


4.3 連桿之種類

就與接頭之搭接情形而言,通常可分為雙接頭、參接頭及肆接頭或多接頭等。雙接頭在前面之連桿已經詳明,並且可以使用Matlab進行繪圖連接。多接頭之連結桿則是一種具有兩個以上之結點在同一個桿上。因此其形狀及外觀也特別與雙接頭桿不同。此節將設法介紹此類的桿型。

下面之程式為quatenary函數,可以輸入同一桿中之連結點位置(x0,y0),dd為各結之直徑,若為負值,表示需要繪圖。各點依序以在二維中最高點依反時針方向選取,因而可以得到正確的方位。

程式內容


此程式與前述之應用可以連結在一起,以獲得不同的機構表示法。不過,這個程式並不是最好的結果,讀者也許可以由此在研發更妥善的程式。此函數理論上可以處理許多點構成的多點桿,其結點之座標由(x0,y0)決定,因此它可以使用行陣表示,但必須一對一相互對應。dd為結孔之直徑。

每一結點之圓孔以dd為直徑,其座標利用linspacce(0,2*pi,60)以60點組成。基本上每一點均需由結之中心點開始,故這些座標必須相連接點出。其次決定外圍及內圍的線座標。


function [x,y]=quatenary(x0,y0,dd)
% Drawing for multiple links
% Inputs:
% x0,y0:coordinates of points,the points
% must be specified in counterclockwise.
% dd:dia. of circle
% Example
% [x,y]=quatenary([0.5 2 3 2],[1 3 2 0.5],-1);

nn=length(x0);d=abs(dd)/2;
th=linspace(0,2*pi,60);
xc=[0 d*cos(th) 0];yc=[0 d*sin(th) 0];
%collect the coordinates of all joint holes
x=[];y=[];
for i=1:nn,
x=[x x0(i)+xc];
y=[y y0(i)+yc];
end
x=[x x0(1)];y=[y y0(1)];
mx=mean(x0);my=mean(y0);r=1.3;
%find the outer coordinates
for i=1:nn,
x=[x (x0(i)-mx)*r+mx];
y=[y (y0(i)-my)*r+my];
end
x=[x (x0(1)-mx)*r+mx];
y=[y (y0(1)-my)*r+my];
if dd<0,patch(x,y,[.5 .5 .5]);axis equal;end


執行例



[x,y]=quatenary([0.5 2 3 2],[1 3 2 0.5],-0.1);




應用例



% demo4_1.m
clf;
quatenary([5 4 9],[6 2 5],-0.5);
anchor(-1,0,0,0);
linkxy([0 0],[2,5],-1);
linkxy([2 5],[5,6],-1);
linkxy([9 5],[8,3],-1);
linkxy([8 3],[10,0],-1);
anchor(-1,4,2,0);
anchor(-1,10,0,0);




Linkplate函數


上述繪製多結桿之函數,可以用另一函數改善,此函數稱為linkplate,其輸入及功能與quaternary函數相同,但外觀較佳。下面為此函數之內容:

程式內容


此程式之設計與前述之quaternary略有不同,它是利用複數型式進行向量之轉換,與前述之links函數之應用法相同。但其輸入項及內容與quaternary函數則完全相同。

function [x,y]=linkplate(x0,y0,dd)
% Drawing for multiple links
% Inputs:
% x0,y0:coordinates of points,the points
% must be specified in counterclockwise.
% dd:dia. of circle
% Example
% [x,y]=linkplate([0.5 2 3 2],[1 3 2 0.5],-1);
nn=length(x0);d=abs(dd)/2;
x0=x0(:);x2=[x0(2:end);x0(1)];
y0=y0(:);y2=[y0(2:end);y0(1)];
V=(x2-x0)+(y2-y0)*i;
D=abs(V);th=angle(V);
t=linspace(pi/2,2.5*pi,20);
Cout=max(d,0.2)*exp(j*t');Cin=Cout/2;
M=[];
for k=1:nn,
P=[Cin;Cout(1:10);D(k)+[Cout(11:20);Cin;Cout(20)]];
xx=real(P);yy=imag(P);tx=th(k);
x=xx*cos(tx)-yy*sin(tx)+x0(k);
y=xx*sin(tx)+yy*cos(tx)+y0(k);
if dd<0, line(x,y);end
end


執行例


前例亦可利用上面之函數執行,其指令及結果如下:

% demo4.m
clf;
linkplate([5 4 9],[6 2 5],-1);
anchor(-1,0,0,0);
linkxy([0 0],[2,5],-1);
linkxy([2 5],[5,6],-1);
linkxy([9 5],[8,3],-1);
linkxy([8 3],[10,0],-1);
anchor(-1,4,2,0);
anchor(-1,10,0,0);




另一執行例


此例包括滑塊、底座及多結桿等之應用,其指令檔及執行結果如下:

% demo5.m
clf;
anchor(-1,0,0,0);
linkxy([0 0],[1 2],-1);
linkxy([1 2],[6 3],-1);
linkplate([3 6 6],[6 3 7],-1);
linkxy([6 3],[8 0],-1);
anchor(-1,8,0,0);
linkxy([6 7],[8,9],-1);
anchor(-1,8,9,180);
linkxy([3 6],[-2,7],-1);
[x,y]=sldblk(-2,1.5,-2,7,-30);
ground(x(end-2),y(end-2),-2,-30);



4.4 四連桿之組成

參考書第四章中已經開始討論到四連桿機構的問題。你能列舉一些週遭事物中有那些使用四連桿機構的嗎?既然討論到機動問題,就必須看到它的運動狀態。而要有 運動更必須有參考座標,或參考桿,如此才能看出其相對運動間之關係。故在討論一個物件是否屬於機構,最大的考慮是必須有參考桿。

在四連桿中,其參考桿為第一桿,或稱地桿。這點在書中已經詳明,請同學自行參考。四連桿組顧名思義是由四支連桿相連而成,其中一桿為固定桿,另一桿(第二桿)為動桿。下面之程式是輸入已知之四桿長度及第二桿與水平面之夾角。實際上下述之程式仍然以很簡單的型式表示,而且必須藉助Matlab繪圖界面中之格點指示出第三桿與第四桿連結點之方位。才能進行繪圖。

第三桿與第四桿連結點通常有兩個位置,一為交叉型、一為閉合型。


function z=draw4links(L1,L2,L3,L4,theta)
%clf;
th=theta*pi/180;
r2x=L2*cos(th);r2y=L2*sin(th);
linkxy([0 0],[L1,0],-2);
linkxy([0 0],[r2x r2y],-2);
t=linspace(0,2*pi,100)';
r3x=L3*cos(t)+r2x;r3y=L3*sin(t)+r2y;
r4x=L4*cos(t)+L1;r4y=L4*sin(t);
line(r3x,r3y,'color','b','linestyle',':');
line(r4x,r4y,'color','r','linestyle',':');
line([r2x 118.5],[r2y 43.25],'linewidth',3);
line([L1 118.5],[0 43.25],'linewidth',3);
line([r2x 38.22],[r2y -56.61],'linewidth',3,'linestyle',':');
line([L1 38.22],[0 -56.61],'linewidth',3,'linestyle',':');





四連桿之位置分析將在另一章中敘述,本節僅能就簡單的方式繪出其相關位置。

4.5 葛拉索準則

葛拉索準則是針對四連桿之長度關係判斷其運轉情形。其第一準則是最長桿與最短桿之和小於其他兩桿和時,稱為葛拉索第一型,此時至少有一桿可完全迴轉。反之,三連桿之活動必屬搖桿,或稱葛拉索第二型。利用MATLAB程式可以進行判斷任一組四連桿是否屬於葛拉索型。

grashof函數第一式



程式4.1之呼叫函數為grashof():


function ans=grashof(ground_no,linkage)


其參數定義如下:
ground: 接地桿在linkage中之桿序。
Linkage: 列矩陣,四連桿長度,其次序不拘。

程式4.1
function ans=grashof(ground_no,linkage)
% Function to test the Grashof linkage
% Inputs:
% ground_no:the ground link number in the order
% linkage: row matrix for lengths of the 4 links
% in original assigned order.
% Example:ans=grashof(4,[4 4.2 2.6 2])
% Revised: March 4, 2006
ground=linkage(ground_no);
link=sort(linkage);% sorting the links
ig=find(linkage==link(1));
if link(1)+link(4)>link(3)+link(2),
ans='Non-Grashof Linkage';
elseif link(1)+link(4)==link(3)+link(2)
ans='Neutral Linkage';
elseif link(1)==ground,
ans='Double-Crank Linkage';
else
switch ig
case 1
im=3;
case 2
im=4;
case 3
im=1;
case 4
im=2;
end
if ground==linkage(im)
ans='Double-Rocker Linkage';
else
ans='Crank-Rocker Linkage';
end
end


範例4.1:四連桿組之長度如下:








L1 L2 L3 L4 接地桿號
a 2 4.5 7 8 1
b 3 5 4 4 2
c 3.5 4 1 5 2
d 4 5 3 7 2

(a) s=2, g=8, p=4.5,q=7
s+g=10<p+q=11.5

故應為葛拉索一型,至少有一桿為曲柄桿,其接地桿為最短桿s,故應為雙曲柄機構。用程式計算得知答案如下:

>> ans=grashof(1,[2 4.5 7 8])
ans =Double-Crank Linkage

(b) s=3, g=5, p=4,q=4
s+g=8 = p+q=8

因為相等,故屬於中立連桿組。用程式計算得知答案如下:

>> ans=grashof(2,[3 5 4 4])
ans =Neutral Linkage

(c) s=1, g=5, p=4,q=4
s+g=6<p+q=8

故應為葛拉索一型,至少有一桿為曲柄桿,其接地桿為第二桿,為最短桿s之側桿,故應為曲柄搖桿機構。用程式計算得知答案如下:

>> ans=grashof(2,[3.5 4 1 5])
ans =Crank-Rocker Linkage

(d) s=3, g=7, p=4,q=5
s+g=10>p+q=9

為非葛拉索型,無論任何桿接地均屬雙搖桿機構。用程式計算得知答案如下:

>> ans=grashof(2,[4 5 3 7])
ans =Non-Grashof Linkage


grashof2函數第二式


本程式係就原grashof函數修正而得。其輸入僅有四支連桿之長度,其第一桿為接地桿,其餘依序為第二、三、四桿。其輸出內容不立即印出,可由參數ans及代碼b作輸出,作為後繼程式呼叫之用途。


function [ans,b]=grashof2(linkage)
% revised from ans=grashof(groundno,linkage)
% Function to test types of Grashof linkages
% in which the first link is ground
% linkage: lengths of 4 links, in row matrix
% b: code(0-6),representing types of answers
% Example:[ans,b]=grashof2([4 4 5 6])
% Author: D S Fon, Bime, NTU, Oct. 25, 2007
link=sort(linkage);
ig=find(linkage==link(1));% find the shortest
if link(4)>=link(3)+link(2)+link(1)|length(linkage)~=4,
b=0;ans='Inputs fail to describe a 4-bar linkage!';
return;
end
if link(1)+link(4)>link(3)+link(2),
b=7;ans='Non-Grashof triple-rocker linkage';
return;
end
if link(1)+link(4)==link(3)+link(2)
if linkage(1)==linkage(2)&linkage(3)==linkage(4),
b=2;ans='Natural diamond double-crank linkage.';
elseif linkage(1)==linkage(3)&linkage(2)==linkage(4),
b=3;ans='Natural parrellel double-crank linkage.';
else
b=4;ans='General Natural linkage.';
end
return;
end
if link(1)==linkage(1),
b=1;ans='Double-Crank linkage with rotating coupler.';
elseif link(1)==linkage(3)
b=5;ans='Double-Rocker linkage with rotating coupler.';
else
b=6;ans='Crank-Rocker linkage with rotating crank.';
end
end


執行例


依據前述使用grashof函數之例子,重新執行如下:

>> [ans,b]=grashof2([2 4.5 7 8])
ans =
Double-Crank linkage with rotating coupler.
b =
1
>> [ans,b]=grashof([3 5 4 4])
??? Undefined command/function 'grashof'.

>> [ans,b]=grashof2([3 5 4 4])
ans =
General Natural linkage.
b =
4
>> [ans,b]=grashof2([3.5 4 1 5])
ans =
Double-Rocker linkage with rotating coupler.
b =
5
>> [ans,b]=grashof2([4 5 3 7])
ans =
Non-Grashof triple-rocker linkage
b =
7
>> [ans,b]=grashof2([4 4 5 6])
ans =
Non-Grashof triple-rocker linkage
b =
7

4.6 四連桿之解(1)

四連桿位置解


有關四連桿位置使用幾何代數解法請參照馮丁樹著之參考書四章三節。此例所考慮之四連桿是較簡單的型式,其前提是將第一桿置於水平位置,因此其角度為零,在公式計算上可以簡化。由於各桿之尺寸為已知,第二桿作為曲桿,其水平角度為已知,因此只要求得第四點的位置及第三與第四桿之水平角度即可。



程式內容


以MATLAB應用程式解四連桿之相關位置實際上有多種用法,可參考機動學之位置分析。本程式函數four_link1之呼叫法如下:

function [theta3,theta4,Cx,Cy]=four_link1(theta2,r,mode)

其中輸入參數:

  theta2: 桿2之輸入角度,可為矩陣資料,其單位為度數。
  r: 列矩陣,各桿之長度,如:[r1 r2 r3 r4]。
  mode:+1或-1,選擇連桿組上下之方位。

輸出參數:

theta3,theta4:桿3及桿4之輸出角度。
Cx,Cy:第四桿C結之座標,若為虛數值表示該位置不存在,其輸出角度不能採用,是否為虛數,則可使用isreal函數進行測試。



程式4.5
function [theta3,theta4,Cx,Cy]=four_link1(theta2,r,mode)
% P4.4 function [theta]=four_link1(theta2,r)
% Find the angles of link 3 and link, given theta2 and [r]
% Example:
% [theta3,theta4, Cx, Cy]=four_link1(80,[4 2 4.2 2.6],1)
% Designed by D.S. Fon, BIME, NTU, Date:February 8,2003.
%
rr=r.*r;
d2g=pi/180;
theta=theta2*d2g;
Bx=r(2)*cos(theta);By=r(2)*sin(theta);
m=(rr(4)-rr(3)-rr(1)+rr(2))./(Bx-r(1))/2;
mm=(m-r(1)).^2;
p=By./(Bx-r(1));pp=p.*p;
rootin=mm.*pp-(pp+1).*(mm-rr(4));
arg=sqrt(rootin);
Cy=((m-r(1)).*p+mode*arg)./(pp+1);
Cx=m-p.*Cy;
theta3=atan2(Cy-By,Cx-Bx)/d2g;
theta4=atan2(Cy,Cx-r(1))/d2g;


範例4.6


如圖4.1之連桿長度分別為70、50、110、65mm,桿2之角度為80度,求其對應之θ3與θ4

[解]:

設theta2=80,r=[70 50 110 65],mode=+1,代入four_link1執行。

>> [theta3,theta4,Cx,Cy]=four_link1(80,[70 50 110 65],1)

theta3 = -3.1199
theta4 = 41.7160
Cx = 118.5194
Cy = 43.2536


範例4.7


有一組四連桿,其尺寸分別為[4 2 4.2 2.6]cm,桿1接地,並與地平行。試(A)以45度為區隔,令桿2迴轉一周,求其對應之θ3與θ4及C結座標。(B)另以MATLAB撰寫一程式,繪出以桿2迴轉間隔為10度之所有四連桿之軌跡。

[解]:

設theta2=[0:60:360] ,總計應有6個角度。r=[4 2 4.2 2.6]cm ,mode=+1,代入four_link1執行。其結果如下:

>> [theta3,theta4,Cx,Cy]=four_link1([0:60:360],[4 2 4.2 2.6],1)
theta3 =
27.6604 8.1593 9.8818 21.5404 48.0950 68.1593 27.6604
theta4 =
48.5827 63.5647 109.3695 143.6226 147.5827 123.5647 48.5827
Cx =
5.7200 5.1575 3.1377 1.9067 1.8052 2.5625 5.7200
Cy =
1.9498 2.3281 2.4528 1.5421 1.3938 2.1665 1.9498


由上面所得的結果可知,桿2可以完全迴轉。而最前與最後一項之結果相同,因為0度與360度是同樣的角度。

4.7 四連桿之解(2)

第二種解(B)


MATLAB的指令中,有一簡單的指令可以繪圖,即line(x,y)指令。只要利用four_link1()計算對應於桿2迴轉角之四點座標(即點B及C,A與B固定),將座標值餵入line(x,y)指令中之座標參數(x,y)內,即可進行繪圖。注意每一個位置之四連桿座標繪完後,必須回到原點,才能繪出閉路的連桿型式,例如A-B-C-D-A的程序。故原點A在座標輸入參數項內需輸入兩次。

程式內容


利用MATLAB撰寫的程式four_link2內容如下:

程式範例4.7
function four_link2(theta2,r,mode)
% Program 4.5.2 to draw sets of 4-bar linkage
%Inputs:
% theata2:angles of link 2, in degrees
% r:row matrix for link lengths={r1 r2 r3 r4]
% mode:+1/-1 for toggling position
% Example:
% four_link2(0:5:360,[4 2 4.2 2.6],1)
% Author:D.S. Fon, BIME, NTU, Date:March 8,2006.
%
rr=r.*r;
d2g=pi/180;
theta=theta2'*d2g;
nn=length(theta);
Bx=r(2)*cos(theta);By=r(2)*sin(theta);
if abs(Bx-r(1))>0.01
m=(rr(4)-rr(3)-rr(1)+rr(2))./(Bx-r(1))/2;
mm=(m-r(1)).^2;
p=By./(Bx-r(1));pp=p.*p;
rootin=mm.*pp-(pp+1).*(mm-rr(4));
arg=sqrt(rootin);
Cy=((m-r(1)).*p+mode*arg)./(pp+1);
Cx=m-p.*Cy;
else
if abs(By)>0.01
Cy=(rr(4)-rr(3)-rr(1)+rr(2))/(2*By);
else
Cy=0;
end
Cx=r(1)+mode*sqrt(rr(4)-Cy^2);
end
%theta3=atan2(Cy-By,Cx-Bx)/d2g;
%theta4=atan2(Cy,Cx-r(1))/d2g;
Bx=r(2)*cos(theta);By=r(2)*sin(theta);
for i=1:nn,
if isreal(Cy(i))
line([0 Bx(i) Cx(i) r(1) 0]',[0 By(i) Cy(i) 0 0]',...
'marker','o');
end
end
axis equal;


執行結果:

>>four_link2(0:10:360,[4 2 4.2 2.6],1)

圖4.15 執行four_link2(0:10:360,[4 2 4.2 2.6],1)之結果

由圖4.15可知,左側桿可以迴轉360度,右側桿則僅能作往返運動,故為曲柄搖桿機構。若使用grashof2函數進行測試,其結果如下:

>> [ans,b]=grashof2([4 2 4.2 2.6])
ans =
Crank-Rocker linkage with rotating crank.
b =
6

其結果與實際繪圖的情況一致。故通常在繪圖之前,可以利用grashop2這個函數進行測試。立即能收到事半功倍之效。

設利用上述之draw4links2指令,設桿二之水平角度為80度,則可執行得到如下之結果:

>> draw4links2(4,2,4.2,2.6,80);

1/24/07

第三章 座標應用與分析

3.1複數之應用



在機動學之分析過程,複數亦可作為向量分析之方法,其形式類似極座標。向量表示法對照如下:




極座標垂直座標
R@∠θRcosθi+Rsinθj
ReRcosθi+Rsinθj


無論垂直座標或極座標,其絕對值均是以其某一個單位向量作為量度之基準。而在複數型式中,虛數用j表示(在matlab中i與j同義)。但此處j僅是一個代表虛數之操作元,表示其後的數值係屬於虛數軸上;在直角座標中亦相等於其Y軸。這個j的操作元在MATLAB中也有其特殊的意義。任何實數R與j相乘後,會在複數平面上反時針旋轉90度,即進入虛數的世界,或為Rj,但若將結果再與j相乘,則又會再旋轉90度,回到實數軸,但是已到負的方向,即Rj²= –R,這是因為j的平方等於-1的原故。故若該數再與j相乘,即Rj3= –Rj,則又回到負的虛數世界中了。





複數通常可用下列型式表示:

Z = x + jy
Z = x + iy

其中x稱為實數部,y稱為虛數部分,虛數部可用i或j代替。在MATLAB中有兩種表示法:

>> z=4+5j
z = 4.0000 + 5.0000i
>> z=4+5i
z = 4.0000 + 5.0000i

若其中之5的位置若為變數,則不能直接在數字後加i或j,必須在兩者間加一乘號,如下所述:

>> b=5
b = 5
>> z=4+bi
??? Undefined function or variable 'bi'.
>> z=4+b*i
z =  4.0000 + 5.0000i

因此上述型式中之實數與虛數部份均可用變數取代。若採用函數型式,複數亦可用complex(4,5)表示,如:

>> complex(4,5)
ans = 4.0000 + 5.0000i

其意義因此與前面之應用一樣。由於實數軸與虛數軸可以分為x 軸與y軸,故複數亦可作為表示平面座標點之方式,其與極座標之關係如下:

x = rcosθ
y = rsinθ

則複數可以應用下式表示:

Z=x+jy =rcosθ+jrsinθ


1-2.2尤拉(Euler)公式


尤拉公式為指數與複數間之轉換關係:

  e=cosθ+jsinθ


在直角座標系中,θ為任意一通過原點之向量r與x軸間之夾角。若將複數中之實數部份作為x軸上之分量,虛數部份作為y軸之分量,則兩分向量相互垂直,成為一組表示特定向量之方式。在尤拉式係使用自然對數為底,以jθ為指數,其關係與直角座標正好可以作轉換。
為使θ之角度旋轉符合極座標之不同方位,尤拉(Euler)公式可以作歸納成通式如下:

e±jθ=cosθ±jsinθ
------(1.4)


故如圖1.2中在四個象限中相互垂直(或間隔90度)之單位向量之關係為:

ej(θ+π/2)=cos(θ+π/2)+jsin(θ+π/2)
=-sinθ+jcosθ
=j(cosθ+jsinθ)=je-------(1.5)

ej(θ+π)=cos(θ+π)+jsin(θ+π)
=-sinθ-jcosθ
=-(cosθ+jsinθ)=-e
=j²e---------------------(1.6)

ej(θ+π3/2)=cos(θ+π3/2)+jsin(θ+π3/2)
=sinθ-jcosθ
=-j(cosθ+jsinθ)=-je
=j3e -------------------(1.7)

例:設θ=30度,試利用尤拉公式證明其轉換值相同,並計算每相隔90度之對應複數值。

>> theta=30*pi/180 %轉換為弧度
theta =
0.5236
>> A=exp(j*theta) %求得左邊之對應複數值
A =
0.86603 + 0.5i
>> a=cos(theta)+j*sin(theta) %求得右邊對應值,與前述相同
a =
0.86603 + 0.5i
>> angle=theta:pi/2:2*pi %設定區間為90度之角度值
angle =
0.5236 2.0944 3.6652 5.236

>> M=exp(j*angle') %對應之複數值
M =
0.86603 + 0.5i
-0.5 + 0.86603i
-0.86603 - 0.5i
0.5 - 0.86603i

3.2 複數之運算

一般運算


設某一複數Z=x + yj。若 Z=0,表示其實數部份與虛數部份之係數應為0,即 x=0; y=0。而共軛複數(Congugate) 亦應成立。同理若有兩個複數Z1=x1+y1j與Z2=x2+y2j若相等,即表示x1=x2;y1=y2。因此,一個複數之等式,應可分為兩個等式方程式。此外,在複數的國度裡,下面的運算恆真:

Z1 ± Z2=(x1 ± x2)+j(y1 ± y2)
Z1*Z2=(x1*x2-y1*y2)+j(x1*y2+x2*y1) -------(1.8)
Z1/Z2=(x1+jy1)/(x2+jy2)
=(x1+jy1)(x2-jy2)/(x2+jy2)(x2-jy2)
=[x1*x2+y1*y2+j(x2*y1-x1*y2)]/[x2²+y2²]--(1.9)

設Z之共軛複數z=x-jy,則下列式亦恆真:

Z*z=x²+y²
Z+z=2x =real(Z)
Z-z=2jy=imag(Z)
conj(Z1 ± Z2)=z1 ± z2
conj(Z1*Z2)=z1*z2
conj(Z1/Z2)=z1/z2 ------- (1.10)

在MATLAB中有一指令conj可以設定共軛複數,例如:

>> Z=10+6j
Z =
10 + 6i
>> z=conj(Z) %共軛複數
z =
10 - 6i
>> Z2=3+5j
Z2 =
3 + 5i
>> Z*Z2 %複數相乘
ans =
0 + 68i
>> Z/Z2 %複數相除
ans =
1.7647 - 0.94118i
>> Z*z %複數與其共軛複數相乘
ans =
136
>> Z+z %複數與其共軛複數相加
ans =
20
>> Z-z %複數與其共軛複數相減
ans =
0 + 12i
>> conj(Z+Z2) %兩複數和之共軛複數
ans =
13 - 11i
>> conj(Z)+conj(Z2) %與其共軛複數和相同
ans =
13 - 11i
>> conj(Z*Z2) %兩複數積之共軛複數
ans =
0 - 68i
>> conj(Z)*conj(Z2) %等於兩共軛複數之積
ans =
0 - 68i
>> conj(Z/Z2) %兩複數除數之共軛複數
ans =
1.7647 + 0.94118i
>> conj(Z)/conj(Z2) %等於兩共軛複數相除
ans =
1.7647 + 0.94118i

範例


設兩複數分別為 z1=2+j(3)1/2、z2=5ej(pi/4) ,試利用MATLAB求 z1+z2、z1*z2、z1/z2等項之值。

[解]

>> z1=2+sqrt(3)*j
z1 = 2 + 1.7321i
>> z2=5*exp(pi/4*j)
z2 = 3.5355 + 3.5355i
>> a=z1+z2
a = 5.5355 + 5.2676i
>> b=z1*z2
b = 0.94734 + 13.195i
>> c=z1/z2
c = 0.52779 - 0.037894i

複數之微分


上面各式之關係可以在運算過程中善加利用。將複數式進行對θ微分,會具有如下之型式:

d/dθ[e] = je-----(1.11)

證明上式可依尤拉轉換公式,將其轉為正弦與餘弦函數:

  d/dθ[e] =
=d/dθ(cosθ+jsinθ)
=-sinθ+jcosθ=j(cosθ+jsinθ)
=je

設有一向量B=10e,其中θ=60度,經對角度θ微分後,其值應為:

dB/dθ=10je

利用matlab指令可依序計算及比較如下:

>> theta=60*pi/180 %角度轉為弧度
theta =
1.0472
>> B=10*exp(j*theta) %向量B 
B =
5 + 8.6603i
>> dB=10*j*exp(j*theta) %第一微分後之向量dB
dB =
-8.6603 + 5i

>> BB=[real(B) imag(B) 0] %將向量B轉為三維分向量,最後一項為零
BB =
5 8.6603 0
>> DB=[real(dB) imag(dB) 0] %將向量dB轉為三維分向量
DB =
-8.6603 5 0
>> cross(DB,BB) %兩向量應相垂直,其叉積應在Z軸上
ans =
0 0 -100
>> dot(DB,BB) %兩向量相垂直,其點積應為零
ans =
0

因此,一般之向量可以應用複數指數表示,然後依其微分之方式求得對應之梯度值。以圖1.1為例。rp 之位置向量可為下列三個向量之和:

rp = r1+r2+r3
= r1e1+r2e2+r3e3
=r1(cosθ1+jsinθ1)+r2(cosθ2+jsinθ2)+r3(cosθ2+jsinθ2)----(1.12)

故實際的分析仍然回到三角函數之垂直向量的領域。其微分的過程亦循其在直角座標之分向量分別進行。

3.3 MATLAB處理指令

在MATLAB中,複數的表示法如下:

 Z = R.*exp(j*theta) ----------------(1.13)

此與R 值之表示相同,而j與i均代表虛數部份,兩者同義。注意此處theta為其水平夾角,以弧度表示;R後有一點,表示為元素單項對乘,不是矩陣乘法,因此R值亦可用矩陣型式輸入;而j*theta則無點,代表其為虛數處理。不過上式若R值僅單一項,則可不用加點。

範例一


一向量值為10,角度分別為45度、60度及90度,求其對應複數表示形式,且當其陸續乘以j之後,其對應複數值之變化。
[解]

>> a=10*exp(j*pi./[4 3 2])
a = 7.0711 + 7.0711i 5.0000 + 8.6603i 0.0000 +10.0000i
>> b=a*j
b = -7.0711 + 7.0711i -8.6603 + 5.0000i -10.0000 + 0.0000i
>> c=b*j
c = -7.0711 - 7.0711i -5.0000 - 8.6603i -0.0000 -10.0000i
>> d=c*j
d = 7.0711 - 7.0711i 8.6603 - 5.0000i 10.0000 - 0.0000i
>> e=d*j
e = 7.0711 + 7.0711i 5.0000 + 8.6603i 0.0000 +10.0000i

注意其間符號及其代表象限之變化,至e值時,正好環繞一週,故應與a值相同。

若以同一向量旋轉四個象限,即分別在0、π/2、π、3π/2、2π之角度位置,則其相對複數值應為:

>> format long
>> a=10*exp(j*pi.*[0 1/2 1 3/2 2]')
a =
10.00000000000000
0.00000000000000 +10.00000000000000i
-10.00000000000000 + 0.00000000000000i
-0.00000000000000 -10.00000000000000i
10.00000000000000 - 0.00000000000000i

上述結果之最後一項值實際上與第一項值相同。

此外Z值若為複數,其向量夾角及絕對值可用下面指令求得:

THETA = angle(Z)
R = abs(Z) -----------------------(1.14)

設 Z =10+5i,用MATLAB執行如下:

>>Z =10+5i
Z= 10.0000 + 5.0000i
>> theta=angle(Z)
theta = 0.4636
>> r=abs(Z)
r = 11.1803

theta得到的角度是弧度,若使用度數時可用下面指令轉換:

>> angle=theta*180/pi
angle = 26.5651


1-2.6 複數向量之合成


利用複數之實數與虛數兩種不同性質的領域可以在二維向量上作某種程度的搭配,在解題的技巧上可以簡化觀念。平面物體的運動可以分成兩個互相垂直的座標分量,因而可在二度空間中定域向量之變化。下面的程式函數comb_vectors可以用來計算多個複數向量之合向量與夾角,其輸入參數可為各對應向量之複數值,也可以是各之絕對值與其對應之水平夾角,輸出值為總合向量之值與角度。此處之夾角均以水平夾角之度數表示。

程式內容


此程式容許兩種輸入型式。其一是Z為絕對向量值陣列,theta為其對應之水平夾角度數;其二是僅有一項輸入Z,其型式為複數向量陣列。執行後,其輸出值為合向量及水平角,並在圖中顯示。

function [Zp,thetaZ]=comb_vectors(Z,theta)
% Find the final position vector rp consisting
% of complex vectors Z
% using complex vector method
% Zs are considered as complex forms when theta is not used
%
% Example: [rp, thetap]=comb_vectors([3+4j 3+4j 3+4j])
d2g=pi/180;
if nargin==2,
Z=Z.*exp(i*theta*d2g);
end
format short;
Z=[0+0i;Z(:)];
Zp=sum(Z);
Rz=abs(Zp);%開方根求其絕對值
thetaZ=angle(Zp)/d2g;%求得相角,並轉回一般角度
plot(cumsum(Z));hold on;%繪出各組成向量
plot([0+0i;Zp],'r-');%繪出結果紅線向量
xlabel('X-axis');
ylabel('Y-axis');
title(['Resultant:' num2str(Rz) ' Angle(deg.):' num2str(thetaZ)])
axis equal;grid on

執行例:輸入三個向量,長度分別為[3 4 5]公分,其對應水平角度分別為[35 50 60]度,試利用上式comb_vectors函數求得其合向量:


>> [rp, thetap]=comb_vectors([3 4 5],[35 50 60])
rp =
7.5286 + 9.1150i
thetap =
50.4448


若輸入的三個複數向量,分別為[2+1j 2+3j 1+5j]公分,可利用同一函數求得此三個向量之合向量及水平角度:

>> [rp, thetap]=comb_vectors([2+1j 2+3j 1+5j])
rp =
5.0000 + 9.0000i
thetap =
60.9454



輸入向量及角度[3 4 5],[35 50 60]



輸入為複數向量[2+1j 2+3j 1+5j]


由上述之應用可知:若平面向量能轉為複數表示,其運算可分別依實數與虛數部份個別計算。

3.4 極座標轉換指令

除利用複數型式計算座標向量,在MATLAB中也有一個角度座標之轉換指令,可以將極座標系之向量轉化為直角座標:

[X,Y] = pol2cart(THETA,RHO)
[X,Y,Z] = pol2cart(THETA,RHO,Z) ----------------(3.15)

在極座標中,通常有徑向位移與夾角,類似複數指數在尤拉系的轉換。極座標之指令中,利用pol2cart可將角度(THETA,弧度)及位移量(RHO )轉換為直角座標之[X,Y]值。實際上其函數即含有X=RHO*cos(THETA),Y=RHO*sin(THETA) 之運算。




當然有了[x,y]值時亦可利用下式或cart2pol指令還原其徑向位移與角度:

theta = atan2(y,x)
rho = sqrt(x.^2 + y.^2)
[theta, rho]=cart2pol(x,y)---------------(3.16)

注意函數中有用到atan2(y,x)這個函數,在應用上此與atan(x)單參數有區別,atan2(y,x)具有雙參數,可以得到四個象限的位置,依X與Y項之值決定,其範圍為 -pi <= atan2(y,x) <= pi。 例如:X=-10,Y=5,則

>> theta1=atan(5/(-10))*180/pi
theta1 = -26.5651
>> theta1=atan2(5,-10)*180/pi
theta1 = 153.4349
>> cart2pol(-10,5)*180/pi
ans =
153.4349

顯然第一個答案僅給正負值的角度,而第二答案明顯告訴我們答案為第二象限。值得比較的是cart2pol反而可以得到考慮象限的答案,只是其輸入順序與atan2不同,必須注意。

座標轉換實例

範例2:利用上述三個位置向量r1、 r2、 r3分別為5(10°)、8(135°)、9(12°),求其直角座標?

function [rp,thetap]=Xresult(theta,rho)
% Find the resultants from connecting
% vectors in forms of rho(theta)
% theta in degrees
% Example: [rp, thetap]=Xresult([30 40 60],[3 4 5])
d2g=pi/180; %轉換為弧度
format short
[x,y]=pol2cart(theta*d2g,rho); %將角座標換為垂直座標
rpvec=[0 0;cumsum([x',y'])]; %行矩陣之加總
[theta,rp]=cart2pol(rpvec(:,1),rpvec(:,2));
thetap=theta/d2g;
hold on;
for i=2:length(rpvec)
plot([0 rpvec(i,1)],[0 rpvec(i,2)],'r-')
text(rpvec(i,1),rpvec(i,2),...
[' ' num2str(rp(i),3) '/' num2str(thetap(i),3)])
end
plot(rpvec(:,1),rpvec(:,2),'b*-');
xlabel('X-axis')
ylabel('Y-axis')
title('Find Reultants of Vectors')
axis equal

執行例:向量長[3 4 5]公分,對應水平角度[30 40 60]°,求其合向量及累計向量。

>> [rp, thetap]=Xresult([30 40 60],[3 4 5])
rp =
0
3.0000
6.9739
11.7134
thetap =
0
30.0000
35.7161
45.8268




範例3:試用極座標轉換指令作一繪圓程式。
說明:Target程式係利用pol2cart指令將圓之半徑與角度轉換為直角座標,然後進行繪製粗細不等之同心圓。圓之外圍半徑及點數可以輸入設定,沒有輸入參數時,其預設值r=1,n=100。


function h=target0(r,n)
%Draw circles using pol2cart command
% r,n:radius and number of point in a circle.
% Example: h=target
if nargin<1, r=1; n=100; end
theta=linspace(0,360,n)/180*pi;
rr=linspace(0,r,10);
for k=2:10,
[x,y]=pol2cart(theta,rr(k));
h=line(x,y,'linewidth',11-k,'color',[1 k k]/10);
end
axis equal


執行target指令,沒有輸入參數,其結果如圖1.4。

>> target0
ans =
166.0016




圖3.4 同心圓的繪製

上述之同心圓問題若改成繪製同心圓時,具有動畫的氣氛,使同心圓能依序繪製,則上述程式內容可以修改如下:

function h=target1(r,n)
%Draw circles using pol2cart command
% r,n:radius and number of point in a circle.
% Example: h=target1
if nargin<1, r=1; n=100; end
theta=linspace(0,360,n)/180*pi;
rr=linspace(0,r,10);
h=line(0,0);
while 1
for k=2:10,
[x,y]=pol2cart(theta,rr(k));
set(h,'xdata',x,'ydata',y,'linewidth',11-k,...
'erasemode','xor','color',[1 k k]/10);
pause(1);
drawnow;
end
end
axis equal


同學們不妨存檔然後執行,看其變化如何。若要改變其顏色,則應如何修改才能獲得所要的結果?

3.5 速度分析

有了位置向量,其對時間之導數即為速度,故理論上應可由位置向量之表示法進行演繹延伸,得到相關之速度。速度的因變數為時間,即為位移之時間變化率。在動力學方面,機構之轉動則常以轉動之角度為因變數,因此常必須轉換。下面所述即為導數之轉換過程。設某一任意連桿之位置向量為r,其複數表示如下:

 r = re ------(3.17)

3-5.1 速度公式


物體之位移對時間t微分的結果設為V,則速度V為:

V=dr/dt=d/dt(re)=r'e+jrθ'e
=r'e+jrωe
=(r'+jrω)e
=(r'+jrω)(cos(θ)+sin(θ))
=r'e+rωej(θ+π/2)------------------------------------------(3.18)

若連桿組由三桿組成,則就三連桿之位置向量進行微分,可以得到如圖1.1中之最終速度Vp之變化,即:

Vp=rp'=r1'+r2'+r3'
k=1-3(rk'+jrkωk)ejθk
k=1-3(rk'+jrkω)(cos(θk)+sin(θk))
k=1-3(rk'cos(θk)-rkωksin(θk)+j(rk'sin(θk)+rkωkcos(θk))
----------------------------------------(3.19)

上式中若桿長固定或徑向無速度變化時, r'項為零,故可簡化如下:

Vp=rp'=r1'+r2'+r3'
k=1-3jrkωk)ejθk
k=1-3jrkω(cos(θk)+jsin(θk))
k=1-3rkωk(-sin(θk)+jcos(θk))-----------------------------(3.20)

若變換為直角座標時,則對應如下:

Vp=rp'=r1'+r2'+r3'
k=1-3(rk'cos(θk)-rkωksin(θk)+j(rk'sin(θk)+rkωkcos(θk))
------------------------------(3.21)

程式內容


程式velo2.m可以輸入位置向量及對應之角速度,由此計算其合位移及速度,其中包括絕對值及對應角度。在此程式中,假設徑向量無變化,故成為多桿連結的狀況。各桿之長度與角度採一對一的型式。角速度可以單一值亦可配合各桿之對應值。若不輸入則以一單位代替。

輸出部份有sv與thetap,前者中有兩項,一項為位移絕對值,另一項為速度之絕對值。thetap則為對應之水平夾角,以度數表示,亦分為兩項,前者屬合位移之夾角,後者則為合速度之夾角,兩者依原理應相互垂直,或相差九十度。


function [sv,thetap]=velo2(theta,rho,omega)
% Find the velocity Vp for the
% connecting vectors with theta,rho and w
% theta: angle in degrees for each vector
% rho: radial displacement, or link length
% omega: angular velocity of each vector
% Example: [sv,thetap]=velo2([20 40 50],[3 4 5],2)
theta=theta(:);rho=rho(:);
n=length(rho);
if nargin==2,
omega=ones(size(rho));
elseif length(omega)==1
omega=ones(size(rho))*omega;
end
omega=omega(:);
d2g=pi/180; %轉換為弧度
itheta=i*theta*d2g;
ss=rho.*exp(itheta);
vv=i*rho.*exp(itheta).*omega; %將極座標換為複數座標
sp=sum(ss);vp=sum(vv); %行矩陣之加總
sv=[abs(sp) abs(vp)];%位移與速度之絕對值
thetap=[angle(sp) angle(vp)]/d2g;%位移與速度之相角,並轉回一般角度

設有支連桿以旋轉結相連,其位置向量分別為[3 4 5]公分,角度為[20 40 50]度,角速度均為5rad/s,利用velo2函數可以計算其合位移及速度向量,其對應量存於sv(1)及sv(2),其對應向量水平角分別為thetap(1),thetap(2):

>> [sv,thetap]=velo2([20 40 50],[3 4 5],5)
sv =
11.744 58.721
thetap =
39.23 129.23


下面之程式是將各桿之位置圖,速度向量圖及合速度圖繪出。利用函數vectordraw可以得到比較的圖案。其程式結構如下:


function [ss,vv]=vectordraw(r,dg,w)
% Example: [ss,vv]=vectordraw([10 20 30],[30 40 50],[1 2 3])
d2g=pi/180;
r=r(:);dg=dg(:);w=w(:);
ith=exp(i*dg*d2g);
s=r.*ith;v=i*s.*w;s=[0;s];
x0=real(s(1));y0=imag(s(1));
for k=1:length(r)
x1=real(s(k+1))+x0;y1=imag(s(k+1))+y0;
line([x0 x1],[y0 y1],'marker','o','linewidth',2);
line([x1 x1+real(v(k))],[y1 y1+imag(v(k))],'marker','>>','color','r')
x0=x1;y0=y1;
end
ss=[real(sum(s)) imag(sum(s))];
vv=[real(sum(v)) imag(sum(v))];
line([0 ss(1)],[0 ss(2)],'linestyle',':','linewidth',2)
line([ss(1) ss(1)+vv(1)],[ss(2) ss(2)+vv(2)],...
'marker','>','color','g','linewidth',2)
axis equal;grid on

設連桿長為[10 20 30]度,對應水平角度為[30 50 90]cm,其對應角速度為[1 2 3]rad/s,則可依上述函數指令執行如下:

>> vectordraw([10 20 30],[30 50 90],[1 2 3])
ans =
21.5160 50.3209



藍色代表連桿位置,紅色為各桿之速度,虛線及綠色為合速度之方向及大小

3.6 加速度分析

若將速度量再對時間t微分則為加速度A,即:

A=dV/dt=d/dt(r'e+jrθ'e)
=[(r"+jr'θ')+(jr'θ'+jrθ"-rθ'²)]e
=[(r"-rω²)+j(2r'ω+rα)]e
=[(r"-rω²)+j(2r'ω+rα)](cosθ+jsinθ)
=[(r"-rω²)cosθ-(2r'ω+rα)sinθ)]-...
j[(r"-rω²)sinθ+(2r'ω+rα)cosθ]-----------(3.21)

3.6.1加速度公式


上式若桿長為固定時,即 r'=0,則上式可簡化為:

A=dV/dt=d/dt(jrθ'e)
=[(-rω²)+jrα]e
=[-rω²+jrα](cosθ+jsinθ)
=[-rω²cosθ+rαsinθ)]+j[-rω²sinθ+rαcosθ]---(3.22)

若就圖2之轉換象限之觀念而言,上式加速度亦可改變如下:

A=-rω²e+jrαe
= rω²(j²e)+rα(je)
= rω²ej(θ+π)+rα(jej(θ+π/2))---(3.23)

加速度向量具有二兩項分項量,一為指向負的徑向(相差180度),即為向心加速度;一為指向切線方向,成90度方向,與速度相同。

3.6.2 MATLAB加速度函數


任何連桿運動時,均會產生加速度,其值之計算則可利用Matlab函數acc2.m。此函數之輸入項包括各桿之長度(rho)、水平角度(theta)、對應之角速度(omega)及角加速度(alfa)。輸出則為絕對值與合成桿之夾角。每項內有三元素分別說明位置、速度及加速度之向量。

acc2.m程式中,係採用複數表示法計算位移、速度及加速度量。輸入項中若僅輸入前兩項,則角速度設為一單位,若僅三項則角加速度設為零。連桿數則可以無限。


function [sv,thetap]=acc2(rho,theta,omega,alfa)
% Find the acceleration Ap for the
% connecting vectors with theta,rho,omega & alfa
% theta: angle in degrees for each vector
% rho: radial displacement
% omega: angular velocity of each vector
% alfa: angular accerlation of each vector
% Example:
%[sv,th]=acc2([5 8 9],[10 135 12],[10 20 30],[5 10 20])
theta=theta(:);rho=rho(:);
n=length(rho);
%handling the default values
if nargin<4,
alfa=zeros(size(rho));
if nargin==2,
omega=ones(size(rho));
end
end
if length(omega)==1,omega=ones(size(rho))*omega;end
if length(alfa)==1,alfa=ones(size(rho))*alfa;end
%
omega=omega(:);alfa=alfa(:);
d2g=pi/180; %轉換為弧度
xtheta=exp(i*theta*d2g);
ss=rho.*xtheta;
vv=i*rho.*xtheta.*omega; %將極座標換為複數座標
aa=(-rho.*omega.^2+i*rho.*alfa).*xtheta; %將極座標換為複數座標
sp=sum(ss);vp=sum(vv);ap=sum(aa); %行矩陣之加總
sv=[abs(sp) abs(vp) abs(ap)];%位移與速度之絕對值
thetap=[angle(sp) angle(vp) angle(ap)]/d2g;%位移與速度之相角,並轉回一般角度


執行例1:僅有一桿,其長度為5公分,角速度為10rad/s,角加速度為20rad/s²,則執行結果如下:

>> [sv,th]=acc2(5,10,20)
sv =
5 100 2000
th =
10 100 -170

結果其位移量仍為5cm,絕對速度為100cm/s,絕對加速度為2000cm/s²。其對應角度方面,位移仍維持原有之10度,絕對速度之位角為100度,與前者相差正好90度;絕對加度之相差為-170度,與原位移相差正好180度。

執行例2:若有三桿,其位置向量分別為[5 8 9]cm,對應位移角度為[10 135 12]度;角速度分別為10, 20, 30 rad/s,無角加速度,求三桿合併後之總運動狀態。

>> [sv,th]=acc2([5 8 9],[10 135 12],[10 20 30])
sv =
11.646 267.86 7357
th =
46.133 131.63 -146.75


執行例3:同前例,但各桿之角加速度均為3rad/s²,求三桿合併速度及加速度。

>> [sv,th]=acc2([5 8 9],[10 135 12],[10 20 30],3)
sv =
11.646 267.86 7364.9
th =
46.133 131.63 -147.02


3.6.3 應用實例



多桿連結之分析


多桿連結是一個開放型機構,故其端點之位置完全依各桿之長度與角度。實際上其分析仍然依位置、速度及加速度等三種項量加以整合。但在速度與加速度的場合必須有各桿之角度、角速度及角加速度之參數,才獲得整合的答案。當然有些情況可以讓其中一個桿之輸入參數變化,其餘保留為定數,因而可以瞭解某特定桿變動時對整個機構之影響。多桿連結分析在工程機械上應用較多,如怪手、機器人等均屬此類。

下面的MATLAB程式dyad(),其呼叫語法如下:

function [vec, dyadata] = dyad(rho,theta,td,tdd)

讀者不妨自行執行後,觀察其結果。左邊之輸出參數中,vec為位置,速度及加速度之絕對值及其對應水平角度;dyadata則為各桿之對應值。


function [vec,dyadata] = dyad(rho,theta,td,tdd)
%
% function [vec, th, dyadata] = dyad(rho,theta,td,tdd)
% Analyzes a dyad linkage composing a crank and a dyad.
% Inputs: rho:length of links
% theta:incling angles, deg.
% td:angular velocity, rad/s
% tdd:angular acceleration, rad/s^2
% Outputs: vec:absoute length of links
% dyadata:original data, in complex forms
% Example:[vec,dyadata] = dyad([5 10],[30 50],[2 4])

theta=theta(:);rho=rho(:);
n=length(rho);
if nargin<4,
  tdd=zeros(size(rho));
  if nargin==2,
  td=ones(size(rho));
  end
end
if length(td)==1,td=ones(size(rho))*td;end
if length(tdd)==1,tdd=ones(size(rho))*tdd;end
td=td(:);tdd=tdd(:);
d2g=pi/180;
tt=exp(i*theta*d2g);
pp=rho.*tt;vv=i*td.*pp;aa=-pp.*td.^2+i*pp.*tdd;
dyadata=[pp vv aa];
vec=[abs(sum(dyadata));angle(sum(dyadata))/d2g];

執行範例:

>> [vec, dyadata] = dyad([5 10],[30 50],[2 4])
vec =

14.7976 49.5152 178.9247
43.3637 136.0392 -132.1910


dyadata =

1.0e+002 *

0.0433 + 0.0250i -0.0500 + 0.0866i -0.1732 - 0.1000i
0.0643 + 0.0766i -0.3064 + 0.2571i -1.0285 - 1.2257i


如何利用上面dyad函數程式可以用來計算相關連桿之速度與加速度,然後進行繪圖。此時亦須應用到第二章介紹之 linkshape函數,以繪出一般連桿之外形。下面的程式dyad_draw先呼叫dyad取得相關資料,再利用 linkshape繪製。linkshape函數附於此程式之後,可以隨時呼叫。程式內容如下:

function dyad_draw(rho,theta,td,tdd)
% Call function dyad(rho,theta,td,tdd) to
% analyze a dyad linkage composing cranks and a dyad with drawing.
% Inputs: rho:length of links
% theta:incling angles, deg.
% td:angular velocity, rad/s
% tdd:angular acceleration, rad/s^2
% Example:
% dyad_draw([15 10 5],[0 60 90],[.5 .8 .5],0)
clf;
[vec,data] = dyad(rho,theta,td,tdd);
x=[0;cumsum(real(data(:,1)))];y=[0;cumsum(imag(data(:,1)))];
for i=1:length(x)-1
linkshape([x(i) y(i)],[x(i+1) y(i+1)],1);
end
for k=1:length(rho)
x0=x(k+1);y0=y(k+1);
vx=x0+real(data(k,2));vy=y0+imag(data(k,2));
ax=x0+real(data(k,3));ay=y0+imag(data(k,3));
line([x0 vx],[y0 vy],'marker','s','linewidth',2);
line([x0 ax],[y0 ay],'marker','s','color','r','linewidth',3)
end
sdata=sum(data);
ss=[real(sdata(1)) imag(sdata(1))];
vv=[real(sdata(2)) imag(sdata(2))];
aa=[real(sdata(3)) imag(sdata(3))];
line([0 ss(1)],[0 ss(2)],'linestyle',':','linewidth',2)
line([ss(1) ss(1)+vv(1)],[ss(2) ss(2)+vv(2)],...
'marker','d','color','g','linewidth',3)
line([ss(1) ss(1)+aa(1)],[ss(2) ss(2)+aa(2)],...
'marker','d','color','m','linewidth',3)
axis equal;grid on

function linkshape(A,B,dd)
% Draw a link
% Inputs:
% A,B:Initial & final coordinates of link
% d:thickness of link(negative for line link)
% Example: linkshape([0 0],[5,5])
if nargin==2,dd=1;end;
d=abs(dd);
AB=(B(1)+j*B(2))-(A(1)+j*A(2));
D=abs(AB);th=angle(AB);
t=linspace(pi/2,2.5*pi,20);
Cout=max(d/2,0.2)*exp(j*t');Cin=Cout/2;
if dd>0,
P=[0;Cin;Cout(1:10);D+Cout(11:20);D+Cin;D+Cout(20);Cout(1)];
else
P=[Cin;0;D;D+Cin];
end
xx=real(P);yy=imag(P);
x=xx*cos(th)-yy*sin(th)+A(1);
y=xx*sin(th)+yy*cos(th)+A(2);
line(x,y)
axis equal

執行例:
設有三桿相連結成端桿之型式。各桿之長度分別為[15 10 5]cm,[0 60 90]度,[.5 .8 .5]rad/s^2。設各桿之角加速度為0,利用上面之函數執行如下:


dyad_draw([15 10 5],[0 60 90],[.5 .8 .5],0)

其結果如下:


其中之紅色為各桿之加速度,藍色為速度,綠色為合速度

1/22/07

第二章 動畫之製作

機動學中,物件之動作與位置變化是重要的考慮。如何利用Matlab之動畫功能,以瞭解物件之動作過程,也是重要的一環。動畫之中必須配合繪製物件,其中包括圓形、矩形及線形等,這些繪圖方法及語法,已在另外章節加以說明,讀者可以參考。目前在此處比較重要的是繪出物件之示意圖,所用之工具也較為簡單。


線為圖之基礎,適當地組合線群可以形成複雜的圖形。畫線指令line可以完成這項任務。繪線必須有起始點及終止點,但若有連續的點,採用矩陣表示時,則會 自動連線,繪出整個圖形。故若使用line指令,應可以利用點群之座標,一筆畫繪出需要的圖形。多點連線時,參數可以矩陣表示,但每行代表一線,亦亦採用 多行,形成方形矩陣,下面為其一例:

2.1.1 簡單之線圖


實例1.1:若x=[1 3 5]; y=x²

>> x=[1 3 5],y=x.*x
x = 1 3 5
y = 1 9 25

>> line(x',y')
>> line(x,y)
>> grid on


圖2.2 一維向量之繪圖

2.1.2一筆畫線


利用line之功能可繪製一筆畫之圖形。這種畫法可以利用方格紙先將初稿繪於其上,然後點出其對應座標,即可利用上述line之畫線功能繪製圖案。這種圖案未來可作為程式中應用之基本元件,構成更複雜的機構。以三角形為例,可以將座標分為對應之兩行如下:

triangle=[0 0;1 0;0 2;-1 0;0 0];

其中頭尾均為共同點,以保證整個圖為封閉區線。

實例1.5: 繪出三角形



>> triangle=[0 0;1 0;0 2;-1 0;0 0]
triangle =
0 0
1 0
0 2
-1 0
0 0
>> line(triangle(:,1),triangle(:,2))
>> axis equal
>> line(triangle(:,1),triangle(:,2))


圖2.3 繪出一三角形

實例1.6:星形



>> star5=[4 0;-3.2361 2.3511;1.2361 -3.8042;1.2361 3.8042; -3.2361 -2.3511;4 0];
>> line(star5(:,1),star5(:,2))
>> axis equal %調整座標之比例使成1:1



圖2.4 繪出一星形

實例1.7:房屋



>> house=[25 94;8 56;8 16;44 5;44 52;44 5;60 16;72 12;92 28;...
92 62;92 28;72 12;72 60;72 12;60 16;60 50;60 25;...
48 16;56 21;56 40;48 36;48 16;44 16;16 40;36 36;...
36 52;16 56;16 40;8 50;2 46;24 96;48 44;62 52;52 72;...
60 88;76 52;96 64;80 96;60 88;68 92;60 108;24 96];
>> line(house(:,1),house(:,2))
>> line(house(:,1),house(:,2))
>> axis equal



圖2.5 房屋

實例1.8:斗笠



>> banhat=[0 0;5 0;2 1;0 3;-2 1;-5 0;5 0;4 0;2 1;1.5 0.99;1 0.96;...
0 0.94;-1 0.96;-1.5 0.99;-2 1;-4 0;-3 0;-1.5 0.99;0 3;...
-0.5 0.95;-1 0;-2 0;-1 0.96;-2 0;2 0;1 0.96;2 0;1 0;...
0.5 0.95; 0 3;1.5 0.99;3 0;4 0;2 1;1.5 0.99;0 3;0 0];
>> line(banhat(:,1),banhat(:,2))
>> line(banhat(:,1),banhat(:,2))
>> axis off
>> axis equal


圖2.6斗笠


利用上述之方式,應該可以創出許多不同的作品,唯一需要的是個人的創意。

2.2 畫線的指令參數

在matlab中,完整的繪線指令語法如下:

line(X,Y)
line(X,Y,Z)
line(X,Y,Z,'PropertyName',PropertyValue,...)
line('PropertyName',PropertyValue,...) low-level-PN/PV pairs only
h = line(...)

line 指令可以產生線性之物件,也可以設定線之顏色(Color)、寬度(LineWidth)、線型(LineStyle)及線標(Marker)等。線函數具有兩項特性:其一是自動依序排色及線型,在未有設定之前,函數會使用預設值。此函數中,輸入資料前三項為繪線之座標,其大小應該一致。簡化之指令如 line(X,Y,Z),其顯示之顏色則會依ColorOrder,而線型會依 LineStyleOrder兩參數值循環使用。這個功能與一般之plot指令相同,只是其繪線時不像後者會呼叫newplot函數更新圖面。

在最低層次的應用上,這個函數亦可以利用"屬性名稱/屬性值"成對的型式輸入座標資料。其方式如 line('XData',x,'YData',y,'ZData',z)指令,輸入之位置與時間可以不受限制,但三者必須同時出現,且順序也要一致。這 些座標值必須為向量陣列,不能為一般之矩陣型式。

MATLAB 係利用預設之線顏色在目前之軸上繪圖。其顏色設定可依 colordef函數為之。每執行line(X,Y)一次,即會將定義於 X與 Y陣列之連線繪製於目前的座標軸內。若X與Y同為同大小之矩陣,則會依每行一線的方式繪製。line(X,Y,Z)指令則會在三維座標內繪線,若需另加參 數值,則可利用"屬性名/屬性值"之輸入方式,另外加入,如:line(X,Y,Z, 'PropertyName', PropertyValue,...) 是。其他未加入之部份則仍使用預設值。同二維的情況,繪線之資料點亦可利用下列之格式輸入:

line('XData',x,'YData',y,'ZData',z,'PropertyName',PropertyValue,.. .)

但記得此處之x,y,z值必須為陣列,不得為一般矩陣。每次繪線之握把可使用 h = line(...)設定。後面之項目例如顏色與線寬之設定等:

line(X,Y,Z,'Color','r','LineWidth',4)
line('XData',x,'YData',y,'ZData',z,'Color','r','LineWidth',4)

第一種型式是前三或二項(二維)之位置是預留為座標點,其內容可為矩陣,若為多行,則每行繪一線,且X,Y,Z必須同大小。第二種雖然意義相同,但由於使用配對關係,可以置於參數中任何意位置,此時x,y,z必須為行陣列。

2.2.1操作例



t=linspace(0,5*pi,200);
m0=zeros(size(t));
line('xdata',t,'ydata',sin(t),'zdata',m0,'color',...
'b','linestyle','-','linewidth',3)
line('xdata',t,'ydata',m0,'zdata',t.^3,'color','k',...
'linestyle',':','linewidth',4)
line('xdata',m0,'ydata',sin(t),'zdata',t.^3,'color',...
'g','linestyle','-.','linewidth',5)
line('xdata',t,'ydata',sin(t),'zdata',t.^3,'color',...
'r','linewidth',10)
grid on;
view(53,90)




同樣,下列程式也可以產生相同的結果:

t=linspace(0,5*pi,200);
m0=zeros(size(t));
line(t,sin(t),m0,'color','b','linestyle','-','linewidth',3);
line(t,m0,t.^3,'color','k','linestyle',':','linewidth',4);
line(m0,sin(t),t.^3,'color','g','linestyle','-.','linewidth',5)
line(t,sin(t),t.^3,'color', 'r','linewidth',10)
grid on;
view(14,42)


下面程式也可以產生同樣的結果,只是每條線均設有握把(handle)可以呼叫,因此可以在事後增加或改變參數:

t=linspace(0,5*pi,200);
m0=zeros(size(t));
line1=line(t,sin(t),m0);
line2=line(t,m0,t.^3);
line3=line(m0,sin(t),t.^3);
line4=line(t,sin(t),t.^3);
set(line1,'color','b','linestyle','-','linewidth',3);
set(line2,'color','k','linestyle',':','linewidth',4);
set(line3,'color','g','linestyle','-.','linewidth',5)
set(line4,'color', 'r','linewidth',10)
grid on;
view(14,42)


握把名稱也可以使用陣列,相同的特質可以一併處理:

t=linspace(0,5*pi,200);
m0=zeros(size(t));
h(1)=line(t,sin(t),m0);
h(2)=line(t,m0,t.^3);
h(3)=line(m0,sin(t),t.^3);
h(4)=line(t,sin(t),t.^3);
for i=1:4, set(h(i),'linewidth',i+2);end
set(h(1),'color','b','linestyle','-');
set(h(2),'color','k','linestyle',':');
set(h(3),'color','g','linestyle','-.')
set(h(4),'color', 'r');
grid on;
view(14,42)


也可利用三維的繪圖指令plot3,但記得要使用hold on才不會洗掉前一個圖。不過仍然會有些不同,你能發覺不同之處嗎?


t=linspace(0,5*pi,200);
m0=zeros(size(t));
h(1)=plot3(t,sin(t),m0,'bo-');hold on;
h(2)=plot3(t,m0,t.^3,'k:');
h(3)=plot3(m0,sin(t),t.^3,'g-.');
h(4)=plot3(t,sin(t),t.^3,'r-');
for i=1:4, set(h(i),'linewidth',i+2);end
grid on;
view(14,42)

Line函數之操作例


這是另一種動畫型態,可參考youtube.com網站。

2.3 連桿外型之繪製

連桿之位置及其外觀設計如下圖。連桿兩端之中心點分別為AB,由此兩點界定連桿之尺寸。因此其兩點間之距離D為重要的因素。由於繪製連桿之外圍,必須利用line指令採用描線之方式。這裡所提出之描線方法與前節介紹之內容相同,是屬最簡單之方式。首先設連桿之厚度為d,兩端之外圓與內圓之半徑分別為d/2與d/4。其連線方式由左邊之A點開始,先往上走1,再繞左端之內圓,回到原端點時再往上走3。至外圓時往左走半圓4,經過5,與右邊之半圓6 ,至7往下再走內圓8,其次回到9往上,再往右沿10回到3,完成一線畫。



由於連桿係以平放的位置開始,實際上必須計算BA線段之角度,然後進行迴轉到所需要之位置。此時就必須利用第一章之座標轉移與迴轉,配合運用。

程式內容


為使任何連桿均能具有厚度,以及所需要之長度,程式linkshape就是一種函數寫法,並且主要能繪出連桿之外形。此函數之輸入有點A、B與d三項,分別為連桿之起點及終點座標及桿厚度。程式採用複數表示法,因此可以迅速計算其D值與傾斜角度。

在Matlab中,複數表示法有其特殊的性質,可以作為表示兩獨立座標之值。在複數表示法中,實數i的對應值可代表x軸之向量;而虛數j部份之值可代表y軸之向量,因此在二維的座標中,正好可以利用一個複數值代表一個座標點。其優點是Matlab中有指令如real()、imag()、abs()、angle()等函數可以分別處理實數、虛數、絕對值及角度等值,應用上有其方便之處。所以在linkshape之函數中,可以見證其應用的情形。

連桿AB間,其對應關係式可以轉換如下:

x = x'cos(th)-y'sin(th) +A(1)
y = y'sin(th)+y'cos(th) +A(2)

公式裡,th為其水平角度,A(1)與A(2)分別為起點A之座標值。程式內容如下:

function linkshape(A,B,dd)
% Draw a link
% Inputs:
% A,B:Initial & final coordinates of link
% d:thickness of link(negative for line link)
% Example: linkshape([0 0],[5,5])
if nargin==2,dd=1;end;
d=abs(dd);
AB=(B(1)+j*B(2))-(A(1)+j*A(2));
D=abs(AB);th=angle(AB);
t=linspace(pi/2,2.5*pi,20);
Cout=max(d/2,0.2)*exp(j*t');Cin=Cout/2;
if dd>0,
P=[0;Cin;Cout(1:10);D+Cout(11:20);D+Cin;D+Cout(20);Cout(1)];
else
P=[Cin;0;D;D+Cin];
end
xx=real(P);yy=imag(P);
x=xx*cos(th)-yy*sin(th)+A(1);
y=xx*sin(th)+yy*cos(th)+A(2);
line(x,y)
axis equal

執行例



>> linkshape([15 0],[0,0])
>> linkshape([15 0],[5,6])
>> linkshape([0 0],[5,6])




利用同樣的程式,試執行下列指令,看看會有什麼樣的結果。

>> linkshape([15 0],[0,0],-1)
>> linkshape([15 0],[5,6],-1)
>> linkshape([0 0],[5,6],-1)


設有一四連桿,其ABCD四點之座標分別為A(0,0);B(3,4);C(13,4);D(10,0),其單位為cm,若AD為固定桿,AB為第二桿BC與CD分別為第三與第四桿,各桿厚度分別為3,2,1,3cm。試繪出其相關位置。

若AB為主動迴轉桿,則其每間隔60度間之對應位置會如何?這可寫成一個執行檔如下:

linkshape([0 0],[10 0],3)
linkshape([0 0],[3 4],2)
linkshape([3 4],[13 4],1)
linkshape([13 4],[10 0],3)


請用help patch指令,查詢patch指令的用法,如何更改linkshape函數中之畫圖之顏色?

2.4 動畫之製作例

為配合機動學的講授,使用MATLAB的繪圖功能,可以產生物體移動的效果。前文曾討論有關座標的轉移方式,可以先將某一已設定的圖案作平移或旋轉,或者兩者合併,以達到不同的需求與目的。這些座標轉換指令為:

function xprime=trans4(x,delx,mod)

其中x代表全部的座標矩陣,mod則為轉換的模式,當mod=1時表示移動,其移動量依delx決定,此為一列矩陣,對應於座標數。Mod=5時為轉動,若為平面,此時delx代表轉動的角度,以度數表示。若為立體座標,則其內容應為三元素之列矩陣。

當一個動作包含有移動與轉動時,上述指令必須依功能需要呼叫兩次。最好的處理方法是先行處理純轉動座標變換,其後再作移動轉換,如此才不會因旋轉中心改變而產生意外的結果。

連桿之轉動


為取得特定連桿之外圍座標,可以將上節之linkshape函數稍作修改,使其輸出為一組可用之(x,y)座標(註:本程式未使用trans4函數之呼叫,但直接擷取其內容):

function [x,y]=links(A,B,d)
% Find coordinates of a link
% Inputs:
% A,B:Initial & final coordinates of link
% d:thickness of link
% Output:[x,y]:Coordinates of the link
% Example: links([0 0],[5,5])
if nargin==2,d=1;end;
AB=(B(1)+j*B(2))-(A(1)+j*A(2));
D=abs(AB);th=angle(AB);
t=linspace(pi/2,2.5*pi,20);
Cout=(d/2)*exp(j*t');Cin=Cout/2;
P=[0;Cin;Cout(1:10);D+Cout(11:20);D+Cin;D+Cout(20);Cout(1)];
xx=real(P);yy=imag(P);
x=xx*cos(th)-yy*sin(th)+A(1);
y=xx*sin(th)+yy*cos(th)+A(2);


運動之過程



利用MATLAB繪製動畫有兩種方法,其一是先儲存一定數格的圖格,然後再按順序播放,成為連續的電影;其二是在螢幕上連續繪圖及擦拭,以產生連續的動作變化。圖面內容若相當複雜,且有色彩搭配時,以採用第一種方式為宜,如此可以產生一個接近實體之連續影像。這些影像可以事先輸入或繪製,播放時只要將其一幕幕放映即可。這種方式可能會產生閃耀不穩,且所佔的記憶體也較大。

第二種技巧是利用電腦的快速繪圖功能,將影像繪製完成後,再擦去,使其產生視覺暫留的印象。利用MATLAB所提供的模式可以選擇不同的呈現方式。這種方式速度快,但正確性較差。在動力或機動學方面,可以示範組件的動作,故也相當理想。

擦拭模式


MATLAB 指令中,簡單的圖可以使用擦拭重畫法以產生動畫。產生動畫之原理是將某一物體之座標改變之後,擦去前面所繪的圖像,重新依新座標繪製一次,如此連續不斷,產生一系列的動作畫面。這個動作須令'erasemode'的參數值等於'xor'才能完成。典型的畫法是畫出該物體後,然後在迴圈中,不斷地些微改變其位置座標,然後畫出新的畫面。其座標分別以x,y,z表示。在指令中亦可直接設定'erasemode'之值達到程式控制的目的。'erasemode'的參數有三種選擇:

  • none - MATLAB 會在繪出新圖案前不將舊圖清除。
  • background - MATLAB 會用背景顏色將舊圖重新繪製一次。故此種模式除將舊圖清除外,其底下原有的繪製物如格線等亦被清除。
  • xor – 此模式僅擦拭物體部份,最常用於動畫。

這三種模式是MATLAB繪圖指令中之執行速度最快者,但仍然視圖像的細節及運算的時間而定。若變化太過快速,可利用pause指令來延遲經歷之時間。在機動學中之動作中,本節將以此功能配合各項分析結果產生動畫,使同學瞭解其中之變化過程。

程式內容


本程式呼叫links函數,以得到單桿之外形之座標點,然後依與時間相關之迴轉角度delt順序旋轉。擦拭模式可以在程式中依下述的方法在繪製時就先設定:

h1=line(x,y,'erasemode','xor','color','r','linewidth',2);

指令中,h1為其握把,以此認定所需改變的線圖。注意此時之'erasemode' 參數值設定為 'xor',表示使用擦拭模式,亦即在下一次繪製新圖時,會將上次繪製之舊圖塗去,以觀察連續變動之影像。在每次繪製前必須將新的座標指向該握把,可利用此指令為之:

 set(h1,'xdata',x2,'ydata',y2);

採用set指令與直接下 line有點不同, line指令下達後,繪圖動作立即進行,但set指令僅先設定所要的值在緩衝器上,並不立即執行。因此你可以執行set好幾次,等到滿意之後才一次執行。執行繪圖的動作須配合drawnow之下達,才能一次繪出。下面的程式中,係先執行drawnow,表示先執行第一次準備之資料,但這只是程式的應用技巧,到底drawnow事先下或後下,依使用情況而定,一般並不影響結果。故實際上應用上也可置於執行set指令之後。

function link_rot1(r,d,v)
%link_rot1.m will call links.m
%Rotating a link
%Example: link_rot1(10,3,50)
%Author:DSFon, BIME,NTU. Date:Jan. 24, 2007
clf;
if nargin==0,r=10;d=1;v=50;end
[x,y]=links([0 0],[r,0],d);
delt=pi/120; %delaying time for display
axis([-r r -r r]*1.2);
axis image;axis off;
t=0:2*pi/50:2*pi;
h0=line(r*cos(t),r*sin(t),'color','b','linestyle',':');
h1=line(x,y,'erasemode','xor','color','r','linewidth',2);
x1=get(h1,'xdata');y1=get(h1,'ydata');
title('Press Ctl-C to stop');
theta1=0;s1=delt;
while 1,
drawnow;
pause(1/v);
theta1=theta1+s1;
x2=x1*cos(theta1)-y1*sin(theta1);
y2=x1*sin(theta1)+y1*cos(theta1);
set(h1,'xdata',x2,'ydata',y2);
if theta1>2*pi,theta1=theta1-2*pi;end;
end



執行結果




按此觀賞動畫:

2.5 利用linkshape繪製連桿組

前節提出之linkshape函數可以依輸入之兩端點繪出具有厚度之單一連桿,本節則試利用這個函數繪出不定數目之連桿組,所得之連桿必須每桿頭尾相連,下面為輸入端點座標後連成之程式,稱為link_plot,這個函數必須呼叫linkshape,此函數已在本章第2.3節中討論過:

桿座標輸入



function link_plot(x,y,w)
%link_plot.m will call links.m
%draw connecting linkshape function
%Input:(x,y)are serial coordinates
% w:width of links
%Example: link_plot([0 3 5 1],[0 3 2 1],3)
%Author:DSFon, BIME,NTU. Date:Jan. 25, 2007
clf;
if nargin==2,w=1;end
for i=1:length(x)-1
linkshape([x(i) y(i)],[x(i+1) y(i+1)],w);
end


這個程式雖短,但只要有一組座標輸入即可繪出其連桿組外型。例如:連續四桿之座標分別為(0,0), (3,3), (5,2), (5,5)等四點,其輸入型式為:

>> link_plot([0 3 5 4],[0 3 2 5])


執行結果如下圖:



桿長度與角度輸入

前面程式link_plot之輸入為各桿端點之座標。但有些時候,一組連桿常僅知道桿長及桿與前桿所夾之角度。這種表示法普見於機器手臂之討論。本節提出之link_angle程式,就是實現這個想法。此時之輸入參數為各桿之長度(lens),對應水平角度(angle)及起點(start)與桿寬(w)等。

function link_angle(lens,angle,start,w)
%link_plot.m will call links.m
%draw connecting links
%Input:lens:lengths of links
% angle:link angles w/t the former link, degs
% w:width of links
% start:starting point(x0,y0)
%Example: link_angle([5 3 4 6],[30 60 90 145])
%Author:DSFon, BIME,NTU. Date:Jan. 25, 2007
if nargin<4,w=1;end;
if nargin<3,start=[0 0];end
lens=lens(:);angle=angle(:);
th=angle*pi/180;
x=cumsum([start(1);lens.*cos(th)]);
y=cumsum([start(2);lens.*sin(th)]);
for i=1:length(x)-1
linkshape([x(i) y(i)],[x(i+1) y(i+1)],w);
end



若使用桿長及水平角度輸入,則可使用下面之plot_angle程式進行,設桿長為5,3,4,6個單位,其對應水平角分別為30, 60, 90, 145度,則可執行如下:

>>link_angle([3 4 6],[60 90 145],[4.2,2.5])

執行結果如下圖:

2.6 機器人的軸節

利用Matlab可以模擬工業機器人軸節之運動情形。下面為一個執行例,以 robot IRB 140 (ABB) 為模擬對象,可以看到其軸節間之運動情形。



就目前所學之平面二軸之迴轉,亦可仿不同軸節之機器人動作。程式link_robot函數可以依輸入之桿長及對應角度與各桿之迴轉速度繪出其連續性之運動。本程式並呼叫linkxy求得各桿中間之座標,然後進行繪製。其程式內容如下:


function link_robot(lens,angle,omega,w)
%link_robot.m will call linkspace.m
%draw connecting links as a connected mechanism.
%Input:lens:lengths of links
% angle:link angles w/t the former link, degs
% omega:angular velocity, rad/s
% w:width of links
% time:[stop delt]
%Example:
% link_robot([5 3 4 6],[30 60 90 145],2:5)
%Author:DSFon, BIME,NTU. Date:Jan. 25, 2007
if nargin<4,w=1;end
if nargin<3,omega=1;end;
lens=lens(:);angle=angle(:);omega=omega(:);
if length(omega)==1,
omega=ones(size(lens))*omega;
end
t=0;nn=length(lens);clf;
for i=1:nn,
h(i)=line('xdata',[],'ydata',[],'erasemode','xor');
end
mm=sum(lens)-min(lens);
axis([-mm mm -mm mm]);axis equal;axis off;
title('Hit ctr-C to stop');
while 1
th=cumsum(angle+omega*t)*pi/180;
x=[0;cumsum(lens.*cos(th))];
y=[0;cumsum(lens.*sin(th))];
for i=1:nn
[xi,yi]=linkxy([x(i) y(i)],[x(i+1) y(i+1)],w);
set(h(i),'xdata',xi,'ydata',yi);
end
drawnow;
pause(0.1)
t=t+1;if t>1000,t=0;end
end

function [x,y]=linkxy(A,B,d)
% Draw a link
% Inputs:
% A,B:Initial & final coordinates of link
% d:thickness of link
% Example: linkshape([0 0],[5,5])
if nargin==2,d=1;end;
d=abs(d);
AB=(B(1)+j*B(2))-(A(1)+j*A(2));
D=abs(AB);th=angle(AB);
t=linspace(pi/2,2.5*pi,20);
Cout=max(d/2,0.2)*exp(j*t');Cin=Cout/2;
if d>0,
P=[0;Cin;Cout(1:10);D+Cout(11:20);D+Cin;D+Cout(20);Cout(1)];
else
P=[Cin;0;D;D+Cin];
end
xx=real(P);yy=imag(P);
x=xx*cos(th)-yy*sin(th)+A(1);
y=xx*sin(th)+yy*cos(th)+A(2);

執行上述程式可得到下面影像之動作:

>> link_robot([5 3 4 6],[30 60 90 145],2:5)


1/20/07

第一章 向量與座標

電腦輔助軟體


在沒有電腦輔助軟體以前,機構之繪製與解析均必須仰賴手繪,利用尺及圓規等工具繪圖,並以變動之角度瞭解機構之運作情形。有了電腦之協助,這項工作大為減輕。且由於電腦軟體之發展,其繪製之能力也大為提升,有些甚至捨棄大量之計算完全由程式代勞,此時學習一種電腦工具更形重要。唯任何要借助軟體工具者在學習上必須時刻警覺,否則極易被軟體之走向及其方便性迷失個人之分析能力。為此仍必須體認電腦僅是分析工具之事實,培養自己的知識及設計能力仍然相當重要。

圖形之繪製,若佐以手工,因投影幾何及工程畫之背景,自可操作自如,並瞭解點線面之形成及物體間之相互關係,但在電腦分析過程中,圖形必須轉換為某特定之座標系,依此座標系之規範,方能在監視器上適當地顯示。為維持圖形在各方向之可看性,圖形物件之座標也常需轉換,以觀看物件之不同角度。

一般圖形包括二維與三維的空間度,其所需的座標轉換動作有移動、放大縮小、反轉、剪切及旋轉等,這些統稱為座標轉換。在計算上,這種轉換工作可利用矩陣來表示。由於MATLAB專於處理矩陣變數,故應用在這方面的運算,不但效率高,其速度也快。

有關Matlab之應用,本節不作詳細論述。同學必須在能力上自我要求。在本部落格之相關教學網頁中,有"Matlab在工程應用"之專題,其中有比較詳細的介紹。同學可以自行參考。其中第一章Matlab之簡介第二章矩陣之製作第九章繪圖等較為重點,同學在研讀本課程時,不可不讀。

機動學主題在討論元件之運動狀態,故其大小尺寸及形狀都相當重要。繪出某一元件,某一機構,常是一個工程人員必須具備的技術,也是未來要在職場上吃飯的傢伙。正如一個樂鼓手必需熟悉其樂器,畫家必須知道怎樣繪畫一樣。但很多人在學習當中,往往忘記自己所需的吃飯的傢伙的重要性,不肯花工夫學習,若有其所學也不深。所以要他表示一個觀念、構思一張圖、繪出一件機構,往往十分為難。不肯努力下功夫,只好在腦海中自己畫餅充饑,幻想未來的職場可能的風光。

畫家繪圖憑靠靈感,工程人員畫圖也要靠靈感。前者可能依徒手的技巧,後者則需要靠尺、圓規等器具之協助以及運算軟體來運算與規畫。在這種情況下,一個工程師所需要的工具較多,而且幾乎都要有所專。小學生繪圖有小畫家,其功能僅是塗鴉的工具,真正使用作為工程圖的則必須有如AUTOCAD、CADKEY、PROENGINEER等類的軟體。工欲善其事,必先利其器。不過上述的軟體都是昻貴的,耗費相當巨大,學習時間也很長。因此若要繪製機動學合用的圖形,使用MICROSOFT附在OFFICE內之VISIO應該比較合用的。當然利用MATLAB除運算外,也可繪圖,一舉兩得。

座標系之轉換


座標之轉換,因常視覺角度之不同,必須進行。在視覺上,有些雖存在有轉換的動作,但因圖隨心轉,有時也是順理成章,反而認為輕而易舉,但實際之應用仍需藉助電腦之運算。這些運算則需以相關之數學為基礎。以設二維座標系之轉換為例,設原物件之座標群座落於(x,y)之體系中,今令其轉換至另一座標系(x',y'),其間之關係可以利用下面的通式表示:

x'= a1x + b1y + c1
y'= a2x + b2y + c2-----------1.1

三維座標系則可表示如下:

x'= a1x + b1y + c1z + d1
y'= a2x + b2y + c2z + d2
z'= a3x + b3y + c3z + d3------1.2

若改以矩陣表示,則二維及三維者應可分別表示如下:

[x' y' 1]=[x y 1] [ a1 a2 0
b1 b2 0
c1 c2 1 ]

[x' y' z' 1]=[x y z 1] [a1 a2 a3 0
b1 b2 b3 0
c1 c2 c3 0
d1 d2 d3 1 ]

故中間為一個特性矩陣,設為A,則A須為方矩陣,最後欄填以0或 1,以利矩陣運算,如此代表不同的轉換型式。

1.1.1 平移(Translation)

平移是沿x軸或y軸方向分別平移Tx、Ty之量,換言之,經過平移後之新座標x',y'應有下列的關係:

x'= x + Tx
  y'= y + Ty---------1.5

由此一聯立方程式可知,在運算上,只要將(x,y)之對應座標值分別加上Tx,Ty即可,例如:
若有座標點A,B,C,D分別為(2,4),(3,9),(4,16),(5,25),則在MATLAB中,可以將其分成兩組對應之座標值(x,y),則用列矩陣表示如下:

>> x=[2 3 4 5]
x =
2 3 4 5
>> y=[4 9 16 25]
y =
4 9 16 25

這兩組座標是對應的,屬於二維之座標。故若要將此組座標分別移(5,10)之距離,則可以就x,y組合之座標值處理如下:

>> x=x+5
x =
7 8 9 10
>> y=y+10
y =
14 19 26 35

兩組座標顯然分別加上各別的距離(5,10)。

然而在矩陣的表示法中,有時候反而比較麻煩,然而利用矩陣卻能將座標之轉移,作成特定的格式,放諸四海而皆準,故也有其特殊的意義。因此,配合公式2.3,可以得到下面類似型式:

 [x' y' 1] = [x y 1] [1 0 0
0 1 0
Tx Ty 1 ]-----1.6

公式2.6之最後矩陣稱為轉換矩陣,若以A表示,則:

[A] = [1 0 0
0 1 0
Tx Ty 1 ]----1.7

設原座標矩陣[C]= [ x y 1 ]及轉換後之座標矩陣[C’]= [x' y' 1]分別代表轉換之列矩陣值,則:

[C'] = [C] [A] ------1.8

式中,

[C] = [x y 1]
[C'] = [x' y' 1]

同理,由公式1.4得知:三維之平移轉換方程式為:

x' = x + Tx
y' = y + Ty
z' = z + Tz -------------1.9

利用相同的觀念,實際上亦可以將三維座標分成三部份,然後就其各部份加上特定係數,終會獲得結果。但同樣,利用矩陣相乘的模式,亦可處理三維座標值之平移的型式,因此其相對之平移轉換矩陣型式可以改寫如下:

[x' y' z' 1]=[x y z 1][1 0 0 0
  0 1 0 0
  0 0 1 0
  Tx Ty Tz 1 ]------- 1.10

在等式之右側為原來之座標系與一方型矩陣相乘,此轉換矩陣A即為:

[A] = [1 0 0 0
  0 1 0 0
  0 0 1 0
  Tx Ty Tz 1 ]------1.11

設原座標矩陣[C] = [ x y z 1]及轉換後之座標矩陣[C’] = [ x' y' z' 1],分別代表轉換之列矩陣值。則其關係為C'=CA。

將上述的關係寫成MATLAB之函數,其程式如下:

程式:

function xprime=trans4T(C, T)
%
%function xprime=trans4T(C,T)
%To transform the translational coordinates
%Inputs:
% C: the row matrix to be transformed
% T: increment of transformation [Tx Ty Tz]
% xprime: resultant matrix transformed
% Example: xprime=trans4T([2 3 4],[2 4 6])
%Author: DS Fon, Bime, NTU, Date: Jan. 22,2007
[n,m]=size(C);
if length(T)~=m,
disp('The dimension is inconsistent!');
xprime=[];
return;
end
C=[C,ones(n,1)];
A=eye(m+1);
A(m+1,:)=[T 1];
xx=C*A;xprime=xx(:,1:m);
%%%%%%%%%%%%%%%%%%%%%%%%


初次碰上這種程式者,可能要由程式之寫作過程上另外加把勁,但寫程式重在其邏輯思路,除熟悉指令內容外,必須將每一步驟之思路想通才行。這個程式很短,凡是程式行前面加上%者,都是註釋行,程式是不執行的,所以只當作程式說明而已。所以真正程式本身,才十來行而已。記得在執行MATLAB程式寫作時,每一個變數是以矩陣看得,所以其代表的至少為1X1的矩陣。

以上面的例子,C可能為一群座標點,它可能代表數千點的座標點,但執行時,卻可以將整體座標同時改變。MATLAB程式函數trans4T(c,T)之功能可以完成上項應用在平移的情況。

實例一:

>> trans4T([3 5 1],[2 2 2])
ans =
5 7 3

其如果與3+2=5, 5+2=7, 1+2=3之情況相同。事實上如果僅這麼簡單,那何必動刀動槍,如此大費週章?問題是大部份的例子並不如此簡單,正如前文所述,第一項的輸入值往往不僅代表這樣的一點座標而已。

實例二:設c1為一個三角形之三個頂點座標組成,繪出三角形時,必須採用連連看的方法,由起點開始,逐點連線,最後還要連回起點,故座標點應有四點才能連出一個三角形。

>> c1=[0 0;1 3;4 2;0 0]
c1 =
0 0
1 3
4 2
0 0
>> line(c1(:,1),c1(:,2)) %這是連線的指令,其第一項需為X座標第二項為Y座標。
>>c2=trans4T(c1,[2 3]) %將原來的點移動一個距離。
c2 =
2 3
3 6
6 5
2 3
>> line(c2(:,1),c2(:,2)) %再次繪出移動後之三角形。
>> grid on %設定為格線。




圖2.1 座標之平移

函數tans4T()因而可以平移各點,並據以繪圖。特別注意的函數中間的參數C是多維向量,其行數依二維或三維而定,不能混用。Line(x,y)則為繪直線函數,可依x,y之對應值作點間連線。

1.1.2 原點放大縮小(scaling)

對於某些點座標若需進行比例放大或縮小,其影響在特性矩陣中,正好是對角線之元素位置,而放大比例依其放大係數S值而定。故S值可以分為在各座標方向如x與y之分別放大量。在A轉換特性矩陣中,其對角元素均為一,若要進行放大縮小則需改變對角項1的位值。因此,只要將對應之1值改為S之對應值Sx, Sy(或Sz)即可;而最後一個1則維持原狀不變化,即:

二維:

x' = (Sx) x
y' = (Sy) y ----1.12

[A] = [Sx 0 0
0 Sy 0
0 0 1 ]----1.13

三維:

x' = (Sx) x
y' = (Sy) y
z' = (Sz) z ----1.14

[A] = [Sx 0 0 0
0 Sy 0 0
0 0 Sz 0
0 0 0 1 ]----1.15


由 12式及14式可知,座標值之放大縮小實際上只要將其對應之座標各乘以對應之Sx,Sy,Sz即可,在正常的運算上也可如此做。例如:有四個空間點之座標需各放大為[0.5 1 2] 倍:

>>x=[2 3 4 5];y=[-1 2 -4 5];z=[10 24 30 5];
>>X=x*0.5; Y=y*1; Z=z*2;
>> X
X =
1.0000 1.5000 2.0000 2.5000
>> Y
Y =
-1 2 -4 5
>> Z
Z =
20 48 60 10


利用MATLAB指令執行結果如上述。
若要採用一致性的矩陣表示法,亦可採用第13及 15式之矩陣相乘之型式,此時可令[C]=[x y z 1]及[C']=[x' y' z' 1],代入公式1.8即可求得對應之放大值。

  [C'] = [C][A]---------------1.8


利用MATLAB程式函數trans4E(c,S)可以完成放大縮小之應用,此時之輸入倍數則以factor表示,其程式內容如下:

function xprime=trans4E(C, factor)
%
%function xprime=trans4E(C,factor)
%To scale up & down the coordinates by a factor
%Inputs:
% C: the matrix to be transformed,[x y] or [x y z]
% factor: scaling ratio, [sx sy] or [sx sy sz]
% xprime: resultant matrix transformed
% Example: xprime=trans4T([2 3 4],[0.1 0.2 0.3])
%Author: D.S. Fon, Bime, NTU, Date: Nov.20,2004
[n,m]=size(C);
if length(factor)~=m,
disp('The dimension is inconsistent!');
xprime=[];return;
end
C=[C,ones(n,1)];
A=eye(m+1); A(m+1,:)=[factor 1];
xx=C*A;xprime=xx(:,1:m);


實例:

>> trans4E([3 5 1],[2 2 2])
ans =
6 10 2

其結果與3*2=6, 5*2=10, 1*2=2的運算相同。
實例:多點作放大時,其結果如下述。

>> c1=[0 0;1 3;4 2;0 0]
c1 = 0 0
1 3
4 2
0 0
>> line(c1(:,1),c1(:,2))
>> c2=trans4E(c1,[2 3])
c2 = 0 0
2 9
8 6
0 0
>> line(c2(:,1),c2(:,2))
>>grid on

圖1.2 座標之放大與縮小

1.1.3 反轉(Inversion)

座標值之反轉即是一種鏡射的原理,利用其對某一軸為參考線,作上下或左右之對稱變化。在計算上,即是將該原矩陣A之相對應項乘以-1即可,不為負值者即對該軸不產生鏡射。此種情形可直接利用放大之功能,將[Sx Sy Sz]之放大系數值改用1代替,而由其正負值決定是否產生反轉之功能。故利用trans4E()可以執行下列之實例:

實例:以Y軸為對稱反轉




>> c1=[0 0;1 3;4 2;0 0]
c1 =
0 0
1 3
4 2
0 0
>> c2=trans4E(c1,[-1 1])
c2 =
0 0
-1 3
-4 2
0 0
>> line(c1(:,1),c1(:,2))
>> line(c2(:,1),c2(:,2))
>> axis equal;grid on



圖1.3 座標之以Y軸鏡射或反轉



實例:以X軸為對稱反轉




>> c1=[0 0;1 3;4 2;0 0]
c1 =
0 0
1 3
4 2
0 0
>> c2=trans4E(c1,[1 -1])
c2 =
0 0
1 -3
4 -2
0 0
>> line(c1(:,1),c1(:,2))
>> line(c2(:,1),c2(:,2))
>> axis equal;grid on



圖1.4 座標之以X軸鏡射或反轉


若三維座標之反轉,程式中之factor參數可以決定如下:
Factor=[ -1 -1 -1]:對原點反轉。
Factor=[ 1 -1 -1]:對X軸反轉。
Factor=[ -1 1 -1]:對Y軸反轉。
Factor=[ -1 -1 1]:對Z軸反轉。

1.1.4 剪切變形(Shearing)

座標值之剪切變形,是因對某特定面為界,進行相對移動。例如在X方面作剪切時,其方程式應為:

x' = x + (Sy) y
y' = y -----------1.16

[A] = [ 1 0 0
Sy 1 0
0 0 1 ]-----1.17

在Y方面作剪切時,方程式為:

x' = x
y' = y + (Sx) x + y -----------1.18

[A] = [ 1 Sx 0
0 1 0
0 0 1 ]-----1.19

三維部份之剪切問題則可由下式表示:

x' = x + (Sxy) y + (Sxz) z
y' = (Syx) y + y + (Syz) z
z' = (Szx) x + (Szy) y + z -----------1.20



其對應之特性矩陣為:

[A] = [ 1 Syz Szx 0
Sxy 1 Szy 0
Sxz Syz 1 0
0 0 0 1 ]----------------1.21

上述公式可利用trans4Sh()函數求得,其程式內容如下。其中shear為剪切之參數,屬行矩陣。其內容在2D時為[Sxy;Syx];3D時為 [Sxy Sxz;Syx Syz;Szx Szy]。但其內容項可為零值,若在3D時,後面項不埴滿時,程式會自動設為零。

function xprime=trans4Sh(C,shear)
%
%function xprime=trans4Sh(C,shear)
%To shear one axis of the coordinates
%Inputs:
% C: the row matrix to be transformed
% shear: shear factor in column,eg [Sxy;Syx]for 2D
% [Sxy Sxz;Syx Syz;Szx Szy]for 3D
% xprime: resultant matrix transformed
% Example: xprime=trans4Sh([2 3],[2;0])
%Author: D.S. Fon, Bime, NTU, Date: Nov.20,2004
[n,m]=size(C);C=[C,ones(n,1)];A=eye(m+1);
shx=zeros(3,2);[k,p]=size(shear);
shx(1:k,1:p)=shx(1:k,1:p)+shear;
A(2,1)=shx(1,1);A(1,2)=shx(2,1);
if m>2
A(3,1)=shx(1,2); %[1 Sxy Sxz]]' for x
A(3,2)=shx(2,2); %[Syx 1 Syz]' for y
A(1:2,3)=shx(3,:)'; %[Szx Szy 1]' for z
end
xx=C*A;xprime=xx(:,1:m);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

實例:以X軸剪切



>> c1=[-4 -3;4 -3;0 3;-4 -3]
c1 =
-4 -3
4 -3
0 3
-4 -3
>> c2=trans4Sh(c1,[2; 0])
c2 =
-10 -3
-2 -3
6 3
-10 -3
>> line(c1(:,1),c1(:,2))
>> line(c1(:,1),c1(:,2))
>> line(c2(:,1),c2(:,2))
>> axis equal;grid on



圖1.5 座標之以X軸進行剪切

實例:以Y軸剪切



>> c2=trans4Sh(c1,[0;1])
c2 =
-4 -7
4 1
0 3
-4 -7
>> line(c2(:,1),c2(:,2))
>> line(c1(:,1),c1(:,2))
>> axis equal;grid on


圖1.6 座標之以Y軸進行剪切


至於3D之實例,讀者可自行印證之。程式中shear之參數決定剪切之方位。

1.1.5 迴轉(Rotating)

對某特定軸作迴轉,在機構學中是常見的例子。若屬平面旋轉,迴轉座標之計算較為簡單,若屬不規則方向之迴轉,其運算較為複雜。因此,迴轉的動作常因迴轉中心的位置而不同,這裡我們僅討論以原點為轉動中心的例子。在二維中,迴轉可僅考慮垂直於xy面的方向,故僅有一種,三維者則分x-、y-與z-及三者兼具的轉動。茲分別說明如下:

二維:





在二維座標之迴轉,以兩座標間之夾角θ決定之。設原座標為(x,y),同一點P經座標旋轉一個角度θ後,其新座標為(x',y')。由圖中可以獲得下列關係式:

x'= xcosθ+ysinθ
y'=-xsinθ+ycosθ

上式為實際可用的關係式,利用MATLAB也可以直接求得解,例如:

>> x=[2 3 4 5];y=[-1 2 -4 5];
>> xprime=x*cosd(30)+y*sind(30)
xprime =
1.2321 3.5981 1.4641 6.8301
>> yprime=-x*sind(30)+y*cosd(30)
yprime =
-1.8660 0.2321 -5.4641 1.8301

上式中,係原座標x=[2 3 4 5];y=[-1 2 -4 5]經旋轉30度後之座標值。三角函數之參數若使用度數,可利用sind,cosd等函數可必經過轉換為弧度;然而若使用sin,cos則必須事先轉化為弧度才能得到正確的結果。

上式只是簡單的二維計算,若屬三維且轉軸的方向較為複雜時,仍然以採用矩陣乘法為佳。其型式可以轉換如下:


[x' y' 1]=[x y 1][A(3x3)];

A(3x3)=[ cosθ sinθ 0
 -sinθ cosθ 0
  0 0   1]

C =[x y 1];
C'=[x' y' 1];

C' = C A;



三維:


三維中之特性矩陣則依旋轉軸不同如下:

以x軸旋轉:



A(4x4)=[1 0 0 0
0 cosθx sinθx 0
0 -sinθx cosθx 0
0 0 0 1 ]

以y軸旋轉:



A(4x4)=[cosθy 0 sinθy 0
0 1 0 0
sinθy 0 cosθy 0
0 0 0 1 ]

以z軸旋轉:



A(4x4)=[ cosθ sinθ 0 0
-sinθ cosθ 0 0
0 0 1 0
0 0 0 1 ]

以y-, x-, x- 軸複合旋轉:


可以利用上述三個方向連續轉動。

計算得A矩陣後,即可逕行計算新座標如下:
C=[ x y z1 ];
C'=[x' y' z' 1];
C' = C A;

在本節後面所附的函數在MATLAB程式trans4函數是可以就上述不同的情況處理平移、放大及迴轉等問題。在此程式trans4(C, r, mod)中,利用其最後輸入項mod作為控制碼,mod=1時為平移;mod=2時為放大縮小;mod=3mod=4時為mod=5時則為迴轉。因此其對應r值為處理其方矩陣之參數。當mod=5時為迴轉問題,其參數r即為欲旋轉之角度[thetax thetay thetaz];若僅為二維為僅能使用一個參數,即[thetaz]:

例一、以一箭號旋轉z軸六次



>> c=[0 0;0 20;-3 4;3 4;0 20] %製作迴轉物件
c =
0 0
0 20
-3 4
3 4
0 20
>> hold on;
>> for theta=1:6,c=trans4(c,[60],5);line(c(:,1),c(:,2));end; %連續迴轉60度
>> axis equal;
>> grid on;


平面迴轉

例二、進行立體圖以Y軸迴轉



clf; %將前面的各圖清除
c1=[0 0 0;1 0 0;0 1 0;0 0 0;0 0 2;1 0 0;0 1 0;0 0 2]; %製造迴轉用的立體物件
hold on; %連續繪圖模式
for theta=0:30:360, %對Y軸旋作等級旋轉
c=trans4(c1,[0 theta 0],5); %轉換座標
line(c(:,1),c(:,2),c(:,3)); %進行繪圖線
end;
grid on;
xlabel('x'); %X-座標名稱
ylabel('y');
zlabel('z');


立體迴轉

trans4程式內容



function xprime=trans4(x,delx,mod)
%
%function xprime=trans4(x,delx,mod)
%For dimensional transformation
%mod indicates the type of transformation
% mod=1 for translation
% mod=2 for scaling(enlarging or shrinking)
% mod=3 for inversion delx=[-1 -1 -1]
% mod=4 for shearing delx=[xa ya za]
% mod=5 for rotation delx=[angx angy angz]
% written by D.S. Fon, Bime, NTU
%
[n,m]=size(x);
x=[x,ones(n,1)];
tt=eye(m+1);
switch mod,
case 1 %translation
tt(m+1,:)=[delx 1];
xx=x*tt;
case {2,3} %enlarging & inversion
tt(tt==1)=[delx 1];
xx=x*tt;
case 4 %shear
if m<3,
tt(2,1)=delx(1);
tt(1,2)=delx(2);
else
if size(delx,1)==1,
delx(2,:)=0;
end;
tt(2:3,1)=delx(1:2,1);%[1 d g 0]' for x
tt(1,2)=delx(1,2); %[b 1 0 0]' for y
tt(3,2)=delx(2,2); %[b 1 h 0]' for y again
tt(1:2,3)=delx(1:2,3);%[c f 1 0]' for z
end;
xx=x*tt;
case 5
d2g=pi/180;
theta=delx*d2g;
if m<3,
t=theta(1);
sint=sin(t);cost=cos(t);
tt(1,1)=cost;tt(1,2)=sint;
tt(2,1)=-sint;tt(2,2)=cost;
xx=x*tt;
else
tt0=tt;
if abs(theta(1))>0, % rotate with x-axis
cost=cos(theta(1));sint=sin(theta(1));
tt0(2,2:3)=[cost sint];
tt0(3,2:3)=[-sint cost];
x=x*tt0;
end;
if abs(theta(2))>0, % rotate with y-axis
cost=cos(theta(2));sint=sin(theta(2));
tt0=tt;
tt0(1,1:3)=[cost 0 -sint];
tt0(3,1:3)=[sint 0 cost];
x=x*tt0;
end;
if abs(theta(3))>0, % rotate with z-axis
cost=cos(theta(3));sint=sin(theta(3));
tt0=tt;
tt0(1,1:2)=[cost sint];
tt0(2,1:2)=[-sint cost];
x=x*tt0;
end;
xx=x;
end;
end
xprime=xx(:,1:m);

1.1.6 內建之rotate指令

MATLAB也有內建之迴轉指令rotate,該指令係針對某一特定方向作物件旋轉。其語法如下:

rotate(h,direction,alpha)
rotate(...,origin)

迴轉函數是將一個三度空間的物件依右手法則進行旋轉。這個指令並無輸出,僅就h的圖形直接改變資料。

rotate(h,direction,alpha) 之敘述是將物件h(握把)作alpha 角度的迴轉,而direction則是以原點為參考點之轉動軸,通常以二維或三維向量表示。

rotate(...,origin)則是特別指出迴轉軸之原點,以三維向量表示。預設原點為繪圖方塊(plot box)之中心點。



須注意的是待轉動的物件應同屬座標軸的子系,物件之資料會因旋轉過程而改變,這與view 及rotate3d兩指令不同,後者僅是改變外在的觀測方位。

迴轉軸是經過原點與另一點P來決定。P點可用球座標[theta phi]或直角座標表示。球座標係以兩項參數[theta phi]決定旋轉方向,theta為xy面上由正x軸算起反時鐘方向為正;phi則是由xy面上算起的高度。

直角座標則以三項參數(X,Y,Z)表示,迴轉方向則是由原點指向這一點的座標。一般而言,rotate這個指令會改變圖形中之 Xdata, Ydata與 Zdata 等之參數。


例一、將一物件圖沿x軸迴轉180º



% demo_trans4_1
subplot(1,2,1); h1 = surf(peaks(20));% original plot
xlabel(’Orginal plot’);
subplot(1,2,2); h2= surf(peaks(20)); %duplicate it for comparison
rotate(h2,[1 0 0],180);
xlabel(‘Rotate(h2,[1 0 0],180)’);



二、將一個長方形的物件沿z軸迴轉45º、90º、135º、180º、225º



%demo_rot3.m
clf;
h = surf(peaks(20));
zdir = [0 0 1];
center = [10 10 0];
hold on;
axis equal;
xlabel('x-axis');
ylabel('y-label');
zlabel('z-label');
for i=0:45:225,
rotate(h,zdir,i,center);
pause;
end;



利用上面程式執行後,可以看到peaks之圖形依z軸方向作迴轉的情形。其中之pause指令則是依執行者按鍵回應而逐步完成整套圖形之轉換。