フラクタル図形を描こう　その3～Excel編～

どうも、フラクタル図形大好きおじさんと化しつつあるチオールです。

前回までは主にペイントを使用したフラクタル図形の書き方を紹介していましたが、今回は全く新しい方法を紹介したいと思います。

hassium277.hatenablog.com
hassium277.hatenablog.com
※☝前回

今回使用するのは、　Microsoft Excel　です。

グラフの章

Excelは入力・計算したデータを様々な形式のグラフに変換して可視化することができます。

で、そのグラフの種類の中に「散布図」というものがあります。

f:id:Hassium277:20210217202000p:plain

散布図はその名の通り平面上に点をバラまく表示法で、2個ずつの実数からなるデータを可視化する方法です。

そして、この散布図とIFS(Iterated function system：反復関数系)というものを組み合わせることで、フラクタル図形を描画することができます。
ja.wikipedia.org
※「IFS」でググって初めて知ったんですが、エクセルの機能として「IFS関数」というものがあるそうです。今回使うものとは無関係です。

Wikipediaの説明を読んでもよくわからない人が多いと思うので、まずは単純な例を一つ紹介します。

まず、A1のセル(左上)に0を入力します。

f:id:Hassium277:20210219194031p:plain

次に、右隣のB列に上から順番に

=IMPRODUCT(A1,1/2)

=IMSUM(IMPRODUCT(A1,1/2),1/2)

=IMSUM(IMPRODUCT(A1,1/2),COMPLEX(0,1/2))

と入力します。

f:id:Hassium277:20210219201819p:plain

これはA1の値に対して

$A1×\frac{1}{2}$

$A1×\frac{1}{2}+\frac{1}{2}$

$A1×\frac{1}{2}+\frac{1}{2}i$

という計算をすることを表しています。

以下、この3つの関数をそれぞれf₁、f₂、f₃と表記します。

つまり、

$\begin{cases}f_1(z)=\frac{1}{2}z\\f_2(z)=\frac{1}{2}z+\frac{1}{2}\\f_3(z)=\frac{1}{2}z+\frac{1}{2}i\end{cases}$

という関数を定義します。

そしてC列以降は、前の列の値をf₁、f₂、f₃に代入した値を配置します。

f:id:Hassium277:20210219201738p:plain
f:id:Hassium277:20210219202431p:plain

これをひたすら繰り返します。

f:id:Hassium277:20210219204616p:plain

気が済むまで繰り返したら、IMREAL関数とIMAGINARY関数を使って計算結果の実部と虚部を分離します。

f:id:Hassium277:20210219205306p:plain

f:id:Hassium277:20210219205332p:plain

そして実部と虚部を散布図としてプロットすると、

f:id:Hassium277:20210219205916p:plain

ようやくフラクタル図形とご対面できます。

さて、ここまでの説明を読んで「シェルピンスキーのギャスケットなんてペイントでも書けるしその方が簡単だろ」と思った人もいるかと思います。

それは本当にその通りなんですが、Excelの底力はこんなもんじゃありません。

「表計算ソフト」という二つ名が表す通り、Excelは計算に秀でたソフトです。

ドットの描画しかできないペイントと異なり、三角関数などの複雑な計算を自動で行ってくれます。

また、特定のセルの内容を変更した場合、そのセルを計算に使っているセル、そしてそれを使ったグラフにも変更が反映されます。

つまりどういうことかというと、

f:id:Hassium277:20210220153259g:plain

こういうことができます。

これは

$\begin{cases}f_1(z)=\frac{1}{2}z(\text{cos}\theta+i\text{sin}\theta)\\f_2(z)=\frac{1}{2}z(\text{cos}\theta+i\text{sin}\theta)+\frac{1}{2}\\f_3(z)=(\frac{1}{2}z+\frac{1}{2}i)(\text{cos}\theta+i\text{sin}\theta)\end{cases}$

として、 $\theta$ を1°ずつ変化させた画像をgif動画にしたものです。

gif動画の作成には、「Giam」というソフトを使用しました。

forest.watch.impress.co.jp

また、使用する関数を変えることで、以下のように全然違う見た目のものをつくることができます。

f:id:Hassium277:20210221150834g:plain

これは以下の関数から生成されます。

$\begin{cases}f_1(z)=z(\frac{1}{2}+\frac{1}{2}i)\\f_2(z)=(z(\text{cos}\theta+i\text{sin}\theta)+1)(\frac{1}{2}+\frac{1}{2}i)\end{cases}$

VBAの章

Excelには作業を自動化するための「マクロ」という機能があり、「VBA」という言語でプログラミングができます。

フラクタル図形の描画に関しては、以下のようなことができるようになります。

数百回単位の繰り返し処理
繰り返し処理の繰り返し処理の繰り返し処理の・・・というような、繰り返しの多重化
条件に応じた処理の分岐*1
セルの自動塗りつぶし

ところで、「マンデルブロ集合」というフラクタル図形をご存知でしょうか。

この記事をここまで読んでいる人はおそらくフラクタル図形が好きな人だと思いますが、そういった方なら一度は目にしたことがあるぐらい有名な名前だと思います。

マンデルブロ集合の定義は以下の通りです。

$z_0=0$ 、 $z_{n+1}=z_n^2+c$ とするとき、
「 $\displaystyle{\lim_{n\rightarrow\infty}}|z_n|$ が無限大に発散しないような複素数 $c$ 全体の集合」
をマンデルブロ集合という。

ja.wikipedia.org

これを実際に描画するには、以下のような計算をします。

(1)：平面上から点 $(a,b)$ を1つ選ぶ。

(2)：以下の数列をn=100～5000ぐらいまで計算する。計算途中で $x_n^2+y_n^2$ が4を超えたら、計算をやめる。*2

$\begin{cases}x_0=0,x_{n+1}=x_n^2-y_n^2+a\\y_0=0,y_{n+1}=2x_ny_n+b\end{cases}$

(3)：中断せずに指定した回数まで計算できたら、平面上の $(a,b)$ の位置に点を打つ。*3

(4)：(1)に戻る。

マクロを使わないとできないことが目白押しですね。

というわけで実際の作業手順の説明をしたいところですが、VBA自体の説明は私の手に負えないと判断したので、他のサイトに丸投げします。

ちなみに私は以下のサイトで勉強しました。

excel-ubara.com

最終的に必要なソースコードはこの記事内で公開するので、とりあえず第1～7回を読んでVBAを使用する準備だけはしておいてください。

というわけで、ここからの説明はExcelの起動からVBAの実行までを自力でできる人向けの説明になります。

まず、準備としてセルの縦横の幅を変更して、全てのセルを正方形にしてください。

Ctrl+Aを押して全てのセルを選択し、セル番号の欄を左クリックすることで全ての行/列の幅を一括で変更できます。

f:id:Hassium277:20210223202144p:plain — ☝「列の幅」を選択するとこうなる

f:id:Hassium277:20210223202224p:plain — ☝列の幅を0.05に変更

f:id:Hassium277:20210223202338p:plain — ☝列の幅が変更され、マス目が視認できない程小さくなった

私のPCでは列の幅を0.05、行の幅を0.5にするとちょうどよくなりますが、どうやらこの組み合わせは環境によって異なるようです。

設定が完了したら、VBEを開いて以下のコードを張り付けてください。

Sub f()
    Dim x, y, px, py As Double
    Dim o As Boolean
    For a = 0 To 200
        For b = 0 To 200
            x = 0
            y = 0
            o = True
            For c = 0 To 100
                px = x
                py = y
                x = px ^ 2 - py ^ 2 + a / 50 - 2
                y = 2 * px * py + b / 50 - 2
                If x ^ 2 + y ^ 2 > 4 Then
                    Cells(b + 1, a + 1).Interior.Color = (1 - (1 - c / 100) ^ 9) * 255
                    o = False
                    Exit For
                End If
            Next c
            If o Then Cells(b + 1, a + 1).Interior.Color = 0
        Next b
    Next a
End Sub

これを実行すると、

f:id:Hassium277:20210224222426p:plain

マンデルブロ集合が描画されます。

ちなみに、描画された範囲のセルをコピーしてペイントに貼り付けることで、描画した図形を画像に変換することができます。

このとき、先述の「列の幅を0.05、行の幅を0.5」にしているとExcelでの1セル=ペイントでの1ドットに変換されます。

f:id:Hassium277:20210223205519p:plain — ペイントに貼り付けた。なお、同じ手順でパワーポイントやワードに貼り付けると大変なことになる。

さて、先程のコードでは「描画サイズ」「描画範囲」「計算回数」は固定でしたが、これらを変えることでマンデルブロ集合の真の姿を見ることができ、美しきフラクタルの世界を堪能できるようになります。

描画サイズの変更

先程のコードで生成される画像のサイズは201×201でしたが、これを変更するには以下の部分の数値を定数倍すればいいです。

f:id:Hassium277:20210224222452p:plain

例えば以下のように変更することで、生成される画像が2001×2001*4になります。

Sub f()
    Dim x, y, px, py As Double
    Dim o As Boolean
    For a = 0 To 2000
        For b = 0 To 2000
            x = 0
            y = 0
            o = True
            For c = 0 To 100
                px = x
                py = y
                x = px ^ 2 - py ^ 2 + a / 500 - 2
                y = 2 * px * py + b / 500 - 2
                If x ^ 2 + y ^ 2 > 4 Then
                    Cells(b + 1, a + 1).Interior.Color = (1 - (1 - c / 100) ^ 9) * 255
                    o = False
                    Exit For
                End If
            Next c
            If o Then Cells(b + 1, a + 1).Interior.Color = 0
        Next b
    Next a
End Sub

描画範囲の変更

最初のコードでは、描画範囲は「実部と虚部がそれぞれ-2～2」となっていました。

この描画範囲を小さくする、つまり細部にズームインしていくことで、マンデルブロ集合の深淵に踏み込んでいくことが可能になります。

まず、コードを以下のものに書き換えてください。

Sub f()
    Dim x, y, px, py, cx, cy, r As Double
    Dim o As Boolean
    r = 1 / 2
    cx = 0
    cy = 0
    For a = 0 To 200
        For b = 0 To 200
            x = 0
            y = 0
            o = True
            For c = 0 To 100
                px = x
                py = y
                x = px ^ 2 - py ^ 2 + a / (100 * r) + cx - 1 / r
                y = 2 * px * py + b / (100 * r) + cy - 1 / r
                If x ^ 2 + y ^ 2 > 4 Then
                    Cells(b + 1, a + 1).Interior.Color = (1 - (1 - c / 100) ^ 9) * 255
                    o = False
                    Exit For
                End If
            Next c
            If o Then Cells(b + 1, a + 1).Interior.Color = 0
        Next b
    Next a
End Sub

2行目で宣言する変数が3つ追加され、その下で3つとも初期化されています。

新しい変数のうち cx と cy は描画範囲の中心座標、 r は拡大率です。

例えば

    r = 2
    cx = 0
    cy = 1 / 2

とすると、

f:id:Hassium277:20210225000025p:plain

$\frac{1}{2}i$ を中心として先程までのものを4倍に拡大した画像が生成されます。

これでマンデルブロの真の姿が拝めますね。

なお、拡大の条件によっては以下のような形が見えることがあります。

f:id:Hassium277:20210225233653p:plain

実はこれ、マンデルブロ集合の真の姿ではありません。

このことには先述の計算手順の(2)が関係しています。

(1)：平面上から点 $(a,b)$ を1つ選ぶ。
(2)：以下の数列をn=100～5000ぐらいまで計算する。計算途中で $x_n^2+y_n^2$ が4を超えたら、計算をやめる。
$\begin{cases}x_0=0,x_{n+1}=x_n^2-y_n^2+a\\y_0=0,y_{n+1}=2x_ny_n+b\end{cases}$
(3)：中断せずに指定した回数まで計算できたら、平面上の $(a,b)$ の位置に点を打つ。
(4)：(1)に戻る。

マンデルブロ集合に属さない点の中には発散が遅いものが存在するため、有限回の計算ではどうしても判別漏れが生じてしまいます。

根本的で完璧な解決策は多分存在しないと思いますが、有限の解像度の画像を生成するだけなら単純に計算回数を増やすだけで対処できます。

具体的には、以下のように書き換えることで真の姿を拝むことができます。

Sub f()
    Dim x, y, px, py, cx, cy, r As Double
    Dim o As Boolean
    r = 2048
    cx = -1 / 8 + 1 / 32 + 1 / 128 - 1 / 512 - 1 / 1024
    cy = 1 - 1 / 32 - 1 / 128 + 1 / 512 - 1 / 2048
    For a = 0 To 200
        For b = 0 To 200
            x = 0
            y = 0
            o = True
            For c = 0 To 500
                px = x
                py = y
                x = px ^ 2 - py ^ 2 + a / (100 * r) + cx - 1 / r
                y = 2 * px * py + b / (100 * r) + cy - 1 / r
                If x ^ 2 + y ^ 2 > 4 Then
                    Cells(b + 1, a + 1).Interior.Color = (1 - (1 - c / 500) ^ 9) * 255
                    o = False
                    Exit For
                End If
            Next c
            If o Then Cells(b + 1, a + 1).Interior.Color = 0
        Next b
    Next a
End Sub

f:id:Hassium277:20210226004013p:plain

上記のコードでは

            For c = 0 To 100

が

            For c = 0 To 500

に変わっています。

ここの部分の c という変数は $x_n$ 、 $y_n$ の $n$ に対応する変数で、上限を100から500に書き換えることは計算回数を増やすことを意味します。

なお、変更点はもう1か所あって、

                    Cells(b + 1, a + 1).Interior.Color = (1 - (1 - c / 100) ^ 9) * 255

この部分が

                    Cells(b + 1, a + 1).Interior.Color = (1 - (1 - c / 500) ^ 9) * 255

になっているという違いがあります。

この部分はマンデルブロ集合の外側に色を塗る役割のコードですが、ここにも c が使われているので修正する必要があるというわけです。

というわけで、次は色の塗り方について解説したいと思います。

カラーリングの変更

先述の通り、ソースコードの以下の部分は「cの値から色を決めて塗る」という動作を表しています。

                    Cells(b + 1, a + 1).Interior.Color = (1 - (1 - c / 100) ^ 9) * 255

以下、 (1 - (1 - c / 100) ^ 9) * 255 の部分を彩色関数と呼びます。

VBAでは黒と赤の間の色が256段階に分割され、0～255までの整数で番号が付けられています。

この彩色関数では、0～100までの c の値を0～255までの色に割り当てています。

f:id:Hassium277:20210227133614p:plain — ☝cの値(横軸)を色(縦軸)に対応させている

c の範囲を0～500に増やした場合、色の計算結果が0~255の範囲からハミ出すため、この部分も修正する必要があったというわけです。

ちなみに、彩色関数を (1 - (1 - c / 100) ^ 9) * 255 から (1 - (1 - c / 500) ^ 9) * 255 に変更した場合、全体的に暗い色になります。

f:id:Hassium277:20210227221100p:plain

拡大率が小さいと暗すぎて見づらくなりますが、拡大していくと明るすぎて見辛かった模様が見えるようになります。

f:id:Hassium277:20210227221530p:plain — ☝像が正確になるだけでなく、外側が暗くなったことで枝状の模様が見えるようになった

さて、彩色関数を変えれば見た目が変わるという(ほぼ自明な)事実が確認できたわけですが、実はこれよりも大幅に変更することで印象をガラリと変えてしまうことができます。

というわけで、 (1 - (1 - c / 100) ^ 9) * 255 以外の彩色関数の例を紹介します。

基礎パターンの派生

色に関する関数として、VBAではRGB(r,g,b)という関数が使えます。

これは3つの整数 r 、 g 、 b を使って色を指定する関数で、それぞれの変数は光の3原色である赤、緑、青に対応しています。

そして、これを先程まで説明で使っていた彩色関数と組み合わせることで、全体の色を変えることができます。

例えば

Cells(b + 1, a + 1).Interior.Color = RGB(0, (1 - (1 - c / 100) ^ 9) * 255, 0)

とすると

f:id:Hassium277:20210227232102p:plain

緑になったり、

Cells(b + 1, a + 1).Interior.Color = RGB((1 - (1 - c / 100) ^ 9) * 255, (1 - (1 - c / 100) ^ 9) * 255, (1 - (1 - c / 100) ^ 9) * 255)

に変えることで

f:id:Hassium277:20210227234426p:plain

白にできたりします。

※カラーリングの違いをより楽しんでいただくため、画像の生成サイズを2001×2001にしています。

全分割

f:id:Hassium277:20210223235818p:plain

Cells(b + 1, a + 1).Interior.Color = 16777215 * c / 100

VBAでは使用可能な色に対して0～16777215の番号が振られており、こちらはcの値を0～16777215の間に等間隔に割り振ったものです。

このままだとちょっと気持ち悪い配色ですが、最後に100で割っている所を500で割るようにすると緑と黒といういかにも今風な配色になります。

f:id:Hassium277:20210224002103p:plain

晴天

f:id:Hassium277:20210227235012p:plain

Cells(b + 1, a + 1).Interior.Color = RGB(cr(c * 9), cr(c * 9 + 85), cr(c * 9 + 170))

これを使うには、以下のコード(cr関数の定義)を末尾に追加する必要があります。

Function cr(ByVal c As Long) As Integer
    cr = (c Mod 256) * (255 - (c Mod 256)) / 65
End Function

先程までの彩色方法では、cの上限を上げる度にコードを修正する必要があり、さらに全体の明るさが変わってしまうという弱点がありました。

一方この方法なら、剰余を使っているためcが大きくなっても値が一定の範囲に収まるようになっています。

そのため、この彩色関数は発散の遅い領域の描画に適しています。

f:id:Hassium277:20210228171026p:plain

Sub f()
    Dim x, y, px, py, cx, cy, r As Double
    Dim o As Boolean
    r = 1024
    cx = -1 + 1 / 4 + 1 / 512
    cy = 1 / 16 + 1 / 512
    For a = 0 To 2000
        For b = 0 To 2000
            x = 0
            y = 0
            o = True
            For c = 0 To 10000
                px = x
                py = y
                x = px ^ 2 - py ^ 2 + a / (1000 * r) + cx - 1 / r
                y = 2 * px * py + b / (1000 * r) + cy - 1 / r
                If x ^ 2 + y ^ 2 > 4 Then
                    Cells(b + 1, a + 1).Interior.Color = RGB(cr(c * 9), cr(c * 9 + 85), cr(c * 9 + 170))
                    o = False
                    Exit For
                End If
            Next c
            If o Then Cells(b + 1, a + 1).Interior.Color = 0
        Next b
    Next a
End Sub

Function cr(ByVal c As Long) As Integer
    cr = (c Mod 256) * (255 - (c Mod 256)) / 65
End Function

雪解け

f:id:Hassium277:20210228173108p:plain

Cells(b + 1, a + 1).Interior.Color = RGB(cr(c * 9), cr(c * 9 + 85), cr(c * 9 + 170))

曇天のような灰色、雪のような白、泥水のような茶色という4月の北海道のようなカラーリングですが、拡大していくとカラフルになります。

f:id:Hassium277:20210228155229p:plain

色が計算回数の上限に依存しない等、使い勝手は先ほどのものとあまり変わらないので、好みに応じて使い分けると良いかと思います。

逆鱗

f:id:Hassium277:20210228155258p:plain

こちらは彩色関数を変えるだけでなく、ソースコード全体を書き換えています。*5

Sub f()
    Dim x, y, z
    For a = 1 To 2000
        For b = 1 To 2000
            x = 0
            y = 0
            For c = 0 To 500
                z = x
                x = z ^ 2 - y ^ 2 + a / 500 - 2
                y = 2 * z * y + b / 500 - 2
                If x ^ 2 > 9 Then Exit For
            Next c
            Cells(b, a).Interior.Color = (x + y) * 99
        Next b
    Next a
End Sub

このコードは「なるべく少ない文字数でマンデルブロ集合を描画する」という遊びをしていた際に偶然発見したもので、どういう仕組みでこんな配色になるのかは全く分かりません。

また、このコードのようにxy座標の値を彩色関数に使用すると、以下のようなエラーが出ることが少なからずあります。

f:id:Hassium277:20210228175110p:plain

多分これは使用した色の数が多すぎることが原因で、要するに

Excel「16777215色使えるとは言ったが、全部同時に使っていいとは言ってない」

という事だと思います。

描画対象の変更

ここまで描画サイズ、拡大率、カラーリングの変え方を紹介してきましたが、いよいよ描画するフラクタル図形の種類そのものを変える方法を説明します。

マンデルブロ集合の定義に使われる $z_{n+1}=z_n^2+c$ という計算は、ソースコードのこの部分で行っています。

px = x
py = y
x = px ^ 2 - py ^ 2 + a / 50 - 2
y = 2 * px * py + b / 50 - 2

px と py が $z_n$ 、x と y が $z_{n+1}$ に、 a / 50 - 2 と b / 50 - 2が $c$ に対応しています。

VBAは(Excelは対応しているのに何故か)複素数の直接計算に対応していないため、実部と虚部に分けて計算しています。

この部分を他の式に変えることで、別のフラクタル図形を描画することができます。

マルチブロ

f:id:Hassium277:20210228204115p:plain

px = x
py = y
x = px ^ 3 - 3 * px * py ^ 2 + a / 500 - 2
y = 3 * px ^ 2 * py - py ^ 3 + b / 500 - 2

en.wikipedia.org

マルチブロ集合はマンデルブロ集合の指数を2以外の値に一般化したもの、すなわち「 $z_0=0$ 、 $z_n=z_{n+1}^d+c$ という数列の絶対値が無限大に発散しない複素数 $c$ 全体の集合」です。

上記のものは $d=3$ のもので、 $(x+iy)^d$ を計算して実部と虚部に分けた式を直接入力すれば $d$ を4以上の自然数にしたものも描画できます。

f:id:Hassium277:20210228205209p:plain

px = x
py = y
x = px ^ 4 - 6 * px ^ 2 * py ^ 2 + py ^ 4 + a / 500 - 2
y = 4 * px ^ 3 * py - 4 * px * py ^ 3 + b / 500 - 2

$d$ が3以上の整数のとき、 $d$ 次のマルチブロ集合は実軸について線対称で、なおかつ $d-1$ 回回転対称(正 $d-1$ 角形と同じ対称性)になるようです。

トリコーン

f:id:Hassium277:20210228210311p:plain

px = x
py = -y
x = px ^ 2 - py ^ 2 + a / 500 - 2
y = 2 * px * py + b / 500 - 2

en.wikipedia.org

トリコーンは、マンデルブロ集合の漸化式の $z_n$ を $\bar{z}_n$ に変えたものです。

マンデルブロ集合をマルチブロ集合に一般化するのと同様に、トリコーンも一般化可能です。

$d$ が2以上の整数のとき、 $d$ 次のマルチコーンは実軸について線対称で、なおかつ $d+1$ 回回転対称になるようです。

バーニングシップ・フラクタル

f:id:Hassium277:20210315215126p:plain

px = Abs(x)
py = Abs(y)
x = px ^ 2 - py ^ 2 + a / 500 - 2
y = 2 * px * py + b / 500 - 2

ja.wikipedia.org

「バーニングシップ」は直訳すると「燃える船」です。

Wikipediaの画像は白黒なのでわかりにくいですが、赤で塗るとかなり炎っぽくなりますね。

充填ジュリア集合

f:id:Hassium277:20210228222313p:plain

Sub f()
    Dim x, y, px, py As Double
    Dim o As Boolean
    For a = 0 To 200
        For b = 0 To 200
            x = a / 50 - 2
            y = b / 50 - 2
            o = True
            For c = 0 To 5000
                px = x
                py = y
                x = px ^ 2 - py ^ 2 - 0.1
                y = 2 * px * py + 0.65
                If x ^ 2 + y ^ 2 > 4 Then
                    Cells(b + 1, a + 1).Interior.Color = RGB(cr(c * 7), cr(c * 8), cr(c * 9))
                    o = False
                    Exit For
                End If
            Next c
            If o Then Cells(b + 1, a + 1).Interior.Color = 0
        Next b
    Next a
End Sub

Function cr(ByVal c As Long) As Integer
    cr = (c Mod 256) * (255 - (c Mod 256)) / 65
End Function