非光滑曲線等長(zhǎng)樣條逼近?

何國(guó)良, 張 勇

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

1 引言

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

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

2) 對(duì)于一些特殊問(wèn)題，達(dá)不到需要的精度，如圖1所示．

圖1:曲線上的投影點(diǎn)不唯一

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

1) 在曲線不光滑的地方，因?yàn)榍€方程在這里不可導(dǎo)．對(duì)于角平分線上的任何一點(diǎn)pi來(lái)講，其在曲線C上的投影點(diǎn)是不確定的：Li,Ri都是在曲線上的投影點(diǎn)．所以，在圖1中的點(diǎn)A附近，如果直接進(jìn)行三角形(或四邊形)劃分，其計(jì)算過(guò)程中是不穩(wěn)定的；

2) 通常采取的另一種方法，“邊角光滑化(邊界部分以曲代直)”，在這里也變得不可行—為達(dá)到一定的計(jì)算精度，導(dǎo)致局部網(wǎng)格劃分異常密集，否則近似誤差較大，會(huì)污染鄰接節(jié)點(diǎn)的計(jì)算，甚至導(dǎo)致計(jì)算失?。?/p>

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

由于該問(wèn)題是在實(shí)際計(jì)算中新遇到的需要解決的一個(gè)基本問(wèn)題，所以在本文中將作詳細(xì)的分析討論．在下面的討論過(guò)程中，均假設(shè)曲線在包含不可導(dǎo)點(diǎn)的某一區(qū)間內(nèi)有界且一維可測(cè)，如圖2中的左圖；但對(duì)于圖2的右圖的函數(shù)在x=0處附近震蕩或其他分形曲線等復(fù)雜無(wú)界曲線，則不在本文討論范圍內(nèi)，因?yàn)榍€段的一維測(cè)度(長(zhǎng)度)為無(wú)窮大．

注意，雖然實(shí)際是應(yīng)用中可能不是直接使用長(zhǎng)度的概念[6,10,18]，但為簡(jiǎn)單將基本原理和方法敘述清楚，這里暫時(shí)拋開(kāi)具體應(yīng)用，僅以模型曲線來(lái)作為分析對(duì)象，并將抽象測(cè)度用長(zhǎng)度來(lái)代替．

本文的組織結(jié)構(gòu)如下：首先分析滿足上面所提及的光滑性和等測(cè)度的條件，然后建立包含一組數(shù)學(xué)關(guān)系和方程的數(shù)學(xué)模型；在此基礎(chǔ)上分析該模型的可解性，最后通過(guò)數(shù)值方法來(lái)求解．

圖2: 在x=0附近曲線長(zhǎng)度可測(cè)和不可測(cè)

2 模型建立

假設(shè)給定兩點(diǎn)的坐標(biāo)：p(xl,yl),q(xr,yr)，以及鏈接區(qū)間[xl,xr]內(nèi)的曲線長(zhǎng)度(或其他測(cè)度)Lm．根據(jù)光滑性和等長(zhǎng)度的要求，需構(gòu)造一連續(xù)函數(shù)S(x)，滿足如下三個(gè)條件：

1) 盡可能只改變?cè)€不光滑點(diǎn)附近的曲線；

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

3) 長(zhǎng)度性一致條件：在區(qū)間[xl,xr]上的函數(shù)S(x)，其對(duì)應(yīng)曲線的長(zhǎng)度L需滿足指定的長(zhǎng)度Lm，即

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

設(shè)需求的函數(shù)具有如下形式

其中xm為鏈接區(qū)間內(nèi)任意一點(diǎn)，a1,b1,c1,d1,a2,b2,c2,d2是待定參數(shù)．此外，為滿足長(zhǎng)度相同的要求，這里假設(shè)ym也是未知的．

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

1) 插值條件

2) 函數(shù)一階光滑條件

3) 曲線長(zhǎng)度相等條件

其中s表示弧元，Lm為指定長(zhǎng)度．

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

和

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

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

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

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

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

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

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

3 理論結(jié)果

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

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

1) 插值條件

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

3) 一階導(dǎo)數(shù)單調(diào)下降：滿足k0>P′(x)>?k0,x∈[?δ,δ]；

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

證明 為方便起見(jiàn)，假設(shè)ym>0，對(duì)應(yīng)的函數(shù)形式上為區(qū)間[?δ,δ]三次多次多項(xiàng)式

滿足下面的插值條件

由于上述關(guān)于系數(shù)[a,b,c,d]的線性方程組兩兩之間線性無(wú)關(guān)，所以存在唯一解

在區(qū)間(?δ,0)中，函數(shù)P1(x)的一階，二階導(dǎo)函數(shù)分別為

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

另一方面，利用對(duì)稱性可以證明在區(qū)間(0,δ)內(nèi)，若在右端點(diǎn)x=δ處，曲線的斜率k1=k0，在x=0處的斜率等于0．同樣存在一個(gè)三次多項(xiàng)式P2(x)函數(shù)，滿足下面插值條件

該函數(shù)且在區(qū)間(0,δ)向上凸，且一階導(dǎo)函數(shù)單調(diào)下降，在x=δ取得最小?k0．

將上述兩方面結(jié)合起來(lái)，令

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

所以P(x)在區(qū)間(?δ,δ)上具有連續(xù)的一階導(dǎo)數(shù)，即P(x) ∈ C1[?δ,δ]．命題1中的第2)條和第4)條，已經(jīng)在證明過(guò)程中得到證明．

注2 命題1在k0>0的條件下得到．對(duì)于k0<0，我們同樣也可以得到類(lèi)似的結(jié)論，只需要將命題1第3)條中的導(dǎo)數(shù)變成單調(diào)增加即可．該命題的關(guān)鍵在于能找到一個(gè)滿足條件的函數(shù)即可．

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

證明 為方便起見(jiàn)，下面以k0>0為例來(lái)證明，對(duì)于k0<0，結(jié)論一致．

令ym=k0/3，根據(jù)命題1，存在一個(gè)具有一階連續(xù)光滑的函數(shù)P(x)．該函數(shù)在區(qū)間[?δ,0]單調(diào)上升，一階導(dǎo)數(shù)不變號(hào)，且0

同樣，在區(qū)間[0,δ]上

將上面兩個(gè)積分加起來(lái)，有

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

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

命題3 對(duì)于定義在區(qū)間[a,b]上任意連續(xù)且具有可測(cè)長(zhǎng)度的函數(shù)f(x)，如果f(x)在子區(qū)間(c,d)∈[a,b]以外不存在不可導(dǎo)點(diǎn)，則存在一個(gè)如下形式的函數(shù)

使得h(s)在區(qū)間上和f(x)在區(qū)間[a,b]具有相同的長(zhǎng)度，并且h(s)在點(diǎn)x=δ和x=?δ處導(dǎo)數(shù)連續(xù)，且h′(?δ)= ?h′(δ)=k0?=0, ˉa< ?δ< δ< ˉb．

證明 取實(shí)數(shù)xl,xr滿足a

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

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

其中?S表示變換后曲線上的點(diǎn)，

則對(duì)于定義區(qū)間[xl,xr]上的點(diǎn)集S(x,f(x))，在變換(11)下的像集?S(s,ˉy)成為定義在區(qū)間[?δ,δ]上的連續(xù)曲線，曲線長(zhǎng)度依然為L(zhǎng)m．

令

則上面定義的hm(s)在區(qū)間[?δ,δ]上連續(xù)，對(duì)應(yīng)的曲線長(zhǎng)度依然為L(zhǎng)m，且在左、右端點(diǎn)處的單側(cè)導(dǎo)數(shù)存在，互為相反數(shù)：由于f(x)不光滑點(diǎn)附件的曲線長(zhǎng)度有限，所以總可以適當(dāng)移動(dòng)區(qū)間[xl,xr]，使得k0?=0．

設(shè)aˉ

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

類(lèi)似地，在區(qū)間[δ,b]上，存在函數(shù)p2(s)=a2s2+b2s+c2，使得在該區(qū)間上，函數(shù)p2(s)所對(duì)應(yīng)的曲線長(zhǎng)度為L(zhǎng)R，并且p2(x)=k0．

最后，令

即為所求函數(shù)．

說(shuō)明：1) 根據(jù)h(s)的構(gòu)造方法，對(duì)任意的故對(duì)其作一個(gè)反向旋轉(zhuǎn)變化

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

2) 該命題的主要意義在于：從理論上通過(guò)平移和旋轉(zhuǎn)變換，構(gòu)造一個(gè)簡(jiǎn)單的屋頂函數(shù)h(x)：這個(gè)函數(shù)只在x=[xl+xr]/2處不可導(dǎo)(不光滑點(diǎn))，而在其它地方是光滑的，并且在相應(yīng)區(qū)間上保持曲線長(zhǎng)度不變．構(gòu)造這樣的函數(shù)主要目的是為以后用介值定理做準(zhǔn)備，實(shí)際的數(shù)值計(jì)算中不需要構(gòu)造該函數(shù)．

定理1 設(shè)f(x)為區(qū)間[a,b]上的連續(xù)且長(zhǎng)度可測(cè)函數(shù)，其長(zhǎng)度為L(zhǎng)0，則存在滿足約束關(guān)系(?)的解．

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

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

對(duì)于區(qū)間[?δ,δ]以外的部分，h(s)所對(duì)應(yīng)的曲線長(zhǎng)度和f(x)的在區(qū)間[xl,xr]以外相同，所以這里只需再證明，可用等長(zhǎng)光滑函數(shù)來(lái)代替h(s)在區(qū)間[?δ,δ]上的部分即可．設(shè)曲線(c,f(x))在區(qū)間[xl,xr]上的長(zhǎng)度為L(zhǎng)m．根據(jù)h(s)的構(gòu)造，h(s)在區(qū)間[ˉa,ˉb]上光滑，在區(qū)間[?δ,δ]上的長(zhǎng)度和f(x)在[xl,xr]的長(zhǎng)度相同，即

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

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

中的被積函數(shù)關(guān)于變量s是連續(xù)函數(shù)，從而方程(6)中的兩項(xiàng)定積分值均存在的．令

其中

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

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

另一方面，對(duì)于任意的固定值Lm，當(dāng)

不恒為零時(shí)，H(ym)是關(guān)于ym的增函數(shù)(0

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

將式(17),(18)結(jié)合起來(lái)，根據(jù)介質(zhì)定理，存在ym?，使得H(ym?)=Lm．再將其在區(qū)間[a,b]上進(jìn)行適當(dāng)延拓即可從而證明非線性方程(6)有解．從而定理得證．

圖5:不同的y m對(duì)應(yīng)的光滑曲線

說(shuō)明：1) 定理說(shuō)明對(duì)于任意一個(gè)連接區(qū)間[xl,xr]，只要在這個(gè)區(qū)間上，曲線長(zhǎng)度可測(cè)，則總可以用一條一階光滑曲線來(lái)進(jìn)行等長(zhǎng)替代．

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

3) 由于ym不唯一，實(shí)際應(yīng)用中選擇ym的一般方法是借助于圖形的慣性方向來(lái)選取；比較簡(jiǎn)單的方法是借助重心坐標(biāo)來(lái)算，方法如下：

?計(jì)算曲線在區(qū)間[xl,xr]的部分，連接點(diǎn)(xl,f(x1))和點(diǎn)(xr,f(xr))；對(duì)給定的xm計(jì)

算兩直線

?在區(qū)間[xl,xr]內(nèi)的交點(diǎn)，選ym的初值為所得三角形的重心即可，如圖6中的ym．

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

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

定理2 設(shè)f(x)為區(qū)間[xl,xr]上的連續(xù)且長(zhǎng)度可測(cè)函數(shù)，有意義，則對(duì)任意的x∈[xl,xr]，存在滿足如下約束關(guān)系函數(shù)三次樣條函數(shù)S(x)：

1) 插值條件

2) 函數(shù)光滑條件

3) 曲線長(zhǎng)度相等條件

其中s表示弧元，Lm為原曲線在區(qū)間[xl,xr]上的長(zhǎng)度．

注4 定理2表明，在一個(gè)局部區(qū)間，不僅可以用簡(jiǎn)單曲線來(lái)不光滑曲線，也可以代替復(fù)雜曲線，比如高頻震蕩曲線．這在實(shí)際應(yīng)用中是很有好處的，尤其在對(duì)曲線形狀和長(zhǎng)度比較敏感的場(chǎng)合．

4 算法分析

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

由于一般形式積分關(guān)系(15)的結(jié)果，形式非常復(fù)雜，難于計(jì)算，所以這里考慮用高斯積分?jǐn)?shù)值積分方法來(lái)計(jì)算．下面給出求解恰當(dāng)y?m的牛頓-拉夫遜算法(New ton-Raphson Method)．

4.1 算法

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

2) 計(jì)算幾何區(qū)域的慣性方向(或重心)，選取恰當(dāng)?shù)膞m∈[xl,xr]，計(jì)算ym相應(yīng)的初值

3) 計(jì)算對(duì)應(yīng)的S(x)：

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

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

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

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

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

(V) 若|L?Lm|>ε，則返回3)繼續(xù)計(jì)算；否則，終止跳出循環(huán)，結(jié)束計(jì)算；

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

4.2 算法分析

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

2) 理論上講，對(duì)任意的xm∈[xl,xr]都可以求出ym對(duì)應(yīng)的近似值，不過(guò)從后續(xù)計(jì)算的穩(wěn)定性和效率來(lái)看，在|k1|和|k2|比較接近的情況下，將xm選擇靠近區(qū)間的中點(diǎn)還是不錯(cuò)的方案，這樣得到的曲線，極性小些[14,19]．

4.3 數(shù)值結(jié)果

下面以函數(shù)f(x)=|sin(x)|,x∈[?1.5,1.5]為例來(lái)計(jì)算．該函數(shù)在x=0處導(dǎo)數(shù)不存在，有一個(gè)不光滑點(diǎn).在區(qū)間[?1,1]上面長(zhǎng)度有限．

圖7為選擇不同的插值節(jié)點(diǎn)xm來(lái)計(jì)算的光滑曲線．他們分別對(duì)應(yīng)xm分別為?2/3,0,2/3三個(gè)值的三對(duì)曲線．上面一組點(diǎn)線和下一組點(diǎn)線相對(duì)應(yīng)，粗線為原始非光滑曲線；小圓圈表示插值節(jié)點(diǎn)xm在曲線上的位置．每一組曲線在長(zhǎng)度、光滑程度均滿足要求，都可以作為局部的曲線近似．所以我們很容易根據(jù)實(shí)際需要來(lái)選擇相應(yīng)的曲線，這保證了方法的穩(wěn)定可靠．

圖7:不同x m對(duì)應(yīng)的連續(xù)光滑曲線

表1的結(jié)果在牛頓迭代終止條件為10?5下取得．迭代次數(shù)5～6次．曲線在區(qū)間[?1,1]上的真實(shí)長(zhǎng)度：2.622884996431094．計(jì)算結(jié)果僅保留小數(shù)點(diǎn)后面5位．

從表1可以看出，實(shí)際計(jì)算速度很快；而且對(duì)不同的xm，計(jì)算結(jié)果正確且算法穩(wěn)定，所以本文提出的等長(zhǎng)度曲線近似方法計(jì)算效果很好．為簡(jiǎn)單得到唯一的曲線，只需將插值節(jié)點(diǎn)xm選為區(qū)間的中點(diǎn)即可．此外，本文主要討論了連續(xù)曲線的逼近問(wèn)題，如果工程中需要對(duì)于間斷曲線進(jìn)行近似，可以有兩個(gè)辦法：一個(gè)是逐短光滑，另一種是整體光滑．這兩種問(wèn)題利用本文的處理方法不存在任何困難，因?yàn)楸疚奶岢龅慕品椒ㄖ恍枰涝诮o定區(qū)間上曲線的長(zhǎng)度即可，與曲線的具體幾何形態(tài)無(wú)關(guān)．

表1:不同x m對(duì)應(yīng)的連續(xù)光滑曲線

參考文獻(xiàn)：

[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]吳向紅，葉繼根，馬遠(yuǎn)樂(lè)，等.水平井蒸汽輔助重力驅(qū)油藏模擬方法[J].計(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]陶文銓.數(shù)值傳熱學(xué)[M].西安：西安交通大學(xué)出版社，2001 Tao W Q.Num erical Heat Transfer[M].X i’an:X i’an Jiaotong University Press,2001

[6]范曉光，馬學(xué)虎.微通道內(nèi)蒸汽及混合蒸氣冷凝[D].大連：大連理工大學(xué)，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]李慶揚(yáng)，王能超，易大義.數(shù)值分析[M].武漢：華中科技大學(xué)出版社，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]鐘爾杰，黃廷祝.數(shù)值分析[M].北京：高等教育出版社，2004 Zhong E J,Huang T Z.Num erical Analysis[M].Beijing:H igher Education Press,2004