遞歸(分片)仿射分形插值數(shù)值模擬與盒維數(shù)研究*

袁利國(guó)， 余榮忠， 曠菊紅

(1.華南農(nóng)業(yè)大學(xué) 數(shù)學(xué)與信息學(xué)院,廣州 510640；2.九江學(xué)院 理學(xué)院,江西 九江 332005； 3.五邑大學(xué) 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院,廣東 江門(mén) 529020)

0 引 言

傳統(tǒng)的歐氏幾何主要分析規(guī)則的幾何體，如直線(xiàn)、雙曲線(xiàn)、圓、立方體等。而對(duì)不規(guī)則的對(duì)象，如海岸線(xiàn)的長(zhǎng)度、自然界的各種形狀，傳統(tǒng)的幾何很難進(jìn)行處理，而分形幾何的創(chuàng)立，對(duì)大自然的研究提供了新的思路與方法，分形幾何是描述大自然的幾何，其分形維數(shù)可以更好地衡量不規(guī)則體。分形幾何引起廣泛的研究及應(yīng)用[1-10]。迭代函數(shù)系統(tǒng)(IFS)在分形幾何中扮演重要的作用，IFS可以生成各種分形圖，如:Koch曲線(xiàn)、Sierpinski三角形、Cantor集、各種植物等。IFS也可應(yīng)用到插值與擬合極不規(guī)則的數(shù)據(jù)，相比傳統(tǒng)的Newton插值與Lagrange插值，有很強(qiáng)的優(yōu)勢(shì)。基于IFS的分形插值函數(shù)與初等函數(shù)一樣有自身的幾何特征，能用一組IFS表示出來(lái)，迭代這組IFS，得到其吸引子，它就是分形插值函數(shù)，通常它具有分?jǐn)?shù)維數(shù)，是非整數(shù)的。在插值已知數(shù)據(jù)集時(shí)，由于吸引子函數(shù)(即分形插值函數(shù))有分?jǐn)?shù)維，因此可以插值或擬合不規(guī)則的震蕩的數(shù)據(jù)集，比傳統(tǒng)的插值方法有優(yōu)勢(shì)。常規(guī)的分形插值是把整個(gè)區(qū)間上的數(shù)據(jù)I=[a,b]，在IFS的某個(gè)壓縮仿射變換Wk作用下，分別映射到其每個(gè)子區(qū)間Ik=[xk-1,xk]?I上，從而構(gòu)造一組仿射變換Wk，再基于確定的迭代算法或隨機(jī)的迭代算法，得到此組IFS的吸引子曲線(xiàn)，此曲線(xiàn)會(huì)穿過(guò)已知的插值數(shù)據(jù)點(diǎn)，其迭代函數(shù)系中仿射變換的系數(shù)公式為

1 遞歸仿射分形插值

給定區(qū)間I=[0,1](或其他的一般區(qū)間I=[a,b]也可以)，令0=x0xi-xi-1。設(shè)K=[0,1]×R,wi:K→K是仿射映射:

其中,|di|<1,wi把(xl,yl)與(xm,ym)分別映射到(xi-1,yi-1)與(xi,yi)，或把(xl,yl)與(xm,ym)分別映射到(xi,yi)與(xi-1,yi-1)。定義關(guān)聯(lián)矩陣C=(Cij)，如果Jj?Ji′，則Cij=1，否則Cij=0。設(shè)I(i)={j:Cij=1}，定義映射W=(wij)=(Cijwi)。

(1)

遞歸分形插值要求滿(mǎn)足xm-xl>xi-xi-1，此時(shí)有以下3種情形：(1)Ji?Ji′ ；(2)Ji?Ji′且Ji∩Ji′=Φ；(3)Ji?Ji′，但是Ji∩Ji′≠Φ。如果屬于第(1)種情形，則遞歸分形插值也是通常的分片(或稱(chēng)分段、分塊)分形插值。如果是第(2)或第(3)種情形，即在文獻(xiàn)[1]末尾給出的算例情形，將逐一給出數(shù)值分析。

2 遞歸仿射分形插值數(shù)值模擬

遞歸仿射分形插值本質(zhì)上也是分塊方法[3-4]。文獻(xiàn)[5]給出一類(lèi)多參數(shù)遞歸分形插值函數(shù)的數(shù)學(xué)模型,并進(jìn)行數(shù)值模擬，但給出的代碼是仿射分形插值的情形。文獻(xiàn)[1]末尾給出的算例的遞歸結(jié)構(gòu)如圖1。此處采用區(qū)間I=[0,1.2]，采取等距分割，令

x0=0,x1=0.2,x2=0.4,x3=0.6,x4=0.8,x5=1.0,x6=1.2。設(shè)對(duì)應(yīng)y值y0=0,y1=0.8,y2=0.8,y3=0.6,y4=1.0,y5=1.6,y6=0.8。Ji=[xi-1,xi](i=1,2,…,6)，3個(gè)Ji′分別為[x0,x2]，[x2,x4]，[x4,x6]。對(duì)應(yīng)的6個(gè)仿射變換:w1(x0,y0)=(x0,y0),w1(x2,y2)=(x1,y1),與w1類(lèi)似[1],由圖1可得其他的仿射變換wi(i=2,3，…，6)。

基于圖 1的關(guān)系，可以分別求得6個(gè)仿射變換wi(i=1,2,…,6)的系數(shù)，其中di是垂直比例因子(可自由選取)。此處選取d1=0.5,d2=0.6,d3=0.7,d4=0.8,d5=0.6,d6=0.5。對(duì)應(yīng)6組仿射變換的系數(shù)為a1=0.5,e1=0,c1=2×(1-d1),f1=0a2=0.5,e2=0.4,c2=-0.5-2d2,f2=0.8a3=0.5,e3=0.6,c3=1.5-0.5d3，f3=0.4-0.6d3；a4=0.5,e4=0,c4=-0.5d4，f4=0.8-0.6d4a5=0.5,e5=0.6,c5=-2+0.5d5，f5=3.2-1.4d5；a6=0.5,e6=0.2,c6=1+0.5d6，f6=-0.2-1.4d6，其中，L-1=1/0.4。

圖1 遞歸仿射變換結(jié)構(gòu)Fig. 1 Structure of recurrent affine transformation

此遞歸結(jié)構(gòu)的關(guān)聯(lián)矩陣為

下面提出一類(lèi)隨機(jī)迭代算法，更簡(jiǎn)單且能達(dá)到遞歸分形插值目的。

(1) 初始化，設(shè)定自由參數(shù)dn(n=1,2,…,6)、最大迭代次數(shù)與插值數(shù)據(jù)點(diǎn)(x0,y0),(x1,y1),…,(x6,y6)等。

(2) 依據(jù)數(shù)據(jù)計(jì)算出迭代函數(shù)的系數(shù)an,en,cn,fn，得到wn(n=1,2,…,6)。

(3) 給定初始值(取(x1,y1)=(0,0))，判定x1落在哪個(gè)區(qū)間Ji上。因此，每個(gè)Ji上對(duì)應(yīng)有兩個(gè)仿射變換，產(chǎn)生一個(gè)隨機(jī)數(shù)，決定由其中的哪一個(gè)wi作用到(x1,y1)，生成新的點(diǎn)(x(k),y(k))，同時(shí)，把(x(k),y(k))賦值給(x1,y1)，再重復(fù)上面的過(guò)程。

(4) 直至達(dá)到最大迭代次數(shù)，則輸出所得數(shù)據(jù)，即遞歸分形插值曲線(xiàn)。詳細(xì)程序代碼見(jiàn)http://ylgdiy.bokee.com/507551074.html。所得插值結(jié)果如圖2所示。文獻(xiàn)[1]中的遞歸分形插值構(gòu)造，其中的所有變換wi,i=1,2,…,6滿(mǎn)足要求xm-xl>xi-xi-1，且是情形(1)Ji?Ji′與情形(2)Ji?Ji′且Ji∩Ji′=Φ的混合形式。

圖2 遞歸仿射分形插值曲線(xiàn)Fig. 2 Recurrent affine fractal interpolation curve

3 分片仿射分形插值

分片(或稱(chēng)分塊、分段)形插值是遞歸分形插值的特殊情形[2,8-9]，有聯(lián)系也有一定的區(qū)別。下面分兩種情況討論，一是分片區(qū)間之間是沒(méi)有重疊情形。另一種是分片區(qū)間有重疊情形，現(xiàn)構(gòu)造如下有重疊情形的分片分形插值，如圖3所示。它符合遞歸分形插值的要求:xm-xl>xi-xi-1，且有Jj?Ji′，是第(1)種情形。Ji′,i=1,2,3，三區(qū)間無(wú)重疊區(qū)域。對(duì)應(yīng)的6個(gè)仿射變換都采用把(xl,yl)與(xm,ym)分別映射到(xi-1,yi-1)與(xi,yi)。此遞歸結(jié)構(gòu)的關(guān)聯(lián)矩陣為

圖3 分片仿射分形插值結(jié)構(gòu)(無(wú)重疊)Fig. 3 Structure of piecewise affine fractal interpolation (no overlap)

分別求得6個(gè)仿射變換wi(i=1,2,…,6)的系數(shù)，其中di(i=1,2,…,6)是垂直比例因子，自由選取。d1=0.6,d2=0.8,d3=0.65,d4=0.7,d5=0.66,d6=0.6,6組系數(shù)分別為a1=0.5,e1=0,c1=2×(1-d1),f1=0;a2=0.5,e2=0.2,c2=-2d2,f2=0.8;a3=0.5,e3=0.2,c3=-0.5-0.5d3,f3=1-0.6d3;a4=0.5,e4=0.4,c4=1-0.5d4,f4=0.2-0.6d4;a5=0.5,e5=0.4,c5=1.5+0.5d5,f5=-0.2-1.4d5;a6=0.5,e6=0.6,c6=-2+0.5d6,f6=3.2-1.4d6。圖4為分片仿射分形插值曲線(xiàn)(無(wú)重疊)。采用隨機(jī)迭代算法，Matlab程序見(jiàn)http://ylgdiy.bokee.com/507560391.html。

圖4 分片仿射分形插值曲線(xiàn)(無(wú)重疊)Fig. 4 Piecewise affine fractal interpolation curve (no overlap)

針對(duì)分片區(qū)間有重疊情形，J1′與J2′有部分重疊區(qū)域?，F(xiàn)構(gòu)造如下遞歸(分片)分形插值重疊情形(圖5)，它符合遞歸分形插值的要求:xm-xl>xi-xi-1，且此時(shí)有Jj?Ji′，是第(1)種情形。

圖5 分片仿射分形插值結(jié)構(gòu)(重疊)Fig. 5 Structure of piecewise affine fractal interpolation (overlap)

對(duì)應(yīng)的8個(gè)仿射變換:w1(x0,y0)=(x0,y0),w1(x4,y4)=(x1,y1),其他的類(lèi)似。因此，可以分別求得仿射變換wi(i=1,2,…,8)的系數(shù)，其中di(i=1,2,…,8)是垂直比例因子，自由選取。具體取d1=0.6,d2=0.5,d3=0.4,d4=0.3,d5=0.3,d6=0.6,d7=0.5,d8=0.4,8組系數(shù):a1=0.25,e1=0,c1=1-1.25d1,f1=0；a2=0.25,e2=0.2,c2=-1.25d2,f2=0.8；a3=0.25,e3=0.4,c3=-0.25-1.25d3,f3=0.8;a4=0.25,e4=0.6,c4=0.5-1.25d4,f4=0.6;a5=0.25,e5=0.3,c5=-0.25,f5=0.9-0.8d5;a6=0.25,e6=0.5,c6=0.5,f6=0.4-0.8d6;a7=0.25,e7=0.7,c7=0.75,f7=0.7-0.8d7;a8=0.25,e8=0.9,c8=-1.0,f8=2-0.8d8。相應(yīng)的分形插值曲線(xiàn)如圖6 所示，其中重疊區(qū)域的插值曲線(xiàn)更復(fù)雜。相應(yīng)的程序代碼可參考http://ylgdiy.bokee.com/507560394.html。

圖6 分片仿射分形插值曲線(xiàn)(重疊)Fig. 6 Piecewise affine fractal interpolation curve (overlap)

此遞歸結(jié)構(gòu)的關(guān)聯(lián)矩陣為

4 仿射分形插值的盒維數(shù)

仿射分形插值曲線(xiàn)的盒維數(shù)有如下結(jié)論:

(2)

的解，否則分形維數(shù)是1[3-4,7]。

但現(xiàn)實(shí)中的插值數(shù)據(jù)通常不是等距的。此時(shí)求解式(2)的精確解是很難?，F(xiàn)在提出把分形維數(shù)的求解轉(zhuǎn)化為最優(yōu)化問(wèn)題式(3)，即

(3)

再采用粒子群最優(yōu)化算法處理問(wèn)題式(3),現(xiàn)給出兩個(gè)算例。

最大迭代次數(shù)為80。重復(fù)運(yùn)行10次程序，取10次結(jié)果的平均值，即得到最優(yōu)分形維數(shù)值為D=1.758 146 555 910 46。二者之間的誤差達(dá)到了10-13，其中某一次的收斂過(guò)程如圖7所示，很快可以收斂到真實(shí)值，大概在第10步左右。

圖7 分形維數(shù)收斂過(guò)程Fig. 7 Convergence process of fractal dimension

算例2 設(shè)插值點(diǎn){(0,0),(30,50),(60,40),(100,10)}，選取d1=0.5,d2=-0.5,d3=0.23,則a1=0.3,a2=0.3,a3=0.4。精確的分形維數(shù)D=1.180 2[10]。構(gòu)造目標(biāo)函數(shù)：

最大迭代數(shù)為80。重復(fù)運(yùn)行10次程序，取10次結(jié)果的平均值，即得到最優(yōu)分形維數(shù)值為D=1.180 163 827 574 234。二者之間的誤差達(dá)到了10-5，其中某一次的收斂過(guò)程如圖8所示。Matlab程序詳情見(jiàn)網(wǎng)頁(yè)http://ylgdiy.bokee.com/507560407.html。算例1的程序比這個(gè)簡(jiǎn)單，只需修改相應(yīng)的目標(biāo)函數(shù)即可。

圖8 分形維數(shù)收斂過(guò)程Fig. 8 Convergence process of fractal dimension

5 結(jié)束語(yǔ)

對(duì)于數(shù)據(jù)量較大的插值擬合，如經(jīng)濟(jì)數(shù)據(jù)、農(nóng)業(yè)數(shù)據(jù)、自然現(xiàn)象的長(zhǎng)時(shí)間數(shù)據(jù)等，常采用分片的思想來(lái)處理，應(yīng)用也較廣，其本質(zhì)上也是遞歸分形插值思想。遞歸仿射分形插值應(yīng)用到實(shí)際數(shù)據(jù)量較大的情形比較少見(jiàn)。對(duì)遞歸分形插值與分片分形插值的數(shù)值模擬給出了分析與算法程序，這為支撐相關(guān)理論分析及應(yīng)用到實(shí)際問(wèn)題提供了補(bǔ)充。程序可稍微修改，更具有一般的適用性與通用性。隨后也對(duì)仿射變換的盒維數(shù)求解給出最優(yōu)化求解思路，此方法可取得很好的結(jié)果。