Showing posts with label 向量. Show all posts
Showing posts with label 向量. Show all posts

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

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函數中之畫圖之顏色?

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指令則是依執行者按鍵回應而逐步完成整套圖形之轉換。