平面與空間曲線曲率研究與軌跡動畫模擬
指導教授：陳正宗終身特聘教授 施佑勳 高聖凱 李家瑋
國立台灣海洋大學 河海工程系
關鍵詞：Frenet formula，Mathematica ，弧長參數表示法，曲率半徑，扭率
摘要
本論文主要探討平面與空間曲線。由向量微積分的觀點，利用曲線弧長與時間參數表示法作為切入點，引入Frenet formula 的手法求得空間曲率。針對已知曲線起始點、起始架構與相關曲率之限制條件下，透過Frenet formula轉換成具初位置與初架構(frame)之聯立微分方程組求解問題，此聯立方程組可寫成狀態方程，並用e 計算。最後運用Mathematica 軟體來建構軌跡動畫。相信本文對日後學生學習、老師教學與工程運用都會有相當的幫助。
A t
前言
在日常生活中，一般道路在直線段並無曲率或稱曲率半徑無限大。但在於曲線路段時，將產生離心力，使車輛滑出車道，並造成翻車。所以車輛由直線段進入曲線段在此兩路段間必需逐漸調整方向，而在圓曲線路段內車輛有一定的離心力，在公路設計時亦需要有一路段區間以供調整路面的傾斜度以平衡車輛離心力，再則曲線路段內車輛後輪將向圓心方向偏移，所以圓曲線路段需要加寬，這路線方向的調整、路面傾斜度的漸變以及路面加寬通常會由曲率來決定，超高設計是與曲率有關的。在數學應用上，通常使用微積分技巧求曲線之曲率，而曲率主要是用來描述曲線上某處，曲線彎曲變化的程度。當曲線軌跡推至三維後，描述線段軌跡的因素則多了一項扭率。在材料力學梁橈曲中[1]，計算正向應力時的推導，亦需使用到曲率觀念。拜現代科技進步所賜，目前學習趨向多元化，不僅有前人推導的相關理論，更可藉由更多的輔助學習工具，增進學習速度和效率。本論文使用了Symbolic software 的Mathematica[2-3]，利用其強大的符號運算能力，更可輸出以動畫模擬。不同以往的靜態教學，增添動畫模擬的輔助，更能吸引初學者令其產生極大之學習興趣。
研究方法與步驟
1. 文獻回顧與探討
1.1 Frenet formula
給一空間曲線，其時間參數表示式為(x (t ), y (t ), z (t )) ，藉由弧長關係式可寫成
2222
(ds ) =(dx ) +(dy ) +(dz ) (1)
若將時間參數表示法轉至空間參數表示法，則
r =(x (t ), y (t ), z (t )) =((s ), (s ), (s )) (2)
其中，r 為位置向量。單位切向量τ，定義為
dr
τ= (3)
ds
其中，ds 為微小弧長。單位法向量ν因與τ正交，
故可令
τ
ν= (4)
單位雙法向量(binormal vector)β，則與τ和ν向量正交。因此可得
β=τ⨯ν (5)
則τ、ν與β三向量的關係，如圖一所示
圖一 切向量τ 、法向量ν 與雙法向量β
由(3)式可知，單位切向量可表示為
τ(s ) =
(s ) r (s ) r
=
1
ρ
(17)
將(17)式代入(4)式，可知 又 則
ν=ρτ (18)
β⋅β=1
=0β⋅β
=p τ+q νβ
(19) (20)
(21)
用(20)式可知，β與β 正交關係，故可令
(6)
因
將
τ(s ) ⋅τ(s +ds ) =(s ) (s +ds ) cos(d θ) =cos(d θ) (7)
τ⋅β=0
利用泰勒展開式將τ(s +ds ) 與cos(d θ) 展開，可得
111⎡⎤ (s )(ds ) 2+τ (s )(ds ) 3+ τ(s ) ⋅τ(s +ds ) =τ(s ) ⋅⎢τ(s ) +τ (s )(ds ) +τ⎥1! 2! 3! ⎣⎦
∞
則
(18)式代入(23)式可得
(22) (23)
+τ ⋅β=0τ⋅β
=τ(s ) ⋅∑
n =0
τ(s )
n !
2
n
(ds ) (8)
n
+τ⋅β
cos(d θ) =1-
12! 12
(d θ) +
14!
(d θ) -
4
16! 1
ν⋅β=0 (24) ρ
1
6
(d θ) + (9)
因ν與β正交，可得
=0τ⋅β
將(8)與(9)式代入(7)式，可得
1+
(s )(ds ) =1-τ(s ) ⋅τ
2
(25)
2!
(d θ)
2
(10)
由(21)式等號兩邊與τ作內積，可得 因τ、ν正交，可得
⋅τ=p τ⋅τ+q ν⋅τ (26) β
p =0 (27) =q νβ
(s )(ds ) 2=-(d θ) 2 (11) τ(s ) ⋅τ
又 可推得
τ⋅τ =0 (12)
=-τ ⋅τ (13) τ⋅τ
所以 因此令
(28)
ν
(13)式代入(11)式可得
22
-τ ⋅τ (ds ) =-(d θ)
(14)
其中
=-β
1
則
2
σ
1 β
(29)
111⎛d θ⎫
= === (15) ⎪222
ρ⎝ds ⎭⎛ds ⎫⎛ρd θ⎫
⎪ ⎪⎝d θ⎭⎝d θ⎭
2
σ=
(30)
σ為扭率，而
ρ=
ds d θ
β=τ⨯ν (31)
(16)
又ν、β與τ相互正交，故
圖二 曲率半徑ρ
由(15)式與圖二可得
ν=β⨯τ
ν =
1
(32)
τρ
1
σ
β-
(33)
整理(18)、(29)與(33)式得狀態空間表示式如下：
⎡⎢0
⎧τ ⎫⎢⎪ ⎪⎢1⎨ν ⎬=⎢-
ρ ⎪⎪β⎢
⎩⎭⎢
⎢0⎣
1
ρ
0-1
σ
⎤0⎥⎥1⎥⎥σ⎥⎥0⎥⎦
2. 矩陣函數技巧與應用[5-17]
⎧τ⎪ ⎨ν ⎪β⎩
⎫⎪
⎬ (34) ⎪⎭
2.1 傳統方法
假設有一矩陣A ，其特徵向量組成的矩陣為V
V =[v 1⋅⋅⋅v n ] (39)
Av j =λj v j ，j =1, ⋅⋅⋅, n (40)
可寫成
AV =VD (41)
令
⎧τ⎪ P =⎨ν ⎪β
⎩
⎫⎪
⎬ (35) ⎪⎭
其中
⎡λ1⎢0
D =diag (λ1, ⋅⋅⋅, λn ) =⎢
⎢0⎢⎣0
00 0
則(34)式可寫成
λ2
00
即為 Frenet formula，其中
⎡⎢0⎢⎢1A =⎢-
ρ⎢⎢⎢0⎣
1
=AP (36) P
⎤
0⎥⎥1⎥
⎥ (37) σ⎥⎥0⎥⎦
0⎤⎥0
⎥ (42) 0⎥⎥λn ⎦0⎤⎥0⎥
(43) ⎥0⎥λt e n ⎦⎥
則
⎡e λ1t
⎢0λn t λ1t
=diag (e , ⋅⋅⋅, e ) =⎢
⎢0⎢⎢0⎣
e
A t
ρ
0-1
0e
λ2t
00 0
e
Dt
00
-1
σ
可看出A 為反對稱矩陣。 1.2 求解過程
計算步驟流程分述如下：
(a) 轉換參數表示法。將時間曲線參數式，轉換成弧 長表示式r (s ) 。
(b) 求單位切向向量τ。
(c) 求曲率半徑ρ。 (d) 求法向量ν。 (e) 求雙法向量β。
(f) 求扭率σ。
(g) 列出聯立微分方程並求解。
1⎧
=ντ⎪
ρ ⎪
11⎪
β⎨ν =-τ+
ρσ ⎪
⎪1
⎪β=-ν
σ ⎩
, τ(0)=τ0
, ν(0)=ν0 (38) , β(0)=β0
所以
=V e V
D t
(44)
2.2 Jordan 正則式(Canonical Form)
假設有一n ⨯n 矩陣A ，其特徵值有二重根，可利用
Jordan Canonical Form解決重根問題，矩陣特徵值有二重根時Jordan Form表示為
⎡λ1⎢0⎢⎢0⎢⎢0⎢0⎣
00 00
000
λ2
000
A =V
λi
0⎤⎥0⎥0⎥V ⎥1⎥λi ⎥⎦
-1
(45)
其中λi 為特徵值，以下為λ四重根
⎡λ
⎢0J =⎢
⎢0⎢⎣0
1
01
λ
00
λ
2
0⎤⎥0
⎥ (46) 1⎥⎥λ⎦
t
3t ⎤⎥3! ⎥2t ⎥⎥2! ⎥t ⎥⎥1⎦
所以
⎡⎢1⎢⎢⎢0⎢⎢0⎢⎣0
t 100
2! t 10
給定初位置(x (0),y (0),z (0))，與初架構(frame) (h) 利用 Frenet formula 求得的解，並將τ、ν與β座標架構隨參數變化以動畫呈現。
(i) 利用Mathematica 符號運算軟體作出三維軌跡動
畫。
(τ(0),ν(0),β(0))
e
Jt
=e
λt
(47)
推導過程如下：
喬登基本矩陣定義如下
⎡λ⎢0J =⎢
⎢0⎢⎣0
1
01
λ
00
λ
0⎤⎥0
⎥ (48) 1⎥⎥λ⎦
⎤1
f '''(λ) ⎥3! ⎥⎥1
f ''(λ) ⎥2! ⎥⎥1
f '(λ) ⎥1! ⎥
則
⎡
⎢f (λ) ⎢⎢⎢0=⎢⎢⎢0⎢⎢⎣0
⎡⎢1⎢⎢⎢⎢0⎢⎢0⎢⎢0⎣
1
f '(λ) 1! f (λ) 00t 100
t 22! t 10
1
f ''(λ) 2! 1
f '(λ) 1! f (λ) 0t 3⎤3! ⎥t 2⎥
⎥2! ⎥t ⎥
⎥⎥⎥
f (J ) =e Jt
f (λ)
⎥⎦
=e
λt
(49)
1⎥⎦
2.3 矩陣餘式定理[18]
給一矩陣A ，其特徵方程式為f (λ) =0，特徵
值為λ1⋅⋅⋅λn 。由Cayley-Hamilton 定理知f (A ) =0。透過餘式定理可知
e λ=f (λ) q (λ) +a n -1λn -1+a n -2λn -2+⋅⋅⋅
+a 1λ+a 0
4. 狀態方程求解
軌跡方程可寫成狀態方程，因此可用e A t 求解。
簡單列出其解題過程，並為接下來的曲線動畫做準備。
例題1. 解如下狀態方程，求切向量τ(s ) 、法向量
ν(s ) 與雙法向量β(s ) 。
(50)
將λ1⋅⋅⋅λn 分別代入(50)式，可求得a n ⋅⋅⋅a 0。再由實數和矩陣可互換性質，將實數換為矩陣可寫成
e At =f (A ) q (A ) +a n -1(t ) A n -1+a n -2(t ) A n -2+⋅⋅⋅
+a 1(t ) A +a 0I
(51)
代入矩陣A 得e 。
A t
已知控制方程
⎡⎢0⎢⎢1⎢-⎢2⎢⎢0⎣
研究結果
3. 比較二維各種曲率表示式
在 Frenet formula 中，扭率若為零的特例情況 下，可將空間曲線簡化為二維平面曲線。各種二維曲線之曲率計算方法整理比較於表一。
⎧
τ (s ) ⎫⎪⎪⎪ ⎪⎨ν (s ) ⎬⎪ ⎪ (s ) ⎪⎪β⎩⎭
1
20-12
=
⎥⎧τ(s ) ⎫⎪ ⎪1⎥⎪⎪ ⎥⎨ν(s ) ⎬2⎥⎪ ⎪
β(s ) ⎪⎥⎪⎩⎭0⎥ ⎦
0⎥
⎤
(52)
初架構(frame)
⎧
⎪τ(0)= ⎪⎪
⎨ν(0)=(1,0, 0) ⎪
⎪β(0)=(0,⎪⎩
(53)
初位置
⎧x (0)=1
⎪
⎨y (0)=0⎪
⎩z (0)=0
(54)
設
⎧τ(s ) ⎫⎪ ⎪
P =⎨ν(s ) ⎬ (55)
⎪β(s ) ⎪⎩⎭ ⎧τ (s ) ⎫ ⎪ =⎪νP ⎨ (s ) ⎬ (56) ⎪β (s ) ⎪
⎩⎭
將(54)式帶入，則空間曲線(x (s ), y (s ), z (s )) 點式表示如下
(x (s ), y (s ), z (s )) =(- (65)
而
5. 曲線動畫與Frenet formula參數研究
利用例題1. 的曲線針對初位置、初架構、曲率與扭率做參數研究。 5.1 曲線動畫
黑線為空間曲線，紅線表示τ(s ) ，綠線表示ν(s ) ，藍線表示β(s ) ，如圖三所示
則 可解得
=AP (57) P ⎧τ(0)⎫
⎪ ⎪
⋅⎨ν(0)⎬ (58) ⎪β(0)⎪⎩⎭
P =e
As
其中，
⎡-1⎢
=⎢-⎢1⎣
-11
⎡⎢1⎤⎢⎥
0⎥. ⎢0
⎢1⎥⎢0⎦⎢⎢⎣
e
As
e 0
⎤⎡-1
0⎥⎢4⎥⎢⎢10⎥. ⎢-⎥40s
e ⎥⎢⎥⎢1⎥⎦⎢⎣2
1⎤
4⎥⎥1⎥4⎥⎥1⎥⎥2⎦
(59)
代入(53)式
⎡⎢2
⎢cos ⎢
⎢⎢⋅C =⎢⎢⎢⎢⎢2
⎢sin ⎢
⎢⎣
⎤
⎥
2sin ⎥
⎥
⎡⎥0⎢
⎥
⎥⎢
. ⎢1⎥⎢⎥⎢⎥0⎢⎣⎥2
cos ⎥
⎥
⎥⎦
圖三 空間曲線圖
(60)
5.2 改變初位置
0e
As
0⎥⎥改變初位置，僅平移改變其初始點如下：
⎧x (0)=1⎧x (0)=2
⎪⎪
(a ) ⎨y (0)=0與(b ) ⎨y (0)=0 ⎪⎪⎩z (0)=0⎩z (0)=0
得
⎡
⎢⎢
⎧τ(s ) ⎫⎢
⎪ ⎪⎢⎨ν(s ) ⎬=⎢⎢ ⎪β(s ) ⎪
⎩⎭⎢ ⎢⎢
⎢⎣
-⎥0
⎥⎥ (61)
黑線為(a ) 與藍線為(b ) 。發現曲線形狀並未改變。
如圖四所示
再把τ(s ) 對s 做積分得
=-x (s ) =
⎰
+C 1
(62)
圖四 改變初位置表示圖
5.3 改變曲率半徑ρ
改變曲率半徑則會影響曲線投影於xy 平面上圓之大小，如圖五(a)與(b)所示
y (s ) =
⎰
s 0
ds =+C 2
(63)
z (s ) =
⎰
s 0
ds =
C 3 (64)
①ρ②ρ③ρ④ρ
=2
α⎫α⎫⎛⎛
x =a cos s cos() ⎪, y =a sin s cos() ⎪, z =s sin(α)
a ⎭a ⎭⎝⎝
(66)
=2.4 =2.6
=2.2
a 為半徑，α為角度，s 為變數，如圖七所示
⑤ρ=3.0 ⑥ρ=4.0
圖五(a) 改變曲率半徑之曲線(σ=2)
圖七 空間曲線之a 半徑與α角度
6.1 改變半徑
半徑控制曲線投影於xy 平面上形成圓半徑之大小， 半徑越大圓越大，如圖八(a)與(b)所示，其中固定(角度) α=
π
4
圖五(b) (σ=2)
(τ(s ), ν(s ), β(s ))
，紅線半徑為a =1，綠線半徑為a =2，
藍線半徑為a =4。
5.4 改變扭矩σ
改變σ扭矩常數為影響曲線之除峭程度，σ越
大曲線較為平緩，如圖六(a)與(b)所示
①σ②σ③σ④σ⑤σ⑥σ
=2 =2.2 =2.4 =2.6 =3.0
紅線a =1 綠線a =2 藍線a =4
圖八(a) 改變半徑之曲線(側視圖) (α=
π
4
)
=4.0
圖六(a) 改變σ之曲線(ρ=2)
圖八(b) 改變半徑之曲線(俯視圖) (α=
π
4
)
圖六 (ρ=2)
6.2 改變角度
角度控制曲線除峭程度，為一曲線之水平角，角度越大曲線越除，如圖九(a)與(b)所示，其中固定(半徑) a =1，紅線角度為α=藍線角度為α=
π
16
(τ(s ), ν(s ), β(s ))
π
4
，綠線角度為α=
π
8
，
6. 空間曲線之時間與弧長參數表示法[19]
曲線時間與弧長參數表示
。
,
所以
紅線α=綠線α=藍線α=
π
4π8
∂R
=r sin(ϕ) ∂θ
(71)
π
16
圖九(a) (a =1)
由於
d
e r =sin(ϕ) cos(θ) i +sin(ϕ) sin(θ) j +cos(ϕ) k (72)
e ϕ=cos(ϕ) cos(θ) i +cos(ϕ) sin(θ) j -sin(ϕ) k (73)
e θ=-sin(θ) i +cos(θ) j +o k (74)
d ϕd θ
(e r ) =cos(ϕ) cos(θ) i -sin(ϕ) sin(θ) i dt dt dt +cos(ϕ) sin(θ)
d θd ϕ
j +sin(ϕ) cos(θ) j -sin(ϕ) k dt dt dt
d d ϕd θ(e ϕ) =-sin(ϕ) cos(θ) i -cos(ϕ) sin(θ) i dt dt dt
j +cos(ϕ) cos(θ) j -cos(ϕ) k
dt dt dt
d d θd θ
(e θ) =-cos(θ) i -sin(θ) j +0k
dt dt dt
d ϕ
d θ
d ϕ
d ϕ
(75) (76)(77)
-sin(ϕ) sin(θ)
圖九(b) (a =1)
得
⎧e r ⎪ ⎪ ⎨e ϕ ⎪ e ⎪⎩ θ
⎫⎡0⎪⎢⎪
-ϕ⎬=⎢
⎪⎢-sin(ϕ) θ ⎪⎭⎣
ϕ0-cos(ϕ) θ
sin(ϕ) θ ⎤⎧e r
⎥⎪
cos(ϕ) θ ⎥⋅⎨e ϕ
⎥⎪e 0
⎦⎩ θ
⎫⎪
⎬ (78) ⎪⎭
7. 三維空間球座標狀態方程矩陣討論[20]
7.1 球座標
向量R
，長度r 如圖十所示
上式狀態方程矩陣亦是反對稱。 7.2 圓柱座標
空間上任一點(x , y , z ) 可以表示為圓柱座標
(r θ, z , ，如圖十一所示)
y
⎧x =r sin(ϕ) cos(θ)
⎪
⎨y =r sin(ϕ) sin(θ) ⎪z =r cos(ϕ) ⎩
x
圖十 球座標定義 則
(67) (68)
圖十一 圓柱座標定義
R =(x , y , z ) =(r sin(ϕ) cos(θ), r sin(ϕ) sin(θ), r cos(ϕ))
分別對r , θ, ϕ偏微
R =r sin(ϕ) cos(θ) i +r sin(ϕ) sin(θ) j +r cos(ϕ) k
∂R
=sin(ϕ) cos(θ) i +sin(ϕ) sin(θ) j +cos(ϕ) k ∂r
∂R
, =1 (69) ∂r
∂R
=r cos(ϕ) cos(θ) i +r cos(ϕ) sin(θ) j -r sin(ϕ) k ∂ϕ
R =(x , y , z ) =(r cos(θ), r sin(θ), z ) (79)
令位置座標R =(x , y , z ) =(r c o s θ(), r s i θn (), z ，取方)
向導微
∂R ∂R
=cos(θ) i +sin(θ) j +0k , =1 (80) ∂r ∂r ∂R ∂R
=-r sin(θ) i +r cos(θ) j +0k , =r ∂θ ∂θ
∂R
, =r ∂ϕ
(70)
則
(81)
∂R
=-r sin(ϕ) sin(θ) i +r sin(ϕ) cos(θ) j +0k ∂θ
∂R ∂R
=0i +0j +1k , =1 (82) ∂z ∂z
台北市, 文魁 (2002).
[4] 陳正宗, 工程數學(二) 講義, 海大河工力學聲響振
動實驗室 (2008).
[5] Arieh Iserles and Antonella Zanna, “Efficient
所以 computation of the matrix exponential by
generalized polar deco mpositions, ” SIAM Journal r =-sin(θ) θ i +cos(θ) θ j +0k (86) e
on Numerical Analysis, Vol. 42, pp. 2218-2256
θ=-cos(θ) θi -sin(θ) θj +0k (87) [6] J.T. Chen, S.R. Kuo and C.F. Lee, “A new point of e
view for the Householder matrix by using matrix
z =0i +0j +0k (88) e
exponential, ” International Journal of Applied
整理得 Mathematics , Vol.7, No.3, pp.289-308 (2001).
[7] Kenneth G . Wilson and Michael E. Fisher, “Critical r ⎫⎡0θ 0⎤⎧e r ⎫⎧e
Exponents in 3.99 Dimensions, ” Physical Review ⎥⎪ ⎪⎪ ⎪⎢
θ⎬=⎢-θ00⎥⋅⎨e θ⎬ (89) ⎨e Letters , Vol 28, No.4, pp.240-600 (1972). ⎪ ⎪⎪e ⎪⎢⎥[8] J.S. Kole, “Solving seismic wave propagation in 00⎩e z ⎭⎩ z ⎭⎣0⎦ elastic media using the matrix exponential
上式狀態方程矩陣亦是反對稱。
approach, ” W ave Motion, Vol.38, pp.279-293 (2003) . [9] C. Morler and C. V. Loan, “Nineteen Dubious Ways 結論與建議
to Compute the Exponential of a Matrix, ” SIAM Review , Vol.20, No.4, pp.801-836 (1978). 本文先討論e A t 的兩種計算方法(含傳統與矩陣
[10] M. Vidyasagar, “A Novel Method of Evaluating 餘式定理) ，並用於空間曲線中求Frenet formula的曲
in Closed Form,” IEEE Transactions on Automatic
線軌跡，同時也對Frenet formula中的曲率與扭率做
Control , Vol.1, pp.600-611 (1970).
更進一步的分析及探討，包含曲線動畫、改變曲率
[11] H.-K. Hong, H. S. Lan and J.T. Chen, “On Residue
半徑、改變扭率與改變初位置均一併討論。在求解Theorem of Matrix Function and Its Application to A t
e 矩陣函數時，傳統方法較為一般人所知，但當遇Spin Tensor Using Euler Angle Concept,” (1990).
[12] Thomas J.R. Hughes and James Winget, “Finite 到A 矩陣特徵值有重根時，就必須利用到Jordan 正
Rotation Effects in Numerical Integration of Rate 則式(Canonical Form) 的技巧，而矩陣餘式定理較具
Constitutive Equations Arising in Large 一般性，適合解任一之矩陣函數，而遇到重根問題
Deformation Analysis,” Short Comm (1979).
時可運用微分運算求得所缺少的方程式，是個實用
[13] R.Rubinstein and S.N. Atluri, “Objectivity of
面和適用性較高的方法。接下來透過純座標轉換，Incremental Constitutive Relations over Finite Time 展現Frenet formula獨有的特性，找到對應的反對稱Steps in Computational Finite Deformation 矩陣，最後建立空間曲線之時間與弧長參數表示Analysis, ” Comp. Meth. App. Mech. Engng., Vol.36, 法，對於二維曲線則視A 特例亦有探討。由於動畫pp.277-290 (1983).
[14] D.B. Metzger and R. N. Dubey, “Corotational 模擬的方式可將較為抽象的三維空間曲線隨弧長變
Rates in Constitutive Modeling of Elastic-Plastic 化完整呈現於眼前，引起學習之興趣，故本文不僅
Deformation,” Int. J. Plasticity, V ol.4, pp.341-368 可提供Mathematica 動畫之初學者作為參考，還可以
(1987).
當作教學上之一輔助工具，提升教與學之效能。
[15] Juan C. Simo and Jerrold E. Marsden, “On the
Rotated Stress Tensor and Material V ersion of the 誌謝
Doyle - Ericksen Formula” (1983).
[16] D.P. Flanagan and L.M. Taylor, An Accurate 本研究感謝教育部卓越教學計畫與教學中心特
Numerical Algorithm for Stress Integration with 別規劃「大學生暑期學習計畫」，及國立台灣海洋大
Finite Rotations, Comp. Meth. App. Meth. Eng., 學河海工程學系力學聲響振動實驗室之陳正宗終身
V ol.62, pp.305-320 (1987).
特聘教授以及實驗室學長們的辛勤指導與感謝戴志
[17] R.C.Ward, “Numerical Computation of the Matrix
豪同學在二維特例部分的討論。 Exponential with Accuracy Estimate, ” SIAM J.
Number Anal, Vol.14,No.4, pp.600-610 (1977).
參考文獻
[18] 陳正宗, 矩陣函數的新解法–矩陣餘式定理之
[1] J.M. Gere, Mechanics of Materials, 3rd Edition, 應用, 新新季刊, 中山科學研究院, 龍潭,V ol.15,
Prentice Hall (1997). No.4 (1987). [2] 張起豪,Mathematica 套裝程式使用與設計, 台北[19] R. Aris, V ectors Tensors and the Basic Equations
市, 碁峰資訊 (1994). of Fluid Mechanics, 協成書局 (1983). [3] 余家銘 編著,Mathematica 程式設計風格與應用, [20] 林琦焜, 向量分析, 滄海書局 (2007).
e r =cos(θ) i +sin(θ) j +0k (83)
e θ=-sin(θ) i +cos(θ) j +0k (84)
e z =0i +0j +1k (85)
平面與空間曲線曲率研究與軌跡動畫模擬
指導教授：陳正宗終身特聘教授 施佑勳 高聖凱 李家瑋
國立台灣海洋大學 河海工程系
關鍵詞：Frenet formula，Mathematica ，弧長參數表示法，曲率半徑，扭率
摘要
本論文主要探討平面與空間曲線。由向量微積分的觀點，利用曲線弧長與時間參數表示法作為切入點，引入Frenet formula 的手法求得空間曲率。針對已知曲線起始點、起始架構與相關曲率之限制條件下，透過Frenet formula轉換成具初位置與初架構(frame)之聯立微分方程組求解問題，此聯立方程組可寫成狀態方程，並用e 計算。最後運用Mathematica 軟體來建構軌跡動畫。相信本文對日後學生學習、老師教學與工程運用都會有相當的幫助。
A t
前言
在日常生活中，一般道路在直線段並無曲率或稱曲率半徑無限大。但在於曲線路段時，將產生離心力，使車輛滑出車道，並造成翻車。所以車輛由直線段進入曲線段在此兩路段間必需逐漸調整方向，而在圓曲線路段內車輛有一定的離心力，在公路設計時亦需要有一路段區間以供調整路面的傾斜度以平衡車輛離心力，再則曲線路段內車輛後輪將向圓心方向偏移，所以圓曲線路段需要加寬，這路線方向的調整、路面傾斜度的漸變以及路面加寬通常會由曲率來決定，超高設計是與曲率有關的。在數學應用上，通常使用微積分技巧求曲線之曲率，而曲率主要是用來描述曲線上某處，曲線彎曲變化的程度。當曲線軌跡推至三維後，描述線段軌跡的因素則多了一項扭率。在材料力學梁橈曲中[1]，計算正向應力時的推導，亦需使用到曲率觀念。拜現代科技進步所賜，目前學習趨向多元化，不僅有前人推導的相關理論，更可藉由更多的輔助學習工具，增進學習速度和效率。本論文使用了Symbolic software 的Mathematica[2-3]，利用其強大的符號運算能力，更可輸出以動畫模擬。不同以往的靜態教學，增添動畫模擬的輔助，更能吸引初學者令其產生極大之學習興趣。
研究方法與步驟
1. 文獻回顧與探討
1.1 Frenet formula
給一空間曲線，其時間參數表示式為(x (t ), y (t ), z (t )) ，藉由弧長關係式可寫成
2222
(ds ) =(dx ) +(dy ) +(dz ) (1)
若將時間參數表示法轉至空間參數表示法，則
r =(x (t ), y (t ), z (t )) =((s ), (s ), (s )) (2)
其中，r 為位置向量。單位切向量τ，定義為
dr
τ= (3)
ds
其中，ds 為微小弧長。單位法向量ν因與τ正交，
故可令
τ
ν= (4)
單位雙法向量(binormal vector)β，則與τ和ν向量正交。因此可得
β=τ⨯ν (5)
則τ、ν與β三向量的關係，如圖一所示
圖一 切向量τ 、法向量ν 與雙法向量β
由(3)式可知，單位切向量可表示為
τ(s ) =
(s ) r (s ) r
=
1
ρ
(17)
將(17)式代入(4)式，可知 又 則
ν=ρτ (18)
β⋅β=1
=0β⋅β
=p τ+q νβ
(19) (20)
(21)
用(20)式可知，β與β 正交關係，故可令
(6)
因
將
τ(s ) ⋅τ(s +ds ) =(s ) (s +ds ) cos(d θ) =cos(d θ) (7)
τ⋅β=0
利用泰勒展開式將τ(s +ds ) 與cos(d θ) 展開，可得
111⎡⎤ (s )(ds ) 2+τ (s )(ds ) 3+ τ(s ) ⋅τ(s +ds ) =τ(s ) ⋅⎢τ(s ) +τ (s )(ds ) +τ⎥1! 2! 3! ⎣⎦
∞
則
(18)式代入(23)式可得
(22) (23)
+τ ⋅β=0τ⋅β
=τ(s ) ⋅∑
n =0
τ(s )
n !
2
n
(ds ) (8)
n
+τ⋅β
cos(d θ) =1-
12! 12
(d θ) +
14!
(d θ) -
4
16! 1
ν⋅β=0 (24) ρ
1
6
(d θ) + (9)
因ν與β正交，可得
=0τ⋅β
將(8)與(9)式代入(7)式，可得
1+
(s )(ds ) =1-τ(s ) ⋅τ
2
(25)
2!
(d θ)
2
(10)
由(21)式等號兩邊與τ作內積，可得 因τ、ν正交，可得
⋅τ=p τ⋅τ+q ν⋅τ (26) β
p =0 (27) =q νβ
(s )(ds ) 2=-(d θ) 2 (11) τ(s ) ⋅τ
又 可推得
τ⋅τ =0 (12)
=-τ ⋅τ (13) τ⋅τ
所以 因此令
(28)
ν
(13)式代入(11)式可得
22
-τ ⋅τ (ds ) =-(d θ)
(14)
其中
=-β
1
則
2
σ
1 β
(29)
111⎛d θ⎫
= === (15) ⎪222
ρ⎝ds ⎭⎛ds ⎫⎛ρd θ⎫
⎪ ⎪⎝d θ⎭⎝d θ⎭
2
σ=
(30)
σ為扭率，而
ρ=
ds d θ
β=τ⨯ν (31)
(16)
又ν、β與τ相互正交，故
圖二 曲率半徑ρ
由(15)式與圖二可得
ν=β⨯τ
ν =
1
(32)
τρ
1
σ
β-
(33)
整理(18)、(29)與(33)式得狀態空間表示式如下：
⎡⎢0
⎧τ ⎫⎢⎪ ⎪⎢1⎨ν ⎬=⎢-
ρ ⎪⎪β⎢
⎩⎭⎢
⎢0⎣
1
ρ
0-1
σ
⎤0⎥⎥1⎥⎥σ⎥⎥0⎥⎦
2. 矩陣函數技巧與應用[5-17]
⎧τ⎪ ⎨ν ⎪β⎩
⎫⎪
⎬ (34) ⎪⎭
2.1 傳統方法
假設有一矩陣A ，其特徵向量組成的矩陣為V
V =[v 1⋅⋅⋅v n ] (39)
Av j =λj v j ，j =1, ⋅⋅⋅, n (40)
可寫成
AV =VD (41)
令
⎧τ⎪ P =⎨ν ⎪β
⎩
⎫⎪
⎬ (35) ⎪⎭
其中
⎡λ1⎢0
D =diag (λ1, ⋅⋅⋅, λn ) =⎢
⎢0⎢⎣0
00 0
則(34)式可寫成
λ2
00
即為 Frenet formula，其中
⎡⎢0⎢⎢1A =⎢-
ρ⎢⎢⎢0⎣
1
=AP (36) P
⎤
0⎥⎥1⎥
⎥ (37) σ⎥⎥0⎥⎦
0⎤⎥0
⎥ (42) 0⎥⎥λn ⎦0⎤⎥0⎥
(43) ⎥0⎥λt e n ⎦⎥
則
⎡e λ1t
⎢0λn t λ1t
=diag (e , ⋅⋅⋅, e ) =⎢
⎢0⎢⎢0⎣
e
A t
ρ
0-1
0e
λ2t
00 0
e
Dt
00
-1
σ
可看出A 為反對稱矩陣。 1.2 求解過程
計算步驟流程分述如下：
(a) 轉換參數表示法。將時間曲線參數式，轉換成弧 長表示式r (s ) 。
(b) 求單位切向向量τ。
(c) 求曲率半徑ρ。 (d) 求法向量ν。 (e) 求雙法向量β。
(f) 求扭率σ。
(g) 列出聯立微分方程並求解。
1⎧
=ντ⎪
ρ ⎪
11⎪
β⎨ν =-τ+
ρσ ⎪
⎪1
⎪β=-ν
σ ⎩
, τ(0)=τ0
, ν(0)=ν0 (38) , β(0)=β0
所以
=V e V
D t
(44)
2.2 Jordan 正則式(Canonical Form)
假設有一n ⨯n 矩陣A ，其特徵值有二重根，可利用
Jordan Canonical Form解決重根問題，矩陣特徵值有二重根時Jordan Form表示為
⎡λ1⎢0⎢⎢0⎢⎢0⎢0⎣
00 00
000
λ2
000
A =V
λi
0⎤⎥0⎥0⎥V ⎥1⎥λi ⎥⎦
-1
(45)
其中λi 為特徵值，以下為λ四重根
⎡λ
⎢0J =⎢
⎢0⎢⎣0
1
01
λ
00
λ
2
0⎤⎥0
⎥ (46) 1⎥⎥λ⎦
t
3t ⎤⎥3! ⎥2t ⎥⎥2! ⎥t ⎥⎥1⎦
所以
⎡⎢1⎢⎢⎢0⎢⎢0⎢⎣0
t 100
2! t 10
給定初位置(x (0),y (0),z (0))，與初架構(frame) (h) 利用 Frenet formula 求得的解，並將τ、ν與β座標架構隨參數變化以動畫呈現。
(i) 利用Mathematica 符號運算軟體作出三維軌跡動
畫。
(τ(0),ν(0),β(0))
e
Jt
=e
λt
(47)
推導過程如下：
喬登基本矩陣定義如下
⎡λ⎢0J =⎢
⎢0⎢⎣0
1
01
λ
00
λ
0⎤⎥0
⎥ (48) 1⎥⎥λ⎦
⎤1
f '''(λ) ⎥3! ⎥⎥1
f ''(λ) ⎥2! ⎥⎥1
f '(λ) ⎥1! ⎥
則
⎡
⎢f (λ) ⎢⎢⎢0=⎢⎢⎢0⎢⎢⎣0
⎡⎢1⎢⎢⎢⎢0⎢⎢0⎢⎢0⎣
1
f '(λ) 1! f (λ) 00t 100
t 22! t 10
1
f ''(λ) 2! 1
f '(λ) 1! f (λ) 0t 3⎤3! ⎥t 2⎥
⎥2! ⎥t ⎥
⎥⎥⎥
f (J ) =e Jt
f (λ)
⎥⎦
=e
λt
(49)
1⎥⎦
2.3 矩陣餘式定理[18]
給一矩陣A ，其特徵方程式為f (λ) =0，特徵
值為λ1⋅⋅⋅λn 。由Cayley-Hamilton 定理知f (A ) =0。透過餘式定理可知
e λ=f (λ) q (λ) +a n -1λn -1+a n -2λn -2+⋅⋅⋅
+a 1λ+a 0
4. 狀態方程求解
軌跡方程可寫成狀態方程，因此可用e A t 求解。
簡單列出其解題過程，並為接下來的曲線動畫做準備。
例題1. 解如下狀態方程，求切向量τ(s ) 、法向量
ν(s ) 與雙法向量β(s ) 。
(50)
將λ1⋅⋅⋅λn 分別代入(50)式，可求得a n ⋅⋅⋅a 0。再由實數和矩陣可互換性質，將實數換為矩陣可寫成
e At =f (A ) q (A ) +a n -1(t ) A n -1+a n -2(t ) A n -2+⋅⋅⋅
+a 1(t ) A +a 0I
(51)
代入矩陣A 得e 。
A t
已知控制方程
⎡⎢0⎢⎢1⎢-⎢2⎢⎢0⎣
研究結果
3. 比較二維各種曲率表示式
在 Frenet formula 中，扭率若為零的特例情況 下，可將空間曲線簡化為二維平面曲線。各種二維曲線之曲率計算方法整理比較於表一。
⎧
τ (s ) ⎫⎪⎪⎪ ⎪⎨ν (s ) ⎬⎪ ⎪ (s ) ⎪⎪β⎩⎭
1
20-12
=
⎥⎧τ(s ) ⎫⎪ ⎪1⎥⎪⎪ ⎥⎨ν(s ) ⎬2⎥⎪ ⎪
β(s ) ⎪⎥⎪⎩⎭0⎥ ⎦
0⎥
⎤
(52)
初架構(frame)
⎧
⎪τ(0)= ⎪⎪
⎨ν(0)=(1,0, 0) ⎪
⎪β(0)=(0,⎪⎩
(53)
初位置
⎧x (0)=1
⎪
⎨y (0)=0⎪
⎩z (0)=0
(54)
設
⎧τ(s ) ⎫⎪ ⎪
P =⎨ν(s ) ⎬ (55)
⎪β(s ) ⎪⎩⎭ ⎧τ (s ) ⎫ ⎪ =⎪νP ⎨ (s ) ⎬ (56) ⎪β (s ) ⎪
⎩⎭
將(54)式帶入，則空間曲線(x (s ), y (s ), z (s )) 點式表示如下
(x (s ), y (s ), z (s )) =(- (65)
而
5. 曲線動畫與Frenet formula參數研究
利用例題1. 的曲線針對初位置、初架構、曲率與扭率做參數研究。 5.1 曲線動畫
黑線為空間曲線，紅線表示τ(s ) ，綠線表示ν(s ) ，藍線表示β(s ) ，如圖三所示
則 可解得
=AP (57) P ⎧τ(0)⎫
⎪ ⎪
⋅⎨ν(0)⎬ (58) ⎪β(0)⎪⎩⎭
P =e
As
其中，
⎡-1⎢
=⎢-⎢1⎣
-11
⎡⎢1⎤⎢⎥
0⎥. ⎢0
⎢1⎥⎢0⎦⎢⎢⎣
e
As
e 0
⎤⎡-1
0⎥⎢4⎥⎢⎢10⎥. ⎢-⎥40s
e ⎥⎢⎥⎢1⎥⎦⎢⎣2
1⎤
4⎥⎥1⎥4⎥⎥1⎥⎥2⎦
(59)
代入(53)式
⎡⎢2
⎢cos ⎢
⎢⎢⋅C =⎢⎢⎢⎢⎢2
⎢sin ⎢
⎢⎣
⎤
⎥
2sin ⎥
⎥
⎡⎥0⎢
⎥
⎥⎢
. ⎢1⎥⎢⎥⎢⎥0⎢⎣⎥2
cos ⎥
⎥
⎥⎦
圖三 空間曲線圖
(60)
5.2 改變初位置
0e
As
0⎥⎥改變初位置，僅平移改變其初始點如下：
⎧x (0)=1⎧x (0)=2
⎪⎪
(a ) ⎨y (0)=0與(b ) ⎨y (0)=0 ⎪⎪⎩z (0)=0⎩z (0)=0
得
⎡
⎢⎢
⎧τ(s ) ⎫⎢
⎪ ⎪⎢⎨ν(s ) ⎬=⎢⎢ ⎪β(s ) ⎪
⎩⎭⎢ ⎢⎢
⎢⎣
-⎥0
⎥⎥ (61)
黑線為(a ) 與藍線為(b ) 。發現曲線形狀並未改變。
如圖四所示
再把τ(s ) 對s 做積分得
=-x (s ) =
⎰
+C 1
(62)
圖四 改變初位置表示圖
5.3 改變曲率半徑ρ
改變曲率半徑則會影響曲線投影於xy 平面上圓之大小，如圖五(a)與(b)所示
y (s ) =
⎰
s 0
ds =+C 2
(63)
z (s ) =
⎰
s 0
ds =
C 3 (64)
①ρ②ρ③ρ④ρ
=2
α⎫α⎫⎛⎛
x =a cos s cos() ⎪, y =a sin s cos() ⎪, z =s sin(α)
a ⎭a ⎭⎝⎝
(66)
=2.4 =2.6
=2.2
a 為半徑，α為角度，s 為變數，如圖七所示
⑤ρ=3.0 ⑥ρ=4.0
圖五(a) 改變曲率半徑之曲線(σ=2)
圖七 空間曲線之a 半徑與α角度
6.1 改變半徑
半徑控制曲線投影於xy 平面上形成圓半徑之大小， 半徑越大圓越大，如圖八(a)與(b)所示，其中固定(角度) α=
π
4
圖五(b) (σ=2)
(τ(s ), ν(s ), β(s ))
，紅線半徑為a =1，綠線半徑為a =2，
藍線半徑為a =4。
5.4 改變扭矩σ
改變σ扭矩常數為影響曲線之除峭程度，σ越
大曲線較為平緩，如圖六(a)與(b)所示
①σ②σ③σ④σ⑤σ⑥σ
=2 =2.2 =2.4 =2.6 =3.0
紅線a =1 綠線a =2 藍線a =4
圖八(a) 改變半徑之曲線(側視圖) (α=
π
4
)
=4.0
圖六(a) 改變σ之曲線(ρ=2)
圖八(b) 改變半徑之曲線(俯視圖) (α=
π
4
)
圖六 (ρ=2)
6.2 改變角度
角度控制曲線除峭程度，為一曲線之水平角，角度越大曲線越除，如圖九(a)與(b)所示，其中固定(半徑) a =1，紅線角度為α=藍線角度為α=
π
16
(τ(s ), ν(s ), β(s ))
π
4
，綠線角度為α=
π
8
，
6. 空間曲線之時間與弧長參數表示法[19]
曲線時間與弧長參數表示
。
,
所以
紅線α=綠線α=藍線α=
π
4π8
∂R
=r sin(ϕ) ∂θ
(71)
π
16
圖九(a) (a =1)
由於
d
e r =sin(ϕ) cos(θ) i +sin(ϕ) sin(θ) j +cos(ϕ) k (72)
e ϕ=cos(ϕ) cos(θ) i +cos(ϕ) sin(θ) j -sin(ϕ) k (73)
e θ=-sin(θ) i +cos(θ) j +o k (74)
d ϕd θ
(e r ) =cos(ϕ) cos(θ) i -sin(ϕ) sin(θ) i dt dt dt +cos(ϕ) sin(θ)
d θd ϕ
j +sin(ϕ) cos(θ) j -sin(ϕ) k dt dt dt
d d ϕd θ(e ϕ) =-sin(ϕ) cos(θ) i -cos(ϕ) sin(θ) i dt dt dt
j +cos(ϕ) cos(θ) j -cos(ϕ) k
dt dt dt
d d θd θ
(e θ) =-cos(θ) i -sin(θ) j +0k
dt dt dt
d ϕ
d θ
d ϕ
d ϕ
(75) (76)(77)
-sin(ϕ) sin(θ)
圖九(b) (a =1)
得
⎧e r ⎪ ⎪ ⎨e ϕ ⎪ e ⎪⎩ θ
⎫⎡0⎪⎢⎪
-ϕ⎬=⎢
⎪⎢-sin(ϕ) θ ⎪⎭⎣
ϕ0-cos(ϕ) θ
sin(ϕ) θ ⎤⎧e r
⎥⎪
cos(ϕ) θ ⎥⋅⎨e ϕ
⎥⎪e 0
⎦⎩ θ
⎫⎪
⎬ (78) ⎪⎭
7. 三維空間球座標狀態方程矩陣討論[20]
7.1 球座標
向量R
，長度r 如圖十所示
上式狀態方程矩陣亦是反對稱。 7.2 圓柱座標
空間上任一點(x , y , z ) 可以表示為圓柱座標
(r θ, z , ，如圖十一所示)
y
⎧x =r sin(ϕ) cos(θ)
⎪
⎨y =r sin(ϕ) sin(θ) ⎪z =r cos(ϕ) ⎩
x
圖十 球座標定義 則
(67) (68)
圖十一 圓柱座標定義
R =(x , y , z ) =(r sin(ϕ) cos(θ), r sin(ϕ) sin(θ), r cos(ϕ))
分別對r , θ, ϕ偏微
R =r sin(ϕ) cos(θ) i +r sin(ϕ) sin(θ) j +r cos(ϕ) k
∂R
=sin(ϕ) cos(θ) i +sin(ϕ) sin(θ) j +cos(ϕ) k ∂r
∂R
, =1 (69) ∂r
∂R
=r cos(ϕ) cos(θ) i +r cos(ϕ) sin(θ) j -r sin(ϕ) k ∂ϕ
R =(x , y , z ) =(r cos(θ), r sin(θ), z ) (79)
令位置座標R =(x , y , z ) =(r c o s θ(), r s i θn (), z ，取方)
向導微
∂R ∂R
=cos(θ) i +sin(θ) j +0k , =1 (80) ∂r ∂r ∂R ∂R
=-r sin(θ) i +r cos(θ) j +0k , =r ∂θ ∂θ
∂R
, =r ∂ϕ
(70)
則
(81)
∂R
=-r sin(ϕ) sin(θ) i +r sin(ϕ) cos(θ) j +0k ∂θ
∂R ∂R
=0i +0j +1k , =1 (82) ∂z ∂z
台北市, 文魁 (2002).
[4] 陳正宗, 工程數學(二) 講義, 海大河工力學聲響振
動實驗室 (2008).
[5] Arieh Iserles and Antonella Zanna, “Efficient
所以 computation of the matrix exponential by
generalized polar deco mpositions, ” SIAM Journal r =-sin(θ) θ i +cos(θ) θ j +0k (86) e
on Numerical Analysis, Vol. 42, pp. 2218-2256
θ=-cos(θ) θi -sin(θ) θj +0k (87) [6] J.T. Chen, S.R. Kuo and C.F. Lee, “A new point of e
view for the Householder matrix by using matrix
z =0i +0j +0k (88) e
exponential, ” International Journal of Applied
整理得 Mathematics , Vol.7, No.3, pp.289-308 (2001).
[7] Kenneth G . Wilson and Michael E. Fisher, “Critical r ⎫⎡0θ 0⎤⎧e r ⎫⎧e
Exponents in 3.99 Dimensions, ” Physical Review ⎥⎪ ⎪⎪ ⎪⎢
θ⎬=⎢-θ00⎥⋅⎨e θ⎬ (89) ⎨e Letters , Vol 28, No.4, pp.240-600 (1972). ⎪ ⎪⎪e ⎪⎢⎥[8] J.S. Kole, “Solving seismic wave propagation in 00⎩e z ⎭⎩ z ⎭⎣0⎦ elastic media using the matrix exponential
上式狀態方程矩陣亦是反對稱。
approach, ” W ave Motion, Vol.38, pp.279-293 (2003) . [9] C. Morler and C. V. Loan, “Nineteen Dubious Ways 結論與建議
to Compute the Exponential of a Matrix, ” SIAM Review , Vol.20, No.4, pp.801-836 (1978). 本文先討論e A t 的兩種計算方法(含傳統與矩陣
[10] M. Vidyasagar, “A Novel Method of Evaluating 餘式定理) ，並用於空間曲線中求Frenet formula的曲
in Closed Form,” IEEE Transactions on Automatic
線軌跡，同時也對Frenet formula中的曲率與扭率做
Control , Vol.1, pp.600-611 (1970).
更進一步的分析及探討，包含曲線動畫、改變曲率
[11] H.-K. Hong, H. S. Lan and J.T. Chen, “On Residue
半徑、改變扭率與改變初位置均一併討論。在求解Theorem of Matrix Function and Its Application to A t
e 矩陣函數時，傳統方法較為一般人所知，但當遇Spin Tensor Using Euler Angle Concept,” (1990).
[12] Thomas J.R. Hughes and James Winget, “Finite 到A 矩陣特徵值有重根時，就必須利用到Jordan 正
Rotation Effects in Numerical Integration of Rate 則式(Canonical Form) 的技巧，而矩陣餘式定理較具
Constitutive Equations Arising in Large 一般性，適合解任一之矩陣函數，而遇到重根問題
Deformation Analysis,” Short Comm (1979).
時可運用微分運算求得所缺少的方程式，是個實用
[13] R.Rubinstein and S.N. Atluri, “Objectivity of
面和適用性較高的方法。接下來透過純座標轉換，Incremental Constitutive Relations over Finite Time 展現Frenet formula獨有的特性，找到對應的反對稱Steps in Computational Finite Deformation 矩陣，最後建立空間曲線之時間與弧長參數表示Analysis, ” Comp. Meth. App. Mech. Engng., Vol.36, 法，對於二維曲線則視A 特例亦有探討。由於動畫pp.277-290 (1983).
[14] D.B. Metzger and R. N. Dubey, “Corotational 模擬的方式可將較為抽象的三維空間曲線隨弧長變
Rates in Constitutive Modeling of Elastic-Plastic 化完整呈現於眼前，引起學習之興趣，故本文不僅
Deformation,” Int. J. Plasticity, V ol.4, pp.341-368 可提供Mathematica 動畫之初學者作為參考，還可以
(1987).
當作教學上之一輔助工具，提升教與學之效能。
[15] Juan C. Simo and Jerrold E. Marsden, “On the
Rotated Stress Tensor and Material V ersion of the 誌謝
Doyle - Ericksen Formula” (1983).
[16] D.P. Flanagan and L.M. Taylor, An Accurate 本研究感謝教育部卓越教學計畫與教學中心特
Numerical Algorithm for Stress Integration with 別規劃「大學生暑期學習計畫」，及國立台灣海洋大
Finite Rotations, Comp. Meth. App. Meth. Eng., 學河海工程學系力學聲響振動實驗室之陳正宗終身
V ol.62, pp.305-320 (1987).
特聘教授以及實驗室學長們的辛勤指導與感謝戴志
[17] R.C.Ward, “Numerical Computation of the Matrix
豪同學在二維特例部分的討論。 Exponential with Accuracy Estimate, ” SIAM J.
Number Anal, Vol.14,No.4, pp.600-610 (1977).
參考文獻
[18] 陳正宗, 矩陣函數的新解法–矩陣餘式定理之
[1] J.M. Gere, Mechanics of Materials, 3rd Edition, 應用, 新新季刊, 中山科學研究院, 龍潭,V ol.15,
Prentice Hall (1997). No.4 (1987). [2] 張起豪,Mathematica 套裝程式使用與設計, 台北[19] R. Aris, V ectors Tensors and the Basic Equations
市, 碁峰資訊 (1994). of Fluid Mechanics, 協成書局 (1983). [3] 余家銘 編著,Mathematica 程式設計風格與應用, [20] 林琦焜, 向量分析, 滄海書局 (2007).
e r =cos(θ) i +sin(θ) j +0k (83)
e θ=-sin(θ) i +cos(θ) j +0k (84)
e z =0i +0j +1k (85)