Julia言語(v1.0)でJulia集合を描く

## はじめに

この記事に書かれているコードは"[決定版？「Juliaでジュリア集合を描く」](https://qiita.com/tenfu2tea/items/6bcbdd7586ea070cc25d)"のコードを参考に安定版であるJulia v1.0で動くようにしたものです。 

## Julia集合とは

Julia集合(Julia set)とは

>A set of complex numbers which do not converge to any limit when a given mapping is repeatedly applied to them.（オックスフォード新英英辞典）

とのことです。
ここでは最もポピュラーなJulia集合に焦点を当てます。それは、漸化式

```math
z_{n+1}=z_n^2+c

```

において、$c$を固定したとき、$n\rightarrow\infty$においてこの漸化式が無限大に発散しないような初期値 $z_0$ の集合をさします。

ある初期値$z_0$について漸化式が無限大に発散するかどうかを判断するには、$n$を増えるたびに$|z_n|$が2より大きいかどうかを見てやればよいそうです。どこかの$n$で$|z_n|>2$になれば漸化式は発散することが確定で、それ以上の$n$について計算する必要はありません。

前記事によれば計算アルゴリズムは
>
1. 適当な複素数$c=a+bi$を決めて
2. 適当な複素数$z_0$に対して
3. $z_n$が発散する ($|z_n|>2$)まで$z_{n+1}=z_n^2+c$を繰り返し
4. 発散した時の$n$を集める
5. ステップ2から4を十分多く繰り返す
6. 集めた$n$をプロットする

です。フローチャートにおこすと、

![Julia_Set_flow_R.png](https://qiita-image-store.s3.amazonaws.com/0/253018/f1d89b4f-c2df-b2df-34db-f655713387dc.png)


みたいな感じでしょうか。
注１：Julia言語において配列のインデックスは1から始まります。そのためフローチャート図では初期値を$z_0$ではなく$z_1$と表記しています。
注２：実際の計算では$1/(10+n)$を収集しています。

## 計算
6つの複素数$c$についてそれぞれJulia集合を計算し、図示します。
バックエンドにはGR,フロントエンドにはPlotsを使用しています。

```julia
versioninfo()
```

Julia Version 1.0.0
Commit 5d4eaca0c9 (2018-08-08 20:58 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-6500U CPU @ 2.50GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)

```julia:julia_set.jl
using Plots
gr()

function compute_julia_set(zc,ju,zs,nmax)
    ju[:].=1.0/(10.0+nmax)
    
    for i=eachindex(zs)
        zz=zs[i]      
        for n=1:nmax          
            if abs2(zz)>4.0
                ju[i]=1.0/(10.0+n)
                break
            end
            zz=zz*zz+zc
        end
    end
end

function do_compute_julia_set()
    zs=[ComplexF64(x,y) for x in range(-1.5,stop=1.5,length=500), y in range(-1.0,stop=1.0,length=500)]
    xs=real(zs)
    ys=imag(zs)
    
    ju =zeros(Float64,size(zs))
    nmax=500
    result=[]
    
    for zc in [
            ComplexF64(1.0-(1.0+sqrt(5.0))/2.0,0.0),
            ComplexF64( 0.28500,0.0000),
            ComplexF64( 0.28500,0.0100),
            ComplexF64(-0.70176,0.38421),
            ComplexF64(-0.83500,0.2321),
            ComplexF64(-0.80000,0.1560),
         ]
        
        @time compute_julia_set(zc,ju,zs,nmax)
        
        p=plot(legend=false,ticks=false)
        plot!(title="c = $zc",titlefontsize=8)
        plot!(transpose(ju),seriestype=:contour,fill=true)
        push!(result,p)
    end

    plot(result...)
    #savefig("julia_set_using_julialang")
end

do_compute_julia_set()
```
## 出力結果

 0.244310 seconds (2 allocations: 80 bytes)
  0.010541 seconds (2 allocations: 80 bytes)
  0.021119 seconds (2 allocations: 80 bytes)
  0.011475 seconds (2 allocations: 80 bytes)
  0.011199 seconds (2 allocations: 80 bytes)
  0.050636 seconds (2 allocations: 80 bytes)

![julia_set_using_julialang.png](https://qiita-image-store.s3.amazonaws.com/0/253018/4130d6df-a612-ff3b-a054-2c544b387e82.png)

色が明るいほど$1/(10+n)$が大きい（＝発散速度が大きい）ことを示します。きれいなJulia集合が描けていると思います。

## 変更点
* `linspace`から`range`へ変更。
* ローカル変数の`const`宣言はサポートされなくなりました。そのため`zs`を宣言したとき`const`をつけていません
* 描画ライブラリをMatplotlibからPlotsへ変更（単に自分が使い慣れてるので）。
* 実数軸を横軸、虚数軸を縦軸にするために、transposeした行列をplot関数に渡しています。
* `...`演算子について、ここでは`result`が1個の引数ではなくコレクションであることを`plot!`関数に伝えるために使用しています(splattingという機能)。https://docs.julialang.org/en/v1/manual/faq/#What-does-the-...-operator-do?-1


## わからないところ
今は`ticks=false`で隠してありますが、軸刻みラベルを表示させると配列インデックス(100 200 300 400 500)で表示されます。これを対応する$x,y$(-1.5 -0.9 -0.3 0.3 0.9 1.5)で置き換えるにはどう書けばいいのでしょう？

（20181105追記）解決しました。`ticks=(0:100:500, -1.5:0.6:1.5)`と書けばよいです。@goropikariさんアドバイスありがとうございます。

## (20181105追記)改訂版

軸刻みラベルを指定する方法が分かったので表示させるようにしました。ついでに、軸刻みラベルの最小値と最大値の指定には作っただけで使ってなかった`xs`,`ys`を使ってみました。

```julia
using Plots
gr()

function compute_julia_set(zc,ju,zs,nmax)
    
    ju[:].=1.0/(10.0+nmax)
    
    for i=eachindex(zs)
        zz=zs[i]
        
        for n=1:nmax          
            if abs2(zz)>4.0
                ju[i]=1.0/(10.0+n)
                break
            end
            zz=zz*zz+zc
        end
    end
end

function do_compute_julia_set()
    zs=[ComplexF64(x,y) for x in range(-1.5,stop=1.5,length=500), y in range(-1.0,stop=1.0,length=500)] 
    xs=real(zs)
    ys=imag(zs)
    
    ju =zeros(Float64,size(zs))
    nmax=500
    result=[]
    
    for zc in [
            ComplexF64(1.0-(1.0+sqrt(5.0))/2.0,0.0),
            ComplexF64( 0.28500,0.0000),
            ComplexF64( 0.28500,0.0100),
            ComplexF64(-0.70176,0.38421),
            ComplexF64(-0.83500,0.2321),
            ComplexF64(-0.80000,0.1560),
         ]
        
        @time compute_julia_set(zc,ju,zs,nmax)
        
        p=plot(legend=false,
               #ticks=false
               xticks=(0:100:500, xs[1]:0.6:xs[end]),
               yticks=(0:100:500,ys[1]:0.4:ys[end]))
        plot!(title="c = $zc",titlefontsize=8)
        plot!(transpose(ju),seriestype=:contour,fill=true)
        push!(result,p)
    end

    plot(result...)
    #savefig("julia_set_using_julialang")
end

do_compute_julia_set()
```

![julia_set_using_julialang.png](https://qiita-image-store.s3.amazonaws.com/0/253018/6e4b2fec-3149-29c7-8378-4d56fe61054f.png)



## 参考文献
* [@tenfu2tea "決定版？「Juliaでジュリア集合を描く」"](https://qiita.com/tenfu2tea/items/6bcbdd7586ea070cc25d)
* [@anmitsu "Juliaでジュリア集合を描く"](https://qiita.com/anmitsu/items/0f85a6427b4224a2d56a)