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$値(着色ランクあるいは３次元プロット用の $z$ 座標)として扱う．
- 作図は GMT で行うが，10MB を越える eps が作成され取り扱いにくいので，ImageMagick によりこれを jpg 画像に変換する．
- これらの処理は Fortran プログラムの中で一括して行う．( call system('xxx') の活用)
- 発散確定までの繰り返し数の最大値は 10 万回とし，これを越えたら「収束するであろう」，すなわち Mandelbrot 集合に属する点であるとみなす．
- 数値計算の中で，複素数列の絶対値が 2 を越えた時点で発散すると判断する．この時までの繰り返し数に応じて着色を行う．
- 実際の着色は，繰り返し数の常用対数( $z$ 値)を用い，この値に GMT のカラーパレットを定義して適用する．
- GMT には陰線処理機能はないが，３次元プロットでは，奥側から同一 $y$ 座標上の $(x,z)$ 平面の多角形を塗りつぶし，輪郭を描画することにより，それらしく見せている．
プログラム
Programs
| Filename | Description |
|---|---|
| f90_mandelbrot_gmt.txt | 2-dim. grey color plot program |
| f90_mandelbrot_gmt3d.txt | 3-dim. plot program |
２次元グレースケールプロット実行用スクリプト
スクリプトでは，コンパイル後，GMT での着色指定のための cpt ファイルを作り， Mandelbrot 集合の計算と GMT を実行する Fortran プログラムを実行しています．
gfortran -o f90_mandelbrot_gmt f90_mandelbrot_gmt.f90 # ******************************************************** # Make color pallet file # z1 R1 G1 B1 z2 R1 G2 B2 (z1/z2: range of value) #********************************************************* echo "0 0 0 0 2.5 255 255 255" > _col.cpt echo "2.5 255 255 255 5 255 255 255" >> _col.cpt echo "F 0 0 0" >> _col.cpt echo "B 255 255 255" >> _col.cpt ./f90_mandelbrot_gmt 200 200 -0.75 0 1.5 fig_man_a0200.jpg ./f90_mandelbrot_gmt 800 800 -0.75 0 1.5 fig_man_a0800.jpg echo "0 0 0 0 2.5 125 125 125" > _col.cpt echo "2.5 125 125 125 5 255 255 255" >> _col.cpt echo "F 0 0 0" >> _col.cpt echo "B 255 255 255" >> _col.cpt ./f90_mandelbrot_gmt 200 200 -0.743 0.1145 0.003 fig_man_b0200.jpg ./f90_mandelbrot_gmt 800 800 -0.743 0.1145 0.003 fig_man_b0800.jpg
上記スクリプトでは，同じ作図領域にたいし，領域分割数を 200 とした場合と 800 とした場合について計算・作図を行っています．
GMT 用 cpt ファイル
- プログラムでは発散が確定するまでの繰り返し回数の常用対数を記録させており，その値に応じて黒から白に色づけしています．
- 発散が確定するまでの繰り返し回数が1であれば log10(1)=0 となり，これには黒が当てはめられます．
- 発散確定までの繰り返し数の最大値は 10 万回としていますので，10 万回で発散が確定しなければ収束するとみなされ，Mandelbrot 集合とみなされます．
- Mandelbrot 集合に属する点には黒色 (0 0 0) が当てはめられます．
- Mandelbrot 集合に属する点と，発散までの繰り返し数 1 の場合と色がダブりますが，通常 Mandelbrot 集合の近傍は収束し辛く，この色分けの考え方では白に近い色になるので，模様作成上の問題はありません．
fig_man_a0800.jpg 用の GMT 用 cpt ファイル
ここでは，発散確定までの回数 1 から発散確定までの繰り返し回数 316 まで，黒から白に向かい順に色分けされます． 発散確定までの繰り返し回数が 316 より大きい場合は，白で着色されます．
# ******************************************************** # Make color pallet file # z1 R1 G1 B1 z2 R1 G2 B2 (z1/z2: range of value) #********************************************************* echo "0 0 0 0 2.5 255 255 255" > _col.cpt echo "2.5 255 255 255 5 255 255 255" >> _col.cpt echo "F 0 0 0" >> _col.cpt echo "B 255 255 255" >> _col.cpt
fig_man_b0800.jpg 用の GMT 用 cpt ファイル
ここでは，発散確定までの回数 1 から繰り返し上限の 10 万まで，黒から白に向かい順に色分けされます．
echo "0 0 0 0 2.5 125 125 125" > _col.cpt echo "2.5 125 125 125 5 255 255 255" >> _col.cpt echo "F 0 0 0" >> _col.cpt echo "B 255 255 255" >> _col.cpt
Fortran プログラム実行用コマンド
./f90_mandelbrot_gmt jw ih x0 y0 rr filename
| jw | 作図の水平方向分割数(整数値) |
| ih | 作図の鉛直方向分割数(jwと同じ整数値) |
| x0 | 作図の中心座標のx値(実数) |
| y0 | 作図の中心座標のy値(実数) |
| rr | 作図半径(実数) |
| filename | 出力画像ファイル名(jpg) |
Fortran プログラムが吐き出す GMT と ImageMagick 実行用コマンド
- 下記リストのうち，始めの 5 行が GMT 用コマンド，最終行が ImageMagick のコマンドです．
- ImageMagick では，convert コマンドにより GMT で作成された eps ファイルを jpg に変換しています．
- この事例では GMT により _gmt.eps が作成され，最終的に fig_man_b0800.jpg ができあがります．
range=-0.7460000E+00/-0.7400000E+00/0.1115000E+00/0.1175000E+00 scale=0.8000000E+01/0.8000000E+01 inp_dat=_out.txt fig=_gmt.eps psxy $inp_dat -R$range -JX$scale -SC0.01 -C_col.cpt -P > $fig convert -trim -density 300 $fig fig_man_b0800.jpg
３次元プロット実行用スクリプト
３次元プロット実行用スクリプト
コマンドライン引数の意味は，２次元の場合と同じです．
gfortran -o f90_mandelbrot_gmt3d f90_mandelbrot_gmt3d.f90 ./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プログラムが吐き出す GMT 作図用スクリプト
range=-0.7460000E+00/-0.7400000E+00/0.1115000E+00/0.1175000E+00/0/7 scale=10/10 inp_dat=_out.txt fig=_gmt.eps psxyz $inp_dat -R$range -JX$scale -JZ5 -G255 -W1 -E150/30 -m -P > $fig convert -trim -density 300 $fig fig_man_3d_b.jpg
作図事例
| fig_man_a0800.jpg | fig_man_b0800.jpg |
| (x0=-0.75, y0=0, rr=1.5) | (x0=-0.743, y0=0.1145, rr=0.003) |
 |
 |
| fig_man_3d_a.jpg | fig_man_3d_b.jpg |
| (x0=-0.75, y0=0, rr=1.5) | (x0=-0.743, y0=0.1145, rr=0.003) |
 |
 |
