今回のテーマ

Mathematicaショートフラクタルシリーズ、第3回は、フラクタルの説明で最も使われているであろう、コッホ曲線です。

ほぼ完成形

koch_original.nb

generator = {0, Pi/3, -Pi/3, 0};
f[list_] := list + # & /@ generator
inclinationList = Flatten[Nest[f, generator, 5]];
cPoints = Join[{0}, Accumulate[E^(I inclinationList)]];
Graphics[Line[Transpose[{Re[cPoints], Im[cPoints]}]]]

偏角関数

今回のコッホ曲線描画アルゴリズムは、
僕が学部生のときに発表した研究内容に基づきます。

自主課題研究「偏角関数を用いた曲線の研究」(2006) (PPT発表スライド)

この図は、コッホ曲線の第0段階から第2段階までを左から走査して横軸を長さ、縦軸を角度とした偏角関数を示したものです。

偏角関数の
第0段階は $0$
第1段階は $0, \frac\pi3, -\frac\pi3, 0$
第2段階は $0, \frac\pi3, -\frac\pi3, 0, \frac\pi3, \frac{2\pi}3, 0, \frac\pi3, \dots$

という列と捉えられます、$\frac\pi3$ で除算すれば
第1段階 $0$
第2段階 $0, 1, -1, 0$
第3段階 $0, 1, -1, 0, 1, 2, 0, 1, \dots$

となり見やすくなります。
（Mathematicaで実装する際は最初から $\frac\pi3$ が掛かっててもいいし、整数列を作って単位を掛けてもいいです。Cとか一般的な言語だと後者が賢明でしょうが、今回は前者です。）

フラクタルのパーツと偏角関数の関係は以下の図でより伝わるかと思います。

生成過程

generator = {0, Pi/3, -Pi/3, 0};

これが第1段階です。

f[list_] := list + # & /@ generator

これは偏角リストにジェネレーターの各定数を足し、

{
  {0, Pi/3, -Pi/3, 0} + 0,
  {0, Pi/3, -Pi/3, 0} + Pi/3,
  {0, Pi/3, -Pi/3, 0} + -Pi/3,
  {0, Pi/3, -Pi/3, 0} + 0
}

と計算します。（リスト同士の加算は次元を合わせる必要がありますが、リストと定数の加算では各要素に定数が足されます。）

inclinationList = Flatten[Nest[f, generator, 5]];

で第6段階まで進みます。
最終的に一気に平坦化するので、その時点で一連の偏角リストが得られます。

そして次のように偏角から複素平面上の座標を構築します。
オイラーの公式より、 $e^{a i}$ が1辺先となります。
（先頭に原点の0を忘れずに）

cPoints = Join[{0}, Accumulate[E^(I inclinationList)]];

複素数の魔法で曲線をプロットします。

Graphics[Line[Transpose[{Re[cPoints], Im[cPoints]}]]]

ショートコード化

ワンライナー化（2重の無名関数をうまく作れなかったのでやむを得ずFunctionを用意）
AccumulateとJoinを併用するよりPlusをFoldListした方が端的

結果、115文字のプログラムができました。

koch_short.nb

Graphics@Line@Transpose@{Re[#],Im[#]}&@FoldList[Plus,0,E^(Flatten@Nest[Function[l,l+#&/@{0,Pi/3,-Pi/3,0}],{0},5]I)]

@butchi_y (More info: https://t.co/HBUuaZOB3w) #wolframlang pic.twitter.com/qj9AvcqEFN
— Tweet-a-Program (@wolframtap) 2015, 12月 15

デモ

ちなみに、これをもっと一般化したデモをWolframサイトに掲載しています。
これまた有名なフラクタル図形を多数再現できるほか、manualモードで角度も自由に変えられます。
僕のフラクタル研究の集大成のひとつですね。

Angular Function Fractals - Wolfram Demonstrations Project

予告

明日はドラゴン曲線です。