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

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

Python + gnuplot でジュリア集合を描く

かなりまえのエントリーでマンデルブロ集合を描くっていのをやりましたけど、
今度はジュリア集合を描いてみよっていうことで、
前につくったプログラムをちょっと変えてジュリア集合を描いてみました。

また例によって、

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

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

ジュリア集合とは?


z_{n+1} = z_n^2 + c \\
という漸化式で定義されるz_nn\rightarrow \inftyで発散しない複素数zの初期値z_0の集合のこと。
(cは複素数の定数)
つまり、今度はcではなくz_0を変えていってグラフをプロットする、と。

cは適当に複素数を決めておけば、いいんだけども、
このcの値によっていろいろと変わったグラフがプロットされるのがこいつ。

ソースコード

この前のマンデルブロ集合をプロットするときにつかったソースコードを少し変えただけです。
ちょっとした説明は前回をみればたぶん大丈夫。

【julia.py】

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

def juliaSet(z, c, K, LOOPMAX):
	#発散まで回したループ数を返す
	n = 0
	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 != 10):
		#引数がない場合は終了	
		print ('Usage: python %s filename c_real c_image'
			'minReal maxReal minImage maxImage step loop') %argv[0]
		quit()
	return argv

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

	loop = int(argv[9]); K = 100.0

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

	c = float(argv[2]) + (float(argv[3]) * 1j)

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

if __name__ == '__main__':
	main()


【juliaplot.sh】

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

#pythonによるジュリア集合データ出力
python julia.py $1.data $2 $3 $4 $5 $6 $7 $8 $9

# #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

前回同様、シェルスクリプトからpythonのコードとgnuplotコマンドを呼び出しています。

実行結果

まずは、c = -0.74759 + 0.083587 iとしたときのプロット。
プロット範囲は実部・虚部共に1.6~-1.6で、0.06毎に、ループ上限7000回で動かしてみます。

$ ./juliaplot.sh julia01 -0.74759 0.083587 -1.6 1.6 -1.6 1.6 0.006 7000

f:id:sheemaa:20140605211551p:plain

個人的にはジュリア集合っていうとこんな感じのなんだけれども、
パラメータを変えていくとけっこう面白い形が見えてきます。

c = 0.250166 - 0.000004 iでプロットしてみる。

$ ./juliaplot.sh julia02 0.250166 -0.000004 -1.5 1.5 -1.5 1.5 0.005 3000

f:id:sheemaa:20140605211604p:plain
さっきとは全く違った形に。


最後に、c = -0.125475 + 0.7494 iでプロット。

$ ./juliaplot.sh julia03 -0.125475 0.7494 -1.4 1.4 -1.4 1.4 0.006 700

f:id:sheemaa:20140605211503p:plain


まとめ的なものは前回と一緒なので省略。
いい感じの暇つぶしにはなりそう。
わりと時間がかかる処理なので、ベンチマーク代わりにつかったり、CUDAで実装して性能を比較してみたりするにはいいかもしれないですね。