つれづれなる備忘録

日々の発見をあるがままに綴る

gnuplotによるグラフ作成27~再帰定義

 今回はgnuplot再帰定義の方法とその応用として漸化式やフラクタルのプロットについて紹介する。

atatat.hatenablog.com

1. 再帰定義

再帰定義を行うには、関数を定義する際にその関数自体を含むようにする。例えば階乗N!はf(n)=n*f(n-1)と書けるので以下のように再帰定義することができる。 ただし、n=0のときに値が必要になるので、前回紹介した3項演算子を用いてn=0に対して場合分けを行う。

fac(n)=(n==0)?1:n*fac(n-1)
print fac(1)
>1
print fac(4)
>24

プロットする場合は、xが整数以外も含む可能性があるため、int(x)として引数が必ず整数になるようにする。 また再帰定義の計算は時間がかかるため、サンプル点を減らす工夫をする。

set xrange[0:10]
set samples 11
set logscale y
f(x)=(int(x)==0)?1:int(x)*fac(int(x)-1)
plot f(x)

"階乗のプロット(x=10まで)"
階乗のプロット(x=10まで)

2. 漸化式への応用

 フィボナッチ数列は以下の漸化式で与えられる。

\displaystyle F_{n}=F_{n-1}+F_{n-2}

ただしF0=0, F1=1。 2つの項を設定しないといけないため、3項演算子を2回用いる。

fibo(x)=(int(x)==0)?0:((int(x)==1)?1:fibo(int(x)-1)+fibo(int(x)-2))
unset logscale
set xrange[0:11]
set samples 12
plot fibo(x) with points

x=11までのプロットだが、再帰定義関数を用いると計算に若干時間がかかる。

"フィボナッチ数列のプロット"
フィボナッチ数列のプロット

3. フラクタルへの応用

 マンデルブロ集合のときは初期値z0=0.0と固定していたが、ジュリア集合では初期値はその座標値で与えれる。

\begin{eqnarray} z_{0} &=& X+iY \\ z_{n+1} &=& z_{n}^{2} + A \end{eqnarray}

関数mandelはx,yはそのままにしてzの引数を再帰定義している。計算を打ち切る条件としてzの絶対値が2を超えるか、代入回数が200以上になるという条件を設定している。

mandel(x,y,z,n)=(abs(z)>2.0 || n>=200) ? n: mandel(x,y,z*z+complex(x,y),n+1)
set xrange [-0.5:0.5]
set yrange [-0.5:0.5]
set logscale z
set isosample 50
set hidden3d
set contour
a= -0.37
b= -0.612
splot mandel(a,b,complex(x,y),0) notitle

A=[-0.37,-0.612]の場合のプロットを下に示す。

"ジュリア集合(A=-0.37,-0.612)"
ジュリア集合(A=[-0.37,-0.612])

bをわずかに変化させたA=[-0.37,-0.6]の場合のプロットを以下に示すが、bがわずかに異なるだけでプロットが大きく変化していることがわかる。

"ジュリア集合(A=-0.37,-0.6)"
ジュリア集合(A=[-0.37,-0.6])

4. まとめ

今回はgnuplot再帰定義の方法とその応用として漸化式やフラクタルのプロットについて紹介した。再帰定義は漸化式などそのまま扱えるので便利な反面、繰り返し計算回数が増えると時間がかかるというデメリットがある。その場合はCなど別のプログラムで作成して、プロットのみを行ったほうがよい。