***** Ｆｒａｃｔａｌ ＣＧｓ *****
1993/02/10 2000/07/09 中村 雅彦
※各図をクリックするとＪａｖａＳｃｒｉｐｔを実行します。
※使用方法やソースプログラム等は Ｆｒａｃｔａｌ ＣＧｓ 図一覧 をご覧ください。
1975年、B.B.Mandelbrot氏によって拓かれたフラクタル幾何学は、今日、広く市民権を得るようになりました。一口にフラクタル幾何学と言ってもいくつかの分野がありますが、特に「複素力学系」の分野においては、私達は不思議な図形を見ることができます。これはコンピュータに負うところが大きいわけですが、昨今のコンピュータの発展に伴って数々の作品が発表され、多くの人の目に触れるようになりました。また、誰でも簡単にこれらの図形を描くことができます。このレポートは、こういったフラクタル図形のいくつかを実際に描くことを目的として書いたもので、若干の理論とそれに基づく具体的な計算方法を示してみました。
参考文献
- B.B.Mandelbrot 「The fractal geometry of nature」
W.H.Freeman and Company, New York 1982
（広中平祐監訳 「フラクタル幾何学」 日経サイエンス社 1985）
- 宇敷重広 「フラクタルの世界」 日本評論社 1987
- H.-O.Peitgen P.H.Richter 「The Beauty of Fractals」
Springer-Verlag Berlin 1986
（宇敷重広訳 「フラクタルの美」 シュプリンガー・フェアラーク東京 1988）
- 伊藤貢司 伊東直基 「PC-8801/mkII/SR ３Ｄスーパーグラフィクス」
秀和システムトレーディング 1985
- 月刊 「ＰＣマガジン」 ラッセル社
- 数値解析学の講義
私が初めてフラクタル図形を見てから現在に至るまでの間、その他多数の書籍を拝読しました。
ある複素写像 f(z) が与えられたとき、初期値 z に対して
｛ f^0(z)=z, f^1(z)=f(z), f^2(z)=f・f^1(z), …, f^n(z)=f・f^{n-1}(z), … ｝で定められる z の軌道 {f^n(z)} を考えます。関数や初期値によって様々なものが考えられますが、n を大きくしていったときに f^n(z) がどのような値になっていくのかを見ていくことになります。言葉を変えると、初期値 z が f によって繰り返して写像を受けたとき、どのような振舞いをするのかを見ていくわけです。これには３つの場合が考えられます。
1. f^n(z) が ∞ に発散する場合
2. f^n(z) がある値に収束する場合
3. f^n(z) が発散も収束もしない場合
z=x+yi というように、複素数を実数部と虚数部に分けて扱うことにします。例えば z^2 を計算するときは、新たな実数部は x^2-y^2 となり、新たな虚数部は 2xyとなります。この程度だとまだ簡単ですが、変数が増えたり分数が出てきたりすると、だんだんと複雑な式になっていきます。
一般にコンピュータの画面には座標が入っていて、この座標を使って画面の任意の位置に点を表示することができます。さらに、普通は何色かの色を選ぶことができます。そこで、この画面を複素平面のある領域に対応させ、画面の１点につき対応する複素数を計算します。そして、この値を写像 f のパラメータ、または初期値として上記の反復写像を施します。一連の計算を行うと、結果として１つの値が出てきます。これは整数であったり、あるいは実数であったりするので、実数の場合は適当な法則にしたがって整数に変換し、最終的には各整数に応じて色を決定して画面に点を表示します。このような計算を画面のすべての点について行うと、不思議な図形が得られるわけです。
注意しなければならないのは、１つの図形を得るために行う計算の量が膨大になるということです。場合によっては数時間を要することもあります。このために計算時間を短縮するアルゴリズムもありますが、今回はできるだけ簡単なアルゴリズムを採用することにしました。それでも、少しでも速くなるように心がけています。
B.B.Mandelbrot氏の名を冠して呼ばれている集合です。一般には複素２次関数
fc(z)=z^2+cを考えます。氏本人は z→z^2-μ という写像を考えて「μ−ｍａｐ」と呼んでいますが、上式のように ＋ の方が有名になったようです。（ここでは μ の変わりに c を使っています。）
マンデルブロ集合の定義の仕方はいくつかありますが、その１つは
Ｍ ＝｛c in Ｃ｜n→∞ のとき fc^n(0)→∞ とならない.｝です。初期値 0 の fc による軌道が発散しないようなパラメータ c の集合です。ここに出てくる定数 0 は特異点と呼ばれ、fc'(z)=0 となる z の値です。これをもとに図を描きます。
数値計算で ∞ を扱うことはできないので、工夫する必要があります。ここでは、ある k に対して ｜fc^k(0)｜>2 ならば fc^n(0)→∞ となることが証明できるので、これを判定条件とします。実際には各 c につき、有限の n に対して 0<=k<=n の範囲で fc^k(0) を順次計算します。途中ある k について ｜fc^k(0)｜>2 となった場合 fc^n(0) は発散するので、この c は Ｍ に属さないと判定します。逆に最後まで ｜fc^k(0)｜<=2 であった場合、この c は Ｍ に属すると判定するのです。有限の n で計算を打ち切るため、得られる図形は近似的なものになります。n が小さいと正しい判断ができなくなり図形がくずれてしまうので、ある程度大きな値にする必要があります。図１が Ｍ の全体図で、-2.0<Re c<0.5, -1.25<Im c<1.25 の範囲を表示しています。計算の打ち切り回数 n は 50 としました。
さて、Ｍ に属さない部分（図１の白い部分）には、簡単に色付けすることができます。この部分は fc^n(0) が発散するところで、ある k に対して ｜fc^k(0)｜>2 が成り立ちました。そこで、この発散する速さ k に応じて色付けしてやります。例えば、白黒でしか表示できないときは、k の偶奇に応じて白黒に分けることで図２が得られます。カラーで表示できる場合は、綺麗な図形を見ることができるでしょう（タイトル）。
全体図を表示する分には n は 50 くらいで十分ですが、ある部分を拡大してみるとき、例えば「雪だるまの手の指」（図３, -0.25<Re c<0.0, 0.85<Im c<1.1）を表示するときは、n=100 くらいにしないと Ｍ の輪郭がはっきりしません。
マンデルブロ集合は連結で、それ自身と同じ形をしたものがあちこちで見られ（自己相似）、また、どんなに小さな領域を拡大表示してみても複雑で綺麗な図形が見られます。ただし、拡大するにつれて n を大きくしなければならず、それに比例して計算時間がかかるようになります。
氏が「μ−ｍａｐ」と並んで「λ−ｍａｐ」と呼んでいるものがあります。これもマンデルブロ集合と呼ばれてよいものだと思います。氏本人は z→λz(1-z) という写像を考えていますが、ここでは文字の統一をはかって
fc(z)=cz(1-z)を考えます。このマンデルブロ集合の定義は
Ｍ'＝｛c in Ｃ｜n→∞ のとき fc^n(0.5)→∞ とならない.｝です。0.5 は特異点を表します。計算方法は 3.1 の場合とまったく同じで、Ｍ' の全体図は図４のようになります（-2.0<Re c<4.0, -3.0<Im c<3.0）。この図形も、あちこちを拡大表示すると複雑な図形が見えてきます。
上記以外にも、fc(z) は任意の式にすることができます。これらはマンデルブロ集合と呼ばれるものではありませんが、同じ計算方法で様々な図形が得られます。一例として fc(z)=z^d+c を考えてみました。
Ｍd ＝｛c in Ｃ｜n→∞ のとき fc^n(0)→∞ とならない.｝を計算すると、例えば d=6 の場合 Ｍ6 は図５のような「ひとで」になります（-1.4<Re c<1.1, -1.25<Im c<1.25）。一般に Ｍd のひとでの手の数は d-1 本のようです。
マンデルブロ集合 Ｍ は通常べた塗りで表示されますが、これに色付けしようというものです（図２参照）。Ｍ に属さない部分は今までどおり発散する速さ k をもとに色付けしますが、Ｍ に属する部分については
α(c)＝inf｛｜fc^k(0)｜： 0<=k<=n ｝すなわち、各 c における ｜fc^k(0)｜ の最小値をもとに色付けします。α(c) は実数で出てくるので、ここでは大ざっぱに α(c) を 200 倍して、その整数部でグループ分けしてみました。こうしてできたものが図６です。
また、別の方法として
｜fc^k(0)｜＝α(c) のとき index(c)＝kとして、index(c) すなわち、どの k について｜fc^k(0)｜が最小になったかをもとに色付けすると図７が得られます。
G.Julia氏の名を冠して呼ばれている集合です。ジュリア集合は、マンデルブロ集合のように特定の集合を表すのではなく、様々な写像 f についてそれぞれのジュリア集合が存在します。
ジュリア集合の定義の１つは
Ａf(∞) ＝｛z in Ｃ｜n→∞ のとき f^n(z)→∞ ｝としたとき、
Ｊf ＝ δＡf(∞) (…Boundary)です。複素平面の各点を初期値として反復写像を施したとき、軌道が ∞ に発散する点の集合の境界と言うことができます。以下で f の例をいくつか挙げてみます。
マンデルブロ集合のところでも出てきた複素２次関数
f(z)=z^2+c
を考えます。前の場合、c は計算する点の座標だったので各点ごとに式が変化しており、また、初期値は 0 に固定されていました。ここでは、パラメータ c は任意の値に固定されており、初期値は複素平面の各点の座標となります。具体的な計算方法はマンデルブロ集合のところとほとんど同じで、式および初期値が若干変わる程度です。有限の n に対して 1<=k<=n の範囲で f^k(z) を順次計算して、途中ある k について ｜f^k(z)｜>3 となったら z は Ａf(∞) に属すると判定し、そうでなければ z は Ａf(∞) に属さないと判定します。マンデルブロ集合のところでは ２ を判定基準としましたが、ここでは ３ であればよいことが証明できます。図８は c=-1.0+0.0i（-1.8<Re z<1.8, -1.8<Im z<1.8）の場合を表示しており、白い部分が Ａf(∞) を表し、白と黒の境界部分がジュリア集合 Ｊf を表しています。また、この黒い部分（Ｃ＼Ａf(∞)）を充填ジュリア集合と呼ぶこともあります。
さらにここでもマンデルブロ集合のところと同じように、Ａf(∞) の各点をその発散する速さに応じて色付けすることによって、図９を得ることができます。定義からは外れますが、これらの図を単にジュリア集合と言うこともあります。
ここでのパラメータ c をいろいろと変えることで、まったく違った様々な図形を見ることができますが（図10, c=0.3+0.03i, -1.5<Re z<1.5, -1.5<Im z<1.5）、c の値によっては n を非常に大きくしなければならないものもあります。また、このジュリア集合とマンデルブロ集合の関連づけとして、c in Ｍ のとき Ｊf は連結であるということが分かります。さらに、マンデルブロ集合付近のある範囲の拡大図と、そこからパラメータ c を選んで描いたジュリア集合が似ているというのも不思議なことです。（図11, c=-0.1+1.0i, -1.8<Re z<1.8, -1.8<Im z<1.8, 図３と比較）
f の例を４つほど挙げてみます。式によっては判定条件が変わっているものもあります。
f(z)=z^d+c
（図12, d=6, c=-1.0+0.1i, -1.5<Re z<1.5, -1.5<Im z<1.5）
f(z)=cz(1-z)
（図13, c=3.0+0.1i, -0.1<Re z<1.1, -0.6<Im z<0.6）
f(z)=a exp(bz)
（図14, a=0.8, b=0.6, -3.0<Re z<6.0, -4.5<Im z<4.5,
判定条件は exp および sin, cos の計算ができる範囲かどうか）
f(z)=c sin(z)
（図15, c=-1.0+0.5i, -2pie<Re z<2pie, -2pie<Im z<2pie,
判定条件は sin, cos の計算ができる範囲かどうか）
充填ジュリア集合のべた塗りの部分に色付けしようというものです。この部分に含まれる点 z には、ある点 z0 に収束していくものがあります。また、ある程度写像を繰り返すと、…, z1, z2, …, zn, z1, z2, … といった n 周期点にのるものもあります。これから、何回の繰り返しである値になったかをもとに色付けすることができます。実際には誤差を考えて、何回の繰り返しである値の近傍に入ったかを計算します。「ある値」は実際の計算結果をもとにすればよいでしょう。例えば f(z)=z^2-1（図９参照）の場合は、これらの点は …, 0, -1, 0, -1, … という２周期点にのります。図16は、何回目に -1 になったかで色付けしてみました。
他の例として２つほど挙げてみます。
f(z)=z*(z+c)/(1+conj(c)*z) ブラシュケ関数
図17は c=0.2+0.5i（-2.0<Re z<2.0, -2.0<Im z<2.0）の場合で、単位円の外側の点は ∞ に発散し、内側の点は 0 に収束していきます。
図18は c=-0.7+1.2i（-2.0<Re z<2.0, -2.0<Im z<2.0）の場合で、各点は 0.857-0.515i に収束していきます。なんとなく単位円が見えるでしょうか。
f(z)=z^3-3pz+c ミルナー関数
図19は p=0.3, c=0.3+0.1i（-1.7<Re z<1.7, -1.7<Im z<1.7）の場合で、∞ に発散するものと 0.159+0.055i に収束するものとに分かれます。
図20は p=-0.4, c=0.1+0.0i（-1.7<Re z<1.7, -1.7<Im z<1.7）の場合で、∞ に発散するものの他に、0.163+0.529i および 0.163-0.529i に収束するものがあります。ムンクを見ることができるでしょう。
（実数値の）非線形方程式 f(x)=0 を数値的に解く方法の１つに、ニュートン法があります。これは、適当な初期値 x0 から出発して１根 α に収束する列 ｛xk｝を作り、xk が α に十分に近づいたとき計算を打ち切って、xk を近似解とするものです。ここで、
x{k+1} = xk - f(xk) / f'(xk)です。これを複素領域に拡張するとどうなるかという問題です。簡単な例として、f(x)=x^3-1 を考えてみます。実数の範囲で考えると
x{k+1} = xk - (xk^3-1) / 3xk^2 = (2xk^3+1) / 3xk^2となって、x0 として 0 以外の任意の実数を取った場合、xk は 1 に収束していきます。しかし、複素数の範囲で考えて文字 x をあらためて z=x+yi とすると
z{k+1}=x{k+1}+y{k+1}iとなって、z0 として 0 以外の任意の複素数を取った場合、zk が 1 に収束するものの他に、(-1+sqrt(3)i)/2 や (-1-sqrt(3)i)/2 に収束するものがあります。
=(2xk^5+4xk^3yk^2+xk^2+2xkyk^4-yk^2)/(3(xk^4+2xk^2yk^2+yk^4))
+(2yk(xk^4+2xk^2yk^2-xk+yk^4))/(3(xk^4+2xk^2yk^2+yk^4))i
図21は、初期値を -2.0<Re z0<2.0, -2.0<Im z0<2.0 の範囲で取り、zk が 1 に収束したものを白で表示しています。白黒ではできませんが、zk が収束した値によって３色に塗り分けると、各色の境界では２色で接することがなく、必ず３色で接しているという不思議なものです。また、ここでも収束の速さに応じて色付けすることができます。
図22は，３つの解のいずれに収束したかで分類し、さらにその収束の速さで色付けしてあります。同様に、
図23は z^4-1=0 の例で（-2.0<Re z0<2.0, -2.0<Im z0<2.0）、各点は 1,-1,i,-i のいずれかに収束します。
図24は (z-1)(z^2+z+1/2)=0 について、-2.42<Re z0<-1.85, -0.285<Im z0<0.285 の範囲を表示したものですが、前の２例とは違った判定方法で色分けしてみました。（明らかにかおです。）この場合は zk を順次計算して、２点 zkと z{k-1} の距離が小さな値 eps よりも小さくなったところで収束したと判定して、その k に応じて色付けしています。このときの zk を表示すると、それが解であるわけです。こちらが本当のニュートン法による解法であって、前の２例はすでに解が分かっている上で作ったものです。
再びマンデルブロ集合に戻って話をします。図２や図３では Ｃ＼Ｍ の部分を発散する速さに応じて色付けしました。Ｍ を一様に帯電した金属とみなすならば、この各色の境界線は等電位曲線を表します。これはまた地図の等高線と見ることもできるので、立体表示をすることができます。より一般に、
ψ(c)＝lim(1/2^k)log｜fc^k(0)｜とすると、上の等高線は ψ の近似となります。ψ は、Ｍ の付近では 0 になり、Ｍ から離れるにしたがって増加していく、連続した正の値を取ります。これを高さに見立てて立体表示をしてみます。
プログラム「mandelp」は、複素平面の各点について高さを計算して一括出力するものです。実際の計算は、やはり有限回 n で打ち切ります。最終的には上下を逆さまにしてマンデルブロ集合が頂上になるようにしています。また、適当な高さになるように伸び縮みさせています。プログラム「mapolyn」は、入力された一連の高さデータをもとに地図を表示するものです。汎用的なものなので、高さデータを変えることで普通の地図でも表示することができます。プログラムの内容は本題から外れるので説明は省略しますが、この名前は「地図を、ポリゴンを使って、スムーズシェイディングを施して表示する」という意味のつもりです。
任意の関数 f のジュリア集合についても、
ψ(z)＝lim(1/2^k)log｜f^k(z)｜を考えることで立体表示をすることができます。プログラム「juliap」は、f(z)=z^2+c についての例です。この場合は、より立体感を出すために ψ(z) を log(1+log(1+log(ψ(z)))) のように変換してみました。
●フラクタルに関するレポート（PC-9801用原文）
fractcgs.lzh 500,664バイト
readme はこちら
ファイルの内容一覧はこちら
●レポート文書をＰＤＦ化したものはこちらです。（表紙・図なし）