淺析MATLAB環境下分形圖形的生成

馬洪強

引言

19世紀末20世紀初，研究者在探究一些特殊曲線性質時，發現了處處連續但處處不可微曲線，1904年Koch根據這種性質構造出了Koch曲線，該曲線是分形理論初期的研究成果之一。1890年，意大利數學家皮亞諾（Peano G）發現了能填滿一個正方形皮亞諾曲線。20世紀學者對分形理論進行了深入研究，Besicovitch等人研究了曲線的維數、分形集的結構、分形集的局部性質及在數論、測度論等方面的應用，特別是維數理論成果豐富。隨著維數理論、投影理論、隨機技巧等方法先后建立并趨于成熟[1]，使得分形集合的研究形成了自己的方法體系，這一時期分形理論初見雛形。1915年波蘭著名數學家希爾賓斯基（Sierpinski）提出了Sierpinski 三角形及Sierpinski 地毯[2]，1919年法國數學家嘉斯頓·茹利亞（Gaston Julia）研究了迭代保角變換[3]，成為后來人們借助計算機展現一幅幅精彩的圖畫的理論基礎。 20世紀70年代，美籍法裔猶太科學家曼德勃羅（Mandelbrot）在總結前人的基礎上發表了《分形：形狀、機遇和維數》，第一次系統闡述分形幾何的思想、內容、意義和方法，標志著分形幾何作為獨立學科正式誕生。特別是1980年Mandelbrot借助計算機將Mandelbrot集（M-集）展示在熒屏上，標志著分形藝術初步形成，隨著分形軟件的開發應用，分形圖形以其絢麗多彩的色彩感沖擊著人們的視覺感受。現代藝術也漸漸將分形思想應用到創作中，甚至部分藝術家直接創作分形作品進行展覽。

2004年2月3日，分形幾何之父Gaston Julia誕辰111周年，Google網站logo結合Julia集分形圖形做了創新設計，以紀念這位數學家。

本文分四個部分，第一部分介紹了什么是分形；第二部分介紹MATLAB軟件.M文件運行步驟；第三部分在MATLAB（Matrix Laboratory）環境下，以Koch曲線為例介紹分形圖形的產生過程；第四部分介紹Julia集和M集分形實例的MATLAB程序；第五部分闡述Sierpinski地毯的向量化編程技術及其變異，并將向量化編譯程序隨機化，拓寬了分形圖形生成方法；第六部分對向量式隨機分形程序的意義和不足做了說明，并闡述了隨機分形意義所在。

一、 分形

1.1 分形定義

迄今為止，分形還沒有嚴格意義上的定義。1982年，曼德勃羅定義“Hausdorff維數嚴格大于拓撲維數的幾何稱為分形”，四年后又提出“組成部分以某種方式與整體相似的形體叫分形”，但依然不夠確切和嚴格。現階段多數引用法爾內克（K.Falconer）對分形集合F的描述：

1）F具有精細的結構，即在任意小的尺度之下，它總有復雜的細節；

2）F是不規則的，以至它的整體和局部不能用傳統的幾何語言來描述；

3）F通常具有某種自相似性，這種自相似性可以是統計意義上的，也可以是近似的；

4）F在某種意義下的分形維數通常都大于它的拓撲維數；

5）在多數情況下，F以非常簡單的方法定義，或以遞歸或迭代過程產生。

1.2 分形圖形分類

分形可以分為規則分形和不規則分形。一些數學家提出過不少復雜和不光滑的集合，例如Cantor集、Koch曲線、Sierpinski三角、地毯和海綿等，都屬于規則分形，它們具有嚴格的自相似性。但自然界許多事物的不光滑性和復雜性往往是隨機的，像蜿蜒曲折的海岸線，隨機布朗運動等，這類曲線的自相似性是近似的或統計意義上的。

下面以Koch曲線為例，闡述Matlab環境下，分形圖形的生成過程。

二、 Matlab軟件介紹

MATLAB軟件是美國Mathworks公司開發的一款優秀的數學軟件，MATLAB是Matrix Laboratory(矩陣實驗室)的組合縮寫。該款軟件將數據分析、矩陣計算可視化，并集成了非線性動態系統的仿真。MATLAB軟件不僅適用于學術研究，還在模擬仿真等現實工作中得到應用，其版本以其應用廣泛，不斷得到更新，目前已更新到MATLAB 15.0版本。本文以經典版本MATLAB 7.0為例，簡單介紹在分形應用時MATLAB軟件應用步驟。

MATLAB軟件界面主要包含工具欄、命令窗口(Command Windows)、命令歷史窗口(Command History)、工作間管理窗口（Workspace）、當前路徑窗口(Current Directory)、編譯窗口、圖形窗口等。

運行分形.M文件有多種方式，本文簡單介紹打開.M文件運行方式。在Matlab界面中選擇工具欄文件菜單中打開選項，按存儲路徑選中.M文件，打開.M文件后單擊運行按鈕即可得到分形圖。

使用.M文件運行程序比較簡單直接，便于保存、修改設計程序。[4]

三、 Matlab環境下的Koch曲線分形實現

本文以科赫曲線為例，在闡述分形的數學思想基礎上來詳細說明分形原理及在Matlab環境下分形圖形的實現過程。

Koch曲線的構造是將單位為1的一條線段變成四條，且每條線段為1/3，因此變換n次后的Koch曲線總長度為（4/3）^n，當n→∞，則Koch曲線趨于無窮大，Koch曲線的分形維數：D=ln4/ln3=1.2618。

Koch曲線的遞歸過程如圖3.1所示：

圖3 .1

圖3 .2

3.1 算法分析

如圖3.2所示，考慮由直線段（2個點）產生第一個圖形（5個點）的過程。圖中，設p1和p2分別為原始直線段的兩個端點，現需要在直線段的中間依次插入三個點P2，P3，P4；P2，P4分別位于線段三分之一處，線段三分之二處，P3點的位置可看成是由點P4以點P2為軸心，逆時針旋轉60°而得。旋轉由正交矩陣實現。

根據初始數據（P1和P5點的坐標），產生圖中5個結點的坐標。結點的坐標數組形成一個矩陣，矩陣的第一行為P1的坐標，第二行為P2的坐標……，第五行為P5的坐標。矩陣的第一列元素分別為5個結點的x坐標，第二列元素分別為5個結點的y坐標。

進一步考慮Koch曲線形成過程中結點數目的變化規律。設第k次迭代產生的結點數為n(k)，第k+1次迭代產生的結點數為n(k+1)，則n(k)和n(k+1)中間的遞推關系為n(k+1)=4n(k)-3。

3.2 Matlab程序實現

p=[0 0;10 0]; %P為初始兩個點的坐標，第一列為x 坐標，第二列為y坐標

n=2; %n為結點數

A=[cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3)]; %旋轉矩陣

for k=1:9

d=diff(p)/3; %diff計算相鄰兩個點的坐標之差， 得到相鄰兩點確定的向量%則d就計算出每個向量長度的三分之一，與題中將線段三等分對應

m=4*n-3; %迭代公式

q=p(1:n-1，:); %以原點為起點，前n-1個點的坐標為終點形成向量

p(5:4:m，:)=p(2:n，:); %迭代后處于4k+1位置上的點的坐標為迭代前的相應坐標

p(2:4:m，:)=q+d; %用向量方法計算迭代后處于4k+2位置上的點的坐標

p(3:4:m，:)=q+d+d*A′; %用向量方法計算迭代后處于4k+3位置上的點的坐標

p(4:4:m，:)=q+2*d; %用向量方法計算迭代后處于4k位置上的點的坐標

n=m; %迭代后新的結點數目

end

plot(p(:，1)，p(:，2)) %繪出每相鄰兩個點的連線

axis([0 10 0 6]) %設置坐標顯示

3.3 Matlab程序運行結果

Koch曲線Matlab示例見參考文獻[5]。

四、 Julia集和M-集分形MATLAB程序

Julia集和M-集分形MATLAB程序設計有很多種方法，為便于藝術工作者應用MATLAB程序實現分形，此處介紹了較方便的程序設計，在程序運行時間上也較省時。

4.1 Julia集MATLAB程序

a=-0.13;b=0.77;r=2;

for x0=-1:0.01:1

for y0=-1:0.01:1

x=x0;y=y0;

if x0^2+y0^2＜1

for n=1:80

x1=x*x-y*y+a;

y1=2*x*y+b;

x=x1;

y=y1;

end

if (x*x+y*y)＜r

plot(x0，y0);

end

hold on;

end

end

end

4.2 M-集 MATLA程序

r=1; %限界值

for a=-1:0.002:1

for b=-1:0.002:1 %參數a，b取到一個范圍

x=a;y=b; %初始的復數c

for n=1:20

x1=x*x-y*y+a; %復數平方加一個c的運算

y1=2*x*y+b;

x=x1; %迭代

y=y1;

end

if(x*x+y*y)＜r %限界

plot(a，b);

end

hold on;

end

end

五、 Sierpinski地毯

5.1 矩陣法 Sierpinski地毯

圖5 .1.1

Sierpinski地毯是具有代表性的經典分形，可有迭代算法產生。將一個正方形n^2個相等的邊長為1/n的小正方形，依據一定規則挖去k^2個小正方形，得到Sierpinski地毯的生成元，重復上述過程就可得到Sierpinski地毯。

利用MATLAB的imagesc函數，將正方形用矩陣中的方陣模擬，構造迭代矩陣，對方陣進行填充，繪制Sierpinski地毯分形圖，具體程序如下：

function M=Sierpinski(i)

tic

M=0

for j=1:i

M=[M， M， M;

M， ones(3^(j-1))*1， M;

M， M， M];

end

imagesc(M);

colormap(spring);

axis equal;

axis off;

toc

將上述程序保存為Sierpinski.m文件，在MATLAB的command Window中輸入M=Sierpinski(4)得到圖形5.1.1此例詳見參考文獻。[6]

在構造迭代矩陣時，將迭代矩陣構造成不完全對稱矩陣，就得到變異Sierpinski地毯分形圖，如圖5.1.2所示，程序如下：

function M=Sierpinski(i)

tic

M=0

圖5 .1.2

for j=1:i

M=[ones(3^(j-1))*9， M， M;

M， M， M;

M， M， ones(3^(j-1))*9];

end

imagesc(M);

colormap(autumn);

axis equal;

axis off;

toc

雖然利用imagesc函數，通過構造迭代矩陣的方法簡單明了，便于計算，更易通過改變矩陣形式，產生新的分形圖。但在實際運行過程中，受到一定因素的限制，只能進行有限的迭代，對于提高分形圖形的精度及細膩程度存在不足。該方法僅為一般狀態下的規則分形，沒有考慮隨機因素，因此圖形呈現比較單一。

5.2 Sierpinski地毯向量式隨機程序生成分形圖

基于3.1Sierpinski地毯變異程序基礎上，添加隨機因素，使其影響分形生成結果，從而使得分形圖形產生變化。本例應用matlab中rand函數產生隨機數，并依概率采用有差別化向量法程序，程序如下

function M=SS (i)

tic

M=0

for j=1:i

r=rand

if r＜0.5;

N=M;

M=[ones(3^(j-1))*9， M， N;

M， N， M;

N， M， ones(3^(j-1))*9];

else，

N=ones(3^(j-1))*9;

M=[ones(3^(j-1))*9， M， N;

M， N， M;

N， M， ones(3^(j-1))*9];

end

end

imagesc(M);

colormap(autumn);

axis equal;

axis off;

toc

將程序存儲為SS.m函數，在command window中輸入M=SS(4)將得到不同圖形，本例給出部分分形結果，如圖5.2所示:

圖5 .2

六、 展望

雖然向量式隨機分形程序設計解決了向量式隨機分形程序的隨機問題，但是受制于矩陣法缺陷，有限迭代問題沒有解決。此外，對于程序設計等也值得研究，以便得到改進，縮短運行時間。因此也為此類方法的前景及研究提出了新的問題。

從應用領域來說，向量式隨機分形程序設計方法使得分形結果產生隨機變化，為分形圖形的不確定性提供了可能，并使MATLAB環境下分形圖形與藝術設計領域相關學科相互交叉，從而在分形藝術應用方面提供了可能；同時，藝術設計的發展需求也對分形提取更高要求。