剥いだ布団とメモランダム.old

情報系のことをかいてゆく

Python + gnuplot でマンデルブロ集合を描く

Python シェルスクリプト カオス

寒いです。冬ですね。
この前、複雑系の授業の課題でマンデルブロ集合を描いてこいっていう課題がでまして。

最初にPythonのmatplotlibを使ってプロットしようと思ったんだけども、
あれって点数が増えてくるとむちゃくちゃ重くなるんっすね。諦めました。

そんなわけで、

  • Pyrhonでプロットデータ出力
  • gnuplotでグラフ描画＆画像ファイル出力

っていうプロセスでやってみました。

実際は、Pythonのコードを書いて、
シェルスクリプトでPythonコードの実行とgnuplotでのグラフ出力をまとめて実行してます。

マンデルブロ集合とは？

\left\{ \begin{array}{1} z_{n+1} = z_n^2 + c \\ z_0 = 0 \end{array} \right.
という漸化式で定義されるz_nがn\rightarrow \inftyで発散しない複素数cの集合のこと。

つまり、この集合を求めるには、複素数cを色々な値に変えて、
漸化式をたくさん回していって発散するかどうか調べればいいということ。

ググると色付けされたきれいなフラクタルがたくさんでてくる。
これは、cの実部と虚部をそれぞれの軸にとって、
発散した時のnの数によって等高線のように色付けしたものらしい。
（色付けの基準は他にも色々とあるみたいだけど、これがスタンダード）

発散とか\inftyとか実装できないんじゃ…

課題が出て一番最初に思ったのはこれ。

実際にプログラムでは擬似的に、

  • |z_n|がある定数Kを超えたら発散と判断
  • nがあるループ上限を超えたら発散しなかったと判断

とやってあげる。定数Kは割りと小さい値（2や3）でもいい感じでした。

ソースコード

まずはプロットデータを書き出すPythonの方から。
【mandel.py】

# -*- coding: utf-8 -*-
import sys

def mandelbrot(c, K, LOOPMAX):
	#発散まで回したループ数を返す
	n = 0
	z = 0.0 + (0.0 * 1j)
	while (n<LOOPMAX and abs(z)<K):
		z = z**2 + c
		n += 1
	return n

def getArg():
	#コマンドライン引数用関数
	argv = sys.argv
	argc = len(argv)
	if (argc != 9):
		#引数がない場合は終了	
		print ('Usage: python %s filename '
			'minReal maxReal minImage maxImage step loop emission') %argv[0]
		quit()
	return argv

def main():
	argv = getArg()
	
	writer = open(argv[1], 'w')

	loop = int(argv[7]); K = float(argv[8])

	step = float(argv[6])
	x_min = float(argv[2]); x_max = float(argv[3])
	y_min = float(argv[4]); y_max = float(argv[5])

	x = x_min
	while (x < x_max):
		y = y_min
		while (y < y_max):
			c = x + (y * 1j)
			n = mandelbrot(c, K, loop)
			writer.write(str(c.real)+'\t'+str(c.imag)+'\t'+str(n)+'\n')
			y += step
		writer.write('\n')
		x += step
	print x_max

if __name__ == '__main__':
	main()

Pythonは複素数をそのまま扱えるのでいいよね。
複素数使えない言語だと先述の漸化式の二次元実写像を考えればok（ググって）。
途中で改行を入れるのは、gnuplotでデータを扱うときに
1番目の変数（x軸の方）の値が変わるときに改行が必要だからです。


次にPythonのコードを実際に動かして、gnuplotのコマンドを打ち込むシェルスクリプト。
【mandelplot.sh】

#!/bin/sh
#引数チェック
if [ $# -ne 8 ]; then
	echo "Usage: $0 filname minReal maxReal minImage maxImage step loop emission" 1>&2
	exit 1
fi

#pythonによるマンデルブロ集合データ出力
python mandel.py $1.data $2 $3 $4 $5 $6 $7 $8

# #gnuplotによるプロット
echo "set term png; 
set output \"$1\"; 
set grid;
set pm3d map; 
set size square; 
splot \"$1.data\"" | gnuplot

#openはmacで使える
open $1

exit 0

gnuplotのコマンドについてメモしておきます。

set term png          #pngファイルを出力するように設定
set output "filename" #出力ファイル名を設定 
set grid              #グラフにグリッド線を描画するように設定
set pm3d map          #等高線（pm3d）の2次元プロット（mapオプション）
set size square       #グラフエリアを正方形に
splot "plotdata"      #プロットデータのあるファイルを指定＆プロット

実行結果

$ ./mandelplot.sh mtest.png -2 0.6 -1.2 1.2 0.0015 20 3

f:id:sheemaa:20131211141402p:plain

実数部を-2〜0.6の範囲、虚数部を-1.2〜1.2の範囲で0.0015おきにプロット。
（ループ上限は20回、発散判断値は3）


いい感じ。ただ、プロットする座標値の変化数(0.0015)はもうちょい大きい値でもよさそう。
ちょっと拡大したものを出力してみる。


$ ./mandelplot.sh mtest.png zoomed.png -0.16 0.02 0.93 1.03 0.00015 70 3

f:id:sheemaa:20131211141653p:plain
実数部を-0.16〜0.02の範囲、虚数部を0.93〜1.03の範囲で0.0015おきにプロット。
（ループ上限は70回、発散判断値は3）

拡大してみると、さっきの全体像と似た模様が小さく見える。これがフラクタル性。
もう少し拡大してみる。


$ ./mandelplot.sh zoomed2.png -0.018 -0.01 1.016 1.024 0.00001 70 3

f:id:sheemaa:20131211142317p:plain
実数部を-0.018〜-0.01の範囲、虚数部を1.016〜1.024の範囲で0.00001おきにプロット。
（ループ上限は70回、発散判断値は3）

すごい。宇宙を感じる（小学生並みの感想）

まとめ

色付け工夫してみるときれいな図が描けそうですよね。
パステルカラーとかおしゃれになりそう。


f:id:sheemaa:20131211143705p:plain
おまけ。（gnuplotコマンド出力のecho文を下のように変更）

# #gnuplotによるプロット
echo "set term png; 
set output \"$1\"; 
set grid;
set hidden3d;
set pm3d;
set size square;
splot \"$1.data\" with lines;" | gnuplot