uni memo

juliaでjulia集合のプロットとアニメーション

juliaでjulia集合のプロットとアニメーション

juliaでjulia集合を作成した。ついでにmandelbrot集合も。配列系の演算子(.など)の勉強になった

環境

  • julia: 1.4.1
  • Plots: 1.5.9

mantelbrot集合

zn+1=zn2+cz0=0z_{n+1} = z^2_{n}+c \\\\ z_{0} = 0

として漸化式をとく

cは複素数なので実数a,bより以下のように表せる

c=a+ibc = a + ib

このとき zz \neq \infty となるcの集まりをmantelbrot集合という

zを実数(x,y)からなるz=x+iyz=x+iyとして表現すると以下のように計算できる

zn+1=(xn+iyn)2+c(xn+iyn)2+c=xn2+2ixnynyn2+a+ibz_{n+1} = (x_n+iy_n)^2 + c \\\\ (x_n+iy_n)^2 + c = x_n^2 + 2ix_ny_n -y^2_n + a + ib

zn+1=xn+1+iyn+1z_{n+1} = x_{n+1} + iy_{n+1}より実部と虚部を分ける

xn+1=xn2yn2+ayn+1=2xnyn+bx_{n+1} = x_n^2 - y^2_n + a \\\\ y_{n+1} = 2x_ny_n + b

実装

虚部を使用して計算

  • パッケージと関数
mandelbrot.jl
# install package
using Pkg
Pkg.add("Plots")
using Plots

function mandelbrot(c, n=2000)
    z = zero(c)
    for i = 1:n
        z = z^2 + c
        abs2(z) > 4 && return i
    end
    return n
end
  • plot
mandelbrot-plot.jl
# plot
x = range(-2, 1; length=300)
y = range(-1.5, 1.5; length=300)
c = x' .+ y * im

r = mandelbrot.(c, 300)
heatmap(x, y, log10.(r),
    colorbar=false, size=(300, 300), color=reverse(cgrad(:jet1)))

  • アニメーション
mandelbrot-gif.jl
# 始点
xs = -1.45
ys = 0
t = range(0.01, 4, length=200)
ss = @. exp(-2*t)
anim = @animate for s in ss
    x = range(xs-s, xs+s, length=300)
    y = range(ys-s, ys+s, length=300)
    c = x' .+ y * im
    r = mandelbrot.(c, 300)

    # title
    s = @sprintf "x=%.3f~%.3f, y=%.3f~%.3f" xs-s xs+s ys-s ys+s;
    plot(x, y,label=s, size=(250,250))  
    heatmap!(x, y, log.(r),
        colorbar=false, size=(300, 300), color=reverse(cgrad(:jet1)))
end
gif(anim, fps=40, "mandelbrot.gif")

mandelbrot set

julia集合

mandelbrotと同じ漸化式を用いる。初期値と変化させる値が異なる

c=a+ibc = a + ibの(a,b)を変化させる

zz \neq \inftyとなるznz_nの集まりをjulia集合という

以下の初期値についてx, y, cを代入して求めていく

zn=xn+iyn+cz_n = x_n + iy_n + c

実装

  • パッケージと関数
julia.jl
# install package
using Pkg
Pkg.add("Plots")
using Plots

function julia(c, z, n=2000)
    for i = 1:n
        z = z^2 + c
        abs2(z) > 4 && return i
    end
    return n
end
  • plot
julia-plot.jl
x = range(-1.5, 1.5; length=600)
y = range(-1.25, 1.25; length=400)
a = -0.3
b = -0.63
c = complex(a, b)

r = @.julia(c, complex(x', y), 400)
heatmap(x, y, log10.(r),
    colorbar=false, size=(300, 300), color=:jet1)

  • アニメーション
julia-gif.jl
x = range(-1.5, 1.5; length=300)
y = range(-1.25, 1.25; length=300)
a = -0.7
@time anim = @animate for s in 0.5:-0.0025:0.27
    #a = s
    b = s
    c = complex(a, b)
    r = @.julia(c, complex(x', y), 600)
    s = @sprintf "a=%.3f, b=%.3f" a b;
    plot(x, y,label=s,size=(250,250))
    heatmap!(x, y, log10.(r),
        colorbar=false, size=(300, 300), color=:jet1)
end

gif(anim, "julia_set.gif", fps=40)

julia set

演算子

気になったところのメモ

  • 式の定義の前に
    @.
    をつけると全ての演算子が.形式に置換される

Mathematical Operations and Elementary Functions · The Julia Languageより

Furthermore, “dotted” updating operators like a .+= b (or @. a += b) are parsed as a .= a .+ b, where .= is a fused in-place assignment operation (see the dot syntax documentation).

参考

2024, Built with Gatsby. This site uses Google Analytics.