非光滑曲線等長樣條逼近?

何國良, 張 勇

(1-電子科技大學數學科學學院，成都 611731;2-西南醫科大學醫學與信息工程學院，瀘州 646000)

1 引言

在數值傳熱學、計算流體力學、核反應堆模擬、油藏數值模擬、計算機幾何學等學科中常常需要數值求解定義在曲線上的微分方程．在這其中一些精細分析和計算中，為達到比較好的計算效果，需要用光滑曲線來近似已知非光滑曲線，并且還要保持沿曲線的測度不變—保持定義在這上面的重要物理、化學或幾何性質不變，比如濃度、壓力、速度、通量、質量等[1-5]．由于這些問題通常對應著復雜的微分或積分方程，需要用數值方法來求解．此外，由于物理模型要求具備守恒性，這就要求所采取的數值方法也能夠保持這些守恒性質．這樣一來，針對微分或積分算子方程本身的離散方法除要具備守恒特點而外，對網格也提出了很高的要求．比如在生成貼體網格的過程中，需要網格能夠很好地近似曲線(面)隨空間自變量的變化而變化的趨勢[5,6]．對于定義在曲線上的一些方程，通常采用的參數化目標曲線來生成網格的辦法存在兩個方面的不足：

1) 參數化曲線后所導致的積、微分方程異常復雜，所以除簡單問題外，一般在工程上很少使用[7,8]；

2) 對于一些特殊問題，達不到需要的精度，如圖1所示．

圖1:曲線上的投影點不唯一

曲線C代表一條管狀結構(或來自一薄層的簡化，詳細介紹參考文獻[9–12])，在A處有一個不光滑點．實際應用中需要以該曲線為基線，生成貼近曲線走勢的一組三角形(四邊形)劃分：為保證物理守恒，要求這些三角形或四邊形單元節點的值等于相應函數在曲線C上的正交投影[8-10](這樣在數值計算時，能保證我們關心的物理量守恒，從而保證數值方法和物理模型一樣具有守恒性)．對于曲線的光滑部分，劃分沒有障礙，但是在曲線不光滑部分，比如A點，則會遇到下面的問題：

1) 在曲線不光滑的地方，因為曲線方程在這里不可導．對于角平分線上的任何一點pi來講，其在曲線C上的投影點是不確定的：Li,Ri都是在曲線上的投影點．所以，在圖1中的點A附近，如果直接進行三角形(或四邊形)劃分，其計算過程中是不穩定的；

2) 通常采取的另一種方法，“邊角光滑化(邊界部分以曲代直)”，在這里也變得不可行—為達到一定的計算精度，導致局部網格劃分異常密集，否則近似誤差較大，會污染鄰接節點的計算，甚至導致計算失敗．

為克服參數化和局部簡單的以曲代直的不足，本文提出在局部，用等長光滑曲線來代替原始曲線的方法．這里的曲線一方面滿足足夠的光滑性(如曲線關于自變量，比如x是可導的)，另一方面更為重要的是，保證所得曲線在替代范圍內(本文稱為為鏈接區間)和被代替的原始曲線具有相同的測度(如：長度)．

由于該問題是在實際計算中新遇到的需要解決的一個基本問題，所以在本文中將作詳細的分析討論．在下面的討論過程中，均假設曲線在包含不可導點的某一區間內有界且一維可測，如圖2中的左圖；但對于圖2的右圖的函數在x=0處附近震蕩或其他分形曲線等復雜無界曲線，則不在本文討論范圍內，因為曲線段的一維測度(長度)為無窮大．

注意，雖然實際是應用中可能不是直接使用長度的概念[6,10,18]，但為簡單將基本原理和方法敘述清楚，這里暫時拋開具體應用，僅以模型曲線來作為分析對象，并將抽象測度用長度來代替．

本文的組織結構如下：首先分析滿足上面所提及的光滑性和等測度的條件，然后建立包含一組數學關系和方程的數學模型；在此基礎上分析該模型的可解性，最后通過數值方法來求解．

圖2: 在x=0附近曲線長度可測和不可測

2 模型建立

假設給定兩點的坐標：p(xl,yl),q(xr,yr)，以及鏈接區間[xl,xr]內的曲線長度(或其他測度)Lm．根據光滑性和等長度的要求，需構造一連續函數S(x)，滿足如下三個條件：

1) 盡可能只改變原曲線不光滑點附近的曲線；

2) 鏈接區間內[xl,xr]光滑條件：S(x)至少是一階光滑函數；

3) 長度性一致條件：在區間[xl,xr]上的函數S(x)，其對應曲線的長度L需滿足指定的長度Lm，即

根據上面的基本要求，比較容易找到且計算穩定的函數是分段多項式函數[14-16]．由于要求該函數在區間[xl,xr]上至少具有一階連續導數(鏈接區間的端點用單側導數連續來表示)，且在兩端端點滿足相應函數值和導數值的插值條件．這樣，包含原來的不光滑點處的光滑性在內，至少有七個約束條件．基于這些基本分析，自然的選擇是利用分段3次樣條函數來近似．考慮到這里不需要處理太多的插值節點，并且還需要能比較簡單地計算出曲線的長度，達到和后續進一步插值計算或求解微分方程時空離散之間的簡單對接，所以這里采用一般樣條函數[15-18]．

設需求的函數具有如下形式

其中xm為鏈接區間內任意一點，a1,b1,c1,d1,a2,b2,c2,d2是待定參數．此外，為滿足長度相同的要求，這里假設ym也是未知的．

根據S(x)需要滿足的三個條件，我們可以得到如下的約束關系(?)：

1) 插值條件

2) 函數一階光滑條件

3) 曲線長度相等條件

其中s表示弧元，Lm為指定長度．

帶入S(x)的表達式，我們可以得到下面一組約束方程

和

這里k0,k1為原曲線在端點xl,xr處的導數值，ym為中間鏈接點處的函值，其值待定．由于約束方程(4)的引入，使該該問題變得比較復雜．下面分析求解情況．

設A為方程組(3)的系數矩陣，行列順序按照a1,b1,c1,d1,a2,b2,c2,d2的自然順序進行排列，記系數向量X=[a1,b1,c1,d1,a2,b2,c2,d2]T,B=[yl,yr,ym,k1,k2,0,0,0]T．對于方程組(3)，其中前三個為插值條件，只要xl,xm,xr三點互不相同，這三個方程則線性無關；中間兩個是不同端點xl,xr的導數約束；后面兩個為函數在鏈接點的一階光滑性約束．這8個方程構成的線性方程組兩兩線性無關，即A可逆．于是，在給定ym的情形下，我們能夠得到唯一的一組系數a1,b1,c1,d1,a2,b2,c2,d2，即S(x)和ym是一一對應的；并且下面特殊關系式成立

即，若S(x)為直線時，ym定在plpr的連線上(此時自然要求k1=k2)．為方便起見，形式上令

將(5)帶入(4)式，有

下面分析非線性方程(6)解的存在性情況．為方便起見，令

注1 函數f(xm,ym)中的兩個被積函數均隱含著xm,ym，所以是二元函數，其圖像如圖3所示．這里xm并不要求為原曲線的不光滑點．

圖3:f(x m,y m)的二維圖像和等值線圖

3 理論結果

為完整地解決該問題，我們首先從理論上討論方程(6)中解的存在性．下面分幾步來證明方程(6)的解：首先，證明在區間[?1,1]上，存在關于y軸對稱的屋頂函數和長度與給定曲線相同的三次一階光滑函數；其次，引入一個線性變換將一般的區間[c,d]映射至[?1,1]區間；最后，借助于介值定理來完成方程(6)的解的存在性證明．

命題1 設δ∈R，且δ>0．對給定的對稱區間[?δ,δ]，以及在區間端左端點x=?δ處的斜率k0>0，在x= δ處的斜率k1= ?k0．存數ym=和一個三次多項式P(x)，使得該函數同時滿足下面四條：

1) 插值條件

2) P(x)∈ C1[?δ,δ]，其中C1[?δ,δ]表示閉區間[?δ,δ]的一階連續光滑函數；

3) 一階導數單調下降：滿足k0>P′(x)>?k0,x∈[?δ,δ]；

4) 函數P(x)在區間[?δ,δ]凸性不變．

證明 為方便起見，假設ym>0，對應的函數形式上為區間[?δ,δ]三次多次多項式

滿足下面的插值條件

由于上述關于系數[a,b,c,d]的線性方程組兩兩之間線性無關，所以存在唯一解

在區間(?δ,0)中，函數P1(x)的一階，二階導函數分別為

且x∈(?δ,0)，即函數P1(x)在區間(?δ,0)向上凸，所以在該區間(?δ,0)上單調下降，滿足：0

另一方面，利用對稱性可以證明在區間(0,δ)內，若在右端點x=δ處，曲線的斜率k1=k0，在x=0處的斜率等于0．同樣存在一個三次多項式P2(x)函數，滿足下面插值條件

該函數且在區間(0,δ)向上凸，且一階導函數單調下降，在x=δ取得最小?k0．

將上述兩方面結合起來，令

則通過上面構造和分析，我們得到P(x)在區間(?δ,δ)連續的，滿足插值條件1)．由于

所以P(x)在區間(?δ,δ)上具有連續的一階導數，即P(x) ∈ C1[?δ,δ]．命題1中的第2)條和第4)條，已經在證明過程中得到證明．

注2 命題1在k0>0的條件下得到．對于k0<0，我們同樣也可以得到類似的結論，只需要將命題1第3)條中的導數變成單調增加即可．該命題的關鍵在于能找到一個滿足條件的函數即可．

命題2 在區間(?δ,δ)上存在一個具有一階連續光滑的三次多項式函數P(x)∈P3(x)，使得

證明 為方便起見，下面以k0>0為例來證明，對于k0<0，結論一致．

令ym=k0/3，根據命題1，存在一個具有一階連續光滑的函數P(x)．該函數在區間[?δ,0]單調上升，一階導數不變號，且0

同樣，在區間[0,δ]上

將上面兩個積分加起來，有

注3 上述命題的含義是：如圖4所示，在區間[?1,1]上，多項式P(x)所對應的曲線(x,P(x))的長度小于線段AB和BC的長度之和，即

圖4:P(x)的長度小于三角形的兩直角邊的長度和

命題3 對于定義在區間[a,b]上任意連續且具有可測長度的函數f(x)，如果f(x)在子區間(c,d)∈[a,b]以外不存在不可導點，則存在一個如下形式的函數

使得h(s)在區間上和f(x)在區間[a,b]具有相同的長度，并且h(s)在點x=δ和x=?δ處導數連續，且h′(?δ)= ?h′(δ)=k0?=0, ˉa< ?δ< δ< ˉb．

證明 取實數xl,xr滿足a

由于在區間[a,xl]，若LL=e?a，則令?a=a+ε0，可使得LL>xl??a．為簡單起見，以后均假設LL>xl?a,Lm>xr?xl,LR>b?xr．

記S(x,f(x))為曲線上任意一點，引入下面的線性變換：T:(x,f(x))7→(s,ˉy),

其中?S表示變換后曲線上的點，

則對于定義區間[xl,xr]上的點集S(x,f(x))，在變換(11)下的像集?S(s,ˉy)成為定義在區間[?δ,δ]上的連續曲線，曲線長度依然為Lm．

令

則上面定義的hm(s)在區間[?δ,δ]上連續，對應的曲線長度依然為Lm，且在左、右端點處的單側導數存在，互為相反數：由于f(x)不光滑點附件的曲線長度有限，所以總可以適當移動區間[xl,xr]，使得k0?=0．

設aˉ

從方程組(13)的前兩個方程解出b1,c1，有b1=k0+2δa1,c1=k0+2δa1? 3a1δ2．帶入第三個方程得到關于參數ˉa,a1的方程．由于第三個方程左邊關于這兩個參數均是連續函數，左邊積分值的范圍介于?(δ+ˉa)到無窮大之間．所以總可以找到恰當的值，使得在[ˉa,?δ]上，函數p1(s)所對應的曲線長度為LL，并且p1(x)=k0．

類似地，在區間[δ,b]上，存在函數p2(s)=a2s2+b2s+c2，使得在該區間上，函數p2(s)所對應的曲線長度為LR，并且p2(x)=k0．

最后，令

即為所求函數．

說明：1) 根據h(s)的構造方法，對任意的故對其作一個反向旋轉變化

得到原來區間[xl,xr]上的函數，在該區間以外也是光滑的．為避免引入太多記號，這里可將其記為h(x)．

2) 該命題的主要意義在于：從理論上通過平移和旋轉變換，構造一個簡單的屋頂函數h(x)：這個函數只在x=[xl+xr]/2處不可導(不光滑點)，而在其它地方是光滑的，并且在相應區間上保持曲線長度不變．構造這樣的函數主要目的是為以后用介值定理做準備，實際的數值計算中不需要構造該函數．

定理1 設f(x)為區間[a,b]上的連續且長度可測函數，其長度為L0，則存在滿足約束關系(?)的解．

證明 若f(x)在整個區間[a,b]上光滑連續，則自然滿足要求．下面證明不光滑的情況．

由于f(x)為區間[a,b]上的連續函數，我們可以找到一個區間[xl,xr]，滿足[c,d]∈[xl,xr]∈[a,b]．利用命題3，存在一個連續函數h(s)，使得h(s)除x=0以外均是光滑函數．

對于區間[?δ,δ]以外的部分，h(s)所對應的曲線長度和f(x)的在區間[xl,xr]以外相同，所以這里只需再證明，可用等長光滑函數來代替h(s)在區間[?δ,δ]上的部分即可．設曲線(c,f(x))在區間[xl,xr]上的長度為Lm．根據h(s)的構造，h(s)在區間[ˉa,ˉb]上光滑，在區間[?δ,δ]上的長度和f(x)在[xl,xr]的長度相同，即

若Lm=2δ，此時需k0=0,f′(x)≡c0，即f(x)為直線段，方程(6)有唯一的解

若Lm>2δ，f(x)不恒為常數，k0=h′(s)|s=?δ?=0．由于形如

中的被積函數關于變量s是連續函數，從而方程(6)中的兩項定積分值均存在的．令

其中

從式(7)可知，向量X的各分量關于ym均為非線性函數，且h1(s,ym),h2(s,ym)是關于ym連續函數，所以(15)式是關于ym的連續非線性函數．

因k0?=0，由命題2，存在數?y=k0/3和一個具有一階連續光滑的三次多項式函數P(s)∈P3(s)，使

另一方面，對于任意的固定值Lm，當

不恒為零時，H(ym)是關于ym的增函數(0

于是存在一正數M滿足M>Lm，如圖5所示．對該M，根據函數H(ym)的連續性，存在ˉym，使得H(ˉym)=M，即存在ˉym，使得

將式(17),(18)結合起來，根據介質定理，存在ym?，使得H(ym?)=Lm．再將其在區間[a,b]上進行適當延拓即可從而證明非線性方程(6)有解．從而定理得證．

圖5:不同的y m對應的光滑曲線

說明：1) 定理說明對于任意一個連接區間[xl,xr]，只要在這個區間上，曲線長度可測，則總可以用一條一階光滑曲線來進行等長替代．

2) 方程(6)的解一般不唯一．因為對于給定的xm和Lm，由于積分關系中凡是含ym的項整體均為平方項，所以若ym滿足方程(6)，則在?ym附近也存在另一個滿足方程(6)的解，如圖6所示．

3) 由于ym不唯一，實際應用中選擇ym的一般方法是借助于圖形的慣性方向來選取；比較簡單的方法是借助重心坐標來算，方法如下：

?計算曲線在區間[xl,xr]的部分，連接點(xl,f(x1))和點(xr,f(xr))；對給定的xm計

算兩直線

?在區間[xl,xr]內的交點，選ym的初值為所得三角形的重心即可，如圖6中的ym．

圖6: 方程(6)在圖中有上下兩個解

在上面的定理中，我們證明了在含不光滑點的鏈接區間[xl,xr]上存在通過該區間中點xm，且在中點處的導數為0．實際上由關系式(9)可知所以我們同樣可以構造出h(s)，使得從而保證h(s)在區間(xl,xr)內二階光滑．這樣一來，由h(s)衍生的S(x)同樣滿足此外，由于在插值節點xm處，樣條函數S(x),F(xm),H(ym)均是關于xm的連續函數．沿用定理1的證明方法，我們可得到下面更一般的結論．

定理2 設f(x)為區間[xl,xr]上的連續且長度可測函數，有意義，則對任意的x∈[xl,xr]，存在滿足如下約束關系函數三次樣條函數S(x)：

1) 插值條件

2) 函數光滑條件

3) 曲線長度相等條件

其中s表示弧元，Lm為原曲線在區間[xl,xr]上的長度．

注4 定理2表明，在一個局部區間，不僅可以用簡單曲線來不光滑曲線，也可以代替復雜曲線，比如高頻震蕩曲線．這在實際應用中是很有好處的，尤其在對曲線形狀和長度比較敏感的場合．

4 算法分析

下面討論光滑曲線S(x)的具體數值求解辦法．在鏈接區間上[xl,xr]，曲線長度Lm可利用測量數據或用公式積分或數值積分計算而得到．

由于一般形式積分關系(15)的結果，形式非常復雜，難于計算，所以這里考慮用高斯積分數值積分方法來計算．下面給出求解恰當y?m的牛頓-拉夫遜算法(New ton-Raphson Method)．

4.1 算法

1) 給定初始數據：xl,xr,f(xl),f(xr),k1,k2，及曲線的在區間[xl,xr]上的長度Lm，容許誤差ε；

2) 計算幾何區域的慣性方向(或重心)，選取恰當的xm∈[xl,xr]，計算ym相應的初值

3) 計算對應的S(x)：

令 n=1,2,3,···,

(I) 用高斯積分方法計算積分式(15)；

(II) 用牛頓法計算方程(6)的第n此迭代值；

(III) 利用算樣條函數的系數X=[a1,b1,c1,d1,a2,b2,c2,d2]T，從而得到新的?S(x)；

(IV)計算樣條曲線?S(x)在夾在區間[xl,xr]的長度L，并計算|L?Lm|；

(V) 若|L?Lm|>ε，則返回3)繼續計算；否則，終止跳出循環，結束計算；

4) 利用最后的ynm計算曲線整個區間上S(x)，從而得需要的光滑曲線．

4.2 算法分析

1) 雖然牛頓法對初值比較敏感，但從ym的角度上看，如果方程 (6)基本上還是和2次多項式方程接近，所以會很快收斂．只不過由于F(ym)存在最小值點，初值取為最小值點時，算法可能不穩定．所以選取初值時，需注意盡可能避開該點．

2) 理論上講，對任意的xm∈[xl,xr]都可以求出ym對應的近似值，不過從后續計算的穩定性和效率來看，在|k1|和|k2|比較接近的情況下，將xm選擇靠近區間的中點還是不錯的方案，這樣得到的曲線，極性小些[14,19]．

4.3 數值結果

下面以函數f(x)=|sin(x)|,x∈[?1.5,1.5]為例來計算．該函數在x=0處導數不存在，有一個不光滑點.在區間[?1,1]上面長度有限．

圖7為選擇不同的插值節點xm來計算的光滑曲線．他們分別對應xm分別為?2/3,0,2/3三個值的三對曲線．上面一組點線和下一組點線相對應，粗線為原始非光滑曲線；小圓圈表示插值節點xm在曲線上的位置．每一組曲線在長度、光滑程度均滿足要求，都可以作為局部的曲線近似．所以我們很容易根據實際需要來選擇相應的曲線，這保證了方法的穩定可靠．

圖7:不同x m對應的連續光滑曲線

表1的結果在牛頓迭代終止條件為10?5下取得．迭代次數5～6次．曲線在區間[?1,1]上的真實長度：2.622884996431094．計算結果僅保留小數點后面5位．

從表1可以看出，實際計算速度很快；而且對不同的xm，計算結果正確且算法穩定，所以本文提出的等長度曲線近似方法計算效果很好．為簡單得到唯一的曲線，只需將插值節點xm選為區間的中點即可．此外，本文主要討論了連續曲線的逼近問題，如果工程中需要對于間斷曲線進行近似，可以有兩個辦法：一個是逐短光滑，另一種是整體光滑．這兩種問題利用本文的處理方法不存在任何困難，因為本文提出的近似方法只需要知道在給定區間上曲線的長度即可，與曲線的具體幾何形態無關．

表1:不同x m對應的連續光滑曲線

參考文獻：

[1]mehendale S S,Jacobi A M,Shah R K.Fluid fl ow and heat transfer at micro-and meso-scales with app lication to heat ex-changer design[J].App lied mechanics Review s,2000,53(7):175-193

[2]吳向紅，葉繼根，馬遠樂，等.水平井蒸汽輔助重力驅油藏模擬方法[J].計算物理，2002,11(6):549-592 W u X H,Ye J G,Ma Y L,et al.SAGD num erical Simulation with horizontal wells[J].Chinese Jou rnal of Com putational Physics,2002,11(6):549-592

[3]Fang C,Steinbrenner J E,Wang F M,et al.Im pact of wall hyd rophobicity on condensation fl ow and heat transfer in silicon micro-channels[J].Jou rnal of microm echanics and microengineering,2010,20(4):045018

[4]Louah lia-Gualous H,mecheri B.Unsteady steam condensation fl ow patterns inside a miniatu re tube[J].App lied Therm al Engineering,2007,27(8-9):1225-1235

[5]陶文銓.數值傳熱學[M].西安：西安交通大學出版社，2001 Tao W Q.Num erical Heat Transfer[M].X i’an:X i’an Jiaotong University Press,2001

[6]范曉光，馬學虎.微通道內蒸汽及混合蒸氣冷凝[D].大連：大連理工大學，2012 Fan X G,Ma X H.Fluid fl ow and heat transfer of steam and steam-noncondensab le gasm ixtu res condensation in microchannels[D].Dalian:Dalian University of Technology,2012

[7]Floater M,Horm ann K.Su rface param eterization:a tu torial and su rvey[C]//Advances in mu ltiresolution for Geom etric modelling,Sp ringer,2005

[8]Greer JB.An im provem ent of a recent Eu lerian method for solving PDEs on general geom etries[J].Journal of Scientifi c Com puting,2006,29(3):321-352

[9]Salom on D.Cu rves and Surfaces for Com pu ter G raphics[M].New York:Sp ringer Science+Business med ia,2006

[10]Reu ter M,W olter F E,Shenton M,et al.Lap lace-Beltram i eigenvalues and topological featu res of eigenfunctions for statistical shape analysis[J].Com pu ter-A ided Design,2009,41(10):739-755

[11]Bertalm M,Bertozzi A L,Sapiro G.Navier-Stokes,fluid dynam ics,and im age and video inpainting[J].IEEE Com pu ter Society Con ference on Com pu ter V ision&Pattern Recognition,2001,1(1):355-362

[12]Dziuk G,E lliott C M.Eu lerian finite elem entm ethod for parabolic PDEs on im plicit su rfaces[J].Interfaces&Free Boundaries,2008,10(1):119-138

[13]macdonald C B,Ruu th S J.The implicit closest point method for the num erical solution of partial differential equations on su rfaces[J].SIAM Jou rnal on Scientific Computing,2009,31(6):4330-4350

[14]Burden R L,Fairs J D.Num erical Analysis[M].Boston:Cengage Learning Press,2001

[15]Berrut J P,Trefethen L N.Barycentric Lagrange interpolation[J].SIAM Review,2004,46(3):501-517

[16]Liu X D,Osher S,Chan T.W eighted essentially non-oscillatory schem es[J].Journal of Com putational Physics,1994,115(1):200-212

[17]Xu Z F,Shu C W.Anti-d iff usive fl ux corrections for high order finite d iff erenceW ENO schem es[J].Journal of Com putational Physics,2005,205(2):458-485

[18]李慶揚，王能超，易大義.數值分析[M].武漢：華中科技大學出版社，2001 LiQ Y,W ang N C,Y iD Y.Num erical Analysis[M].W uhan:Huazhong University of Science and Technology Press,2001

[19]鐘爾杰，黃廷祝.數值分析[M].北京：高等教育出版社，2004 Zhong E J,Huang T Z.Num erical Analysis[M].Beijing:H igher Education Press,2004