WANtaroHP (f90: Mandelbrot)



Mandelbrot 集合の定義

Mandelbrot集合は,以下の漸化式で定義される複素数列 $\{z_n\}$ が,$n\rightarrow\infty$ の極限において無限大に発散しないという条件を満たす複素数 $c$ が作る集合として定義される.

\begin{equation} z_{n+1}=z_n{}^2+c \qquad\qquad z_0=0 \end{equation}

複素数を用いないで扱うには, $z$ および $c$ を実数部と虚数部にわけて考える. そのためには, $z=x+i\cdot y$$c=a+i\cdot b$ として,以下の式を用いればよい.

\begin{equation} \begin{cases} x_{n+1}=x_n{}^2-y_n{}^2+a \\ y_{n+1}=2 x_n y_n+b \end{cases} \end{equation}

プログラムでの考え方

  • 計算領域の分割数,計算領域,最終的に作成する画像ファイル名はコマンドライン引数で指定する.
  • 作図用の座標として, $x$ (実数部)を横軸, $y$ (虚数部)を縦軸とし,発散確定までの繰り返し数の常用対数を$z$値(着色ランクあるいは3次元プロット用の $z$ 座標)として扱う.
  • 作図は GMT で行うが,10MB を越える eps が作成され取り扱いにくいので,ImageMagick によりこれを jpg 画像に変換する.
  • これらの処理は Fortran プログラムの中で一括して行う.( call system('bash _gmt.sh') の活用)
  • 発散確定までの繰り返し数の最大値は 10 万回とし,これを越えたら「収束するであろう」,すなわち Mandelbrot 集合に属する点であるとみなす.
  • 数値計算の中で,複素数列の絶対値が 2 を越えた時点で発散すると判断する.この時までの繰り返し数に応じて着色を行う.
  • 実際の着色は,繰り返し数の常用対数( $z$ 値)を用い,この値に GMT のカラーパレットを定義して適用する.
  • GMT には陰線処理機能はないが,3次元プロットでは,奥側から同一 $y$ 座標上の $(x,z)$ 平面の多角形を塗りつぶし,輪郭を描画することにより,それらしく見せている.


プログラム

Programs

FilenameDescription
a_f90gmt.txt実行用シェルスクリプト
f90_mandelbrot_gmt.f90Mandelbrot集合2次元描画
f90_mandelbrot_gmt3d.f90Mandelbrot集合3次元描画

プロット実行用スクリプト

スクリプトでは,コンパイル後,2次元プロット用に GMT での着色指定のための cpt ファイルを作り, Mandelbrot 集合の計算と GMT を実行する Fortran プログラムを実行しています.

gfortran -o f90_mandelbrot_gmt f90_mandelbrot_gmt.f90
gfortran -o f90_mandelbrot_gmt3d f90_mandelbrot_gmt3d.f90

# cpt-file for 2d
cat < _col.cpt
  0   0   0   0  2.5 255 255 255
2.5 255 255 255  5.0 255 255 255
F   0   0   0
B 255 255 255
EOT
./f90_mandelbrot_gmt 800 800 -0.75 0 1.5 fig_man_a800.jpg
./f90_mandelbrot_gmt 800 800 -0.743 0.1145 0.003 fig_man_b800.jpg
./f90_mandelbrot_gmt3d 200 200 -0.75 0 1.5 fig_man_3d_a.jpg
./f90_mandelbrot_gmt3d 200 200 -0.743 0.1145 0.003 fig_man_3d_b.jpg

Fortran プログラム実行用コマンド

./f90_mandelbrot_gmt jw ih x0 y0 rr filename   # 2次元プロット
./f90_mandelbrot_gmt3d jw ih x0 y0 rr filename # 3次元プロット
jw 作図の水平方向分割数(整数値)
ih 作図の鉛直方向分割数(jwと同じ整数値)
x0 作図の中心座標のx値(実数)
y0 作図の中心座標のy値(実数)
rr 作図半径(実数)
filename出力画像ファイル名(jpg)


作図事例

(x0=-0.75, y0=0, rr=1.5) (x0=-0.743, y0=0.1145, rr=0.003)
png png
fig_man_a800.jpg fig_man_b800.jpg
png png
fig_man_3d_a.jpg fig_man_3d_b.jpg


inserted by FC2 system