基于MATLAB的空間曲線曲率的數(shù)值計算

孟 彬,韓燕苓

齊魯工業(yè)大學(山東省科學院) 數(shù)學與統(tǒng)計學院,濟南 250353

曲線的幾何特征可以由曲線的彎曲程度刻畫,曲線的彎曲程度稱為曲線的曲率[1]；三維歐氏空間中,曲線的幾何特征不僅可以用彎曲程度刻畫還可以用扭曲程度刻畫,空間曲線的扭曲程度稱為曲線的撓率[1]。曲率在軌道設計、橋梁設計等方面有著重要的應用,陳志剛等[2]介紹了利用曲率類屬性預測儲層裂縫的應用。撓率在生活中也起著不可或缺的作用,C Blankenburg等[3]給出了3D圖像中曲線撓率的估計。除曲線外,曲面幾何特征也可由曲面的曲率刻畫,曲面曲率在圖像處理和地形分析等方面有著重要的作用,文獻[4]至文獻[6]介紹了曲面曲率的求法及應用。

曲率是探索曲線幾何特征最重要的概念之一,本文的主要任務是計算平面、空間曲線的數(shù)值曲率,特別是離散點的數(shù)值曲率。

1 解析方法

對于平面(空間)上任何一條正則曲線C,都可以用顯函數(shù)y=f(x)或者隱函數(shù)F(x,y)=0來表示,顯函數(shù)和隱函數(shù)都可以表示成向量函數(shù)

r(t)=(x(t),y(t))或

r(t)=(x(t),y(t),z(t))。

(1)

定義[1]設曲線C的方程為r(s),其中s是曲線的弧長參數(shù)。

(2)

1.1 平面曲線

解析方法:當曲線(1)的參數(shù)t是弧長參數(shù)s時,曲率的經(jīng)典計算公式[1]為：

k(s)=|r''(s)|。

(3)

弧積分的被積函數(shù)是向量函數(shù)導函數(shù)的模長,若模長不是常數(shù),被積函數(shù)通常含有開方運算,弧積分計算困難。于是,在平面曲線中,(3)式只適用于計算向量函數(shù)模長為常數(shù)的曲率。

例1:計算r(t)=(3sint,3cost),t∈(0,2π)的曲率。

解:上述曲線的參數(shù)不是弧長參數(shù),首先將參數(shù)用弧長參數(shù)表示：

(4)

將(4)式代入原式

由公式(3)

例1的導函數(shù)模長是一個常數(shù),所以參數(shù)t可以表示為弧長參數(shù)s,但是一般情況下弧積分不能表示為關于t的次數(shù)為1的單項式,所以很難用弧長參數(shù)表示任意參數(shù),這使得計算曲率的解析解非常困難。

1.2 空間曲線

由于空間曲線存在法向量,本文首先引入空間曲線的兩個法向量[1]:

2)次法向量:曲線的切向量和主法向量唯一確定了曲線的次法向量γ(s)=α(s)×β(s)。

當曲線r(s)=(x(s),y(s),z(s))的參數(shù)是弧長參數(shù)時,可直接由公式(2)計算曲線的曲率。當k(s)≠0時,對兩個法向量及其切向量做若干次微分運算以及內(nèi)積運算(陳維桓,微分幾何[1])后,可推導出求解空間曲線曲率的計算公式。

解析方法:對空間曲線的任意參數(shù)t,曲線曲率的經(jīng)典計算公式[1]為

(5)

例2:計算曲線r(t)=(5sint,5cost,3t),t∈(0,2π)的曲率。

解:由公式(5)需要計算曲線的一階導函數(shù)、二階導函數(shù)及其向量積:

r'(t)=(5cost,-5sint,3),

r''(t)=(-5sint,-5cost,0),

r'(t)×r''(t)=(15cost,-15sint,-25),

在實際生活中,如果需要計算曲率的對象是一些離散點而不是函數(shù)表達式時,很難用公式(3)和(5)計算曲率的解析解。此外,曲率的計算無論是在橋梁設計還是軌道設計等方面都有很重要的應用,所以計算曲率很有必要。由本節(jié)內(nèi)容得知計算曲率的解析解有很多局限性,于是本文給出一種計算曲線曲率數(shù)值解的方法。

2 數(shù)值方法

定理1[1]:設α(s)是曲線C的單位切向量,s是弧長參數(shù),用Δθ表示α(s)與α(s+△s)之間的夾角,則

由定理1與圖1以及向量夾角與內(nèi)積的關系,本文給出一種曲線曲率的計算公式。

圖1 曲線方向向量與夾角

定理2:設曲線C是一條正則參數(shù)曲線,它的方程為r(t),其中t是任意參數(shù),s是弧長,Δs是弧長增量,若極限

(6)

存在,則曲線C在t0點的曲率為k(t0)。

證明:由定義及定理1

(7)

平移α(s)與α(s+△s)至同一起點(圖1),α(s)與α(s+△s)夾角的余弦值可以由α(s)與α(s+△s)的內(nèi)積表示

(8)

(9)

夾角的大小與參數(shù)無關,所以

(10)

將(10)式代入(7)式

公式(6)是由內(nèi)積與夾角余弦的關系得到的,雖然此式也可以用來計算曲率解析解,但是反三角函數(shù)不易計算具體值,所以此式最主要用于計算曲率數(shù)值解。

特別,若正則曲線是一段圓弧,則曲率可用圓弧所在圓的半徑表示。

推論:設曲線C是一段圓弧,它的方程為r(t),其中t是角度參數(shù),R是圓弧所在圓的半徑,此時曲率可表示為

(11)

證明:由圖2,當曲線是圓弧時Δt=Δθ,所以

圖2 參數(shù)增量與方向向量夾角的關系

化簡后取極限

(12)

設圓弧的向量函數(shù)為r(t)=(Rsint,Rcost),t∈(θ1,θ2),則|r'(t)|=R,所以(12)可寫為

當曲線C是一段圓弧時,計算曲率只需要求曲線C的導函數(shù),這大大降低了計算曲率的難度,并且(11)式既便于計算曲率解析解也便于計算曲率數(shù)值解。

定理3 設曲線C是一條正則參數(shù)曲線,方程為r(t),其中t是任意參數(shù),Δt是參數(shù)t的增量,Δs是弧長增量(t是弧長參數(shù)時Δt=Δs),令k-(t0)、k+(t0)分別為t0的左、右曲率,若左曲率

(13)

與右曲率

(14)

均存在,則

k(t0)=k-(t0)=k+(t0)。

在三維歐氏空間中,(10)式仍然表示曲線切向量的夾角,但圓弧曲線中Δt不再等于切向量的夾角Δθ(實際上等于兩切向量夾角的方向余弦),所以(11)式只適用于求平面圓弧曲線的曲率,而(6)式中切向量的夾角在空間中的計算方法不變,所以(6)式也適用于計算空間曲線的曲率。

3 MATLAB程序?qū)崿F(xiàn)

3.1 驗證公式

本小節(jié)通過計算連續(xù)函數(shù)離散化后點的數(shù)值曲率驗證曲率計算公式(6)和(11)在實際問題中的適用性。

例3:編程計算y=sinx的數(shù)值曲率。

由上表可以看出誤差的階數(shù)都是10-3,而步長取定為10-2,誤差小于步長。

表1中誤差是直接將MATLAB計算出的結果與解析解(精確到小數(shù)點后15位)比較后的誤差,但通過MATLAB計算的數(shù)值結果與精確解(精確到小數(shù)點后15位)的關系可以看出,tk處的數(shù)值解與tk+1處的精確解“幾乎相等”。

表1 y=sinx數(shù)值曲率

選用向后差商公式[7]計算數(shù)值導數(shù)并按照左曲率公式(13)編程得到tk的數(shù)值曲率的值恰為表2中數(shù)值曲率tk+1的值,數(shù)值曲率與誤差如下表所示。

表2 數(shù)值解與精確解

對比表2與表3,數(shù)值導數(shù)公式的選取不同會導致數(shù)值曲率的計算結果不同。

表3 數(shù)值曲率與誤差

例4:編程計算y=x3,x∈(0,3)的數(shù)值曲率。

表4表明,用向前差商公式求數(shù)值導數(shù)和右曲率公式(14)編程求數(shù)值曲率的結果不存在表2中數(shù)值解與解析解“錯位”的問題。同時再用向后差商公式求數(shù)值導數(shù)和左曲率(13)編程求數(shù)值曲率的結果與表4結果相同。

表4 數(shù)值解與精確解

在MATLAB中可以求符號函數(shù)的極限(文獻[8]至文獻[10]給出了用MATLAB計算正則參數(shù)曲線弧長、導數(shù)、極限以及曲率解析解的相關介紹及MATLAB程序)但無法求離散數(shù)據(jù)的極限,在取極限運算時往往以小步長近似極限,所以表2與表3的結果可看作小步長下的左曲率(13)與右曲率(14)。在極限運算下左曲率與右曲率是相等的,但在離散數(shù)據(jù)的MATLAB編程中數(shù)值導數(shù)或左、右曲率的選擇不同數(shù)值曲率的計算結果也不同。

例3和例4表明,數(shù)值導數(shù)選向前差商公式時配合右曲率編程計算數(shù)值曲率的結果與數(shù)值導數(shù)選向后差商公式時配合左曲率編程計算數(shù)值曲率的結果相同,并且計算結果的精度優(yōu)于另外兩種組合。表3與表4的結果表明,步長取Δt=0.01,除個別端點外誤差均小于Δt2,誤差可由步長控制。

例5:編程計算例1的數(shù)值曲率。

解:首先引入求數(shù)值導數(shù)的三點公式[11]:

由公式(11)可知,圓的曲率是半徑的倒數(shù),當半徑是3時,曲率解析解為1/3,取步長0.01,用三點公式求數(shù)值導數(shù),對公式(11)編程計算,得到結果如下表。

表5 例1的數(shù)值曲率

由上表可以看出,數(shù)值曲率的誤差階是10-6,而步長是10-2,誤差可由步長控制。改變步長對比計算結果發(fā)現(xiàn),當誤差達到10-9后很難再減小,原因在于MATLAB軟件本身含有存儲誤差以及計算時產(chǎn)生的截斷誤差。

通過例3、例4和例5驗證了公式(6)和公式(11)適用于計算離散點處的數(shù)值曲率。計算數(shù)值曲率時可以直接用離散公式(5),但由于公式(5)需要計算向量函數(shù)的二階導數(shù),這會增大數(shù)值曲率計算的誤差,而(6)式只需要計算向量函數(shù)的一階導數(shù),理論誤差要小于(5)式。

例6:分別用公式(5)和公式(6)編程計算曲線r(t)=(5sint,5cost,3t),t∈(0,2π)的數(shù)值曲率。

解:由例2可知正螺旋線的曲率是5/34,取步長0.01,用向前差商公式求數(shù)值導數(shù),編程得到結果見表6:

表6 數(shù)值曲率

由上表可以看出,除最后一個點外,公式(5)的誤差都是10-6階,而公式(6)的誤差都是10-7階,所以公式(6)的計算結果要好于公式(5)。除此之外,通過改變本例的步長發(fā)現(xiàn),隨著步長的減小誤差很難再有大幅度的減小,原因是存在截斷誤差。上表中,公式(5)與公式(6)最后一個數(shù)值曲率均為0,這是由于求數(shù)值導數(shù)時使用的是向前差商法,差商法只能求(n-1)個點的數(shù)值導數(shù)(n為離散點的個數(shù)),為使計算時維度相同,可令tn的數(shù)值導數(shù)為向后差商,這使得最后兩個數(shù)值導數(shù)的值是相同的,所以最后一個點的數(shù)值曲率為0。用三點公式計算最后一點的數(shù)值導數(shù)可避免最后一個數(shù)值曲率為0的情況。

由以上四個例子可知,公式(6)與公式(11)適用于計算離散點處的數(shù)值曲率,并且本文例題指出數(shù)值曲率的誤差精度不會低于步長的精度。

3.2 實際應用

曲率半徑與向心力的計算有直接的關系,所以曲率半徑在鐵路、高速等軌道設計中起著至關重要的作用。

例7:給定一組離散數(shù)據(jù),計算各點處的數(shù)值曲率。

表7 離散數(shù)據(jù)

解:根據(jù)公式(6)編程計算可得各點處的數(shù)值曲率如表8。

表8 離散點的數(shù)值曲率

離散點圖像與數(shù)值曲率圖像:

由圖3可以看出:沿著x軸正方向,曲線相鄰點的切線夾角隨著x趨于正無窮逐漸趨于0,與圖4數(shù)值曲率走勢相吻合。

圖3 離散點圖像

圖4 數(shù)值曲率圖像

4 總 結

平面曲線和空間曲線的數(shù)值曲率都可以用公式(6)求解,但是在計算數(shù)值導數(shù)時選用的數(shù)值算法不同會導致計算結果的精度不同。例題3至例題5表明,不管用向前(有的是向后)差商公式還是三點公式求數(shù)值導數(shù),得到結果的精度都不會低于步長的精度,所以數(shù)值曲率的精度可由步長控制。對于公式(6)而言,它需要計算反三角函數(shù),計算曲率的解析解沒有普適性；對于公式(11)而言,它既可以用來計算數(shù)值曲率也可以計算解析曲率。通過編程計算發(fā)現(xiàn):隨著步長的無限減小,誤差精度很難再有提升,原因在于MATLAB的存儲誤差以及在計算時產(chǎn)生的截斷誤差。