Optie研

パソコンで絵や動画を作る方法について

Julia言語と Plots + GR で複素関数のgifアニメーションを作る

1. この記事について

やりたいこと

 最近, 趣味で複素解析の初歩などを勉強しているのですが, ある複素数 $z,w$を

$ \begin{cases} z=x+yi \\ w=u+vi \end{cases} $

とおくと, 複素関数 $w=f(z)$ は, 二変数に二変数が対応する関数

$f:(x,y)\rightarrow(u,v)$

になります。素直にその様子を描こうとすると4次元のグラフが必要ですが, それは不可能です。そのため, 複素関数はそのふるまいを直観的に理解することがやや難しいです。

複素関数の様子を知ろうとすると, 基本的には, 「x-y平面上の点や直線や円などの基本的な図形が, 関数によってu-v平面にどのように写されるか?」を考えることになるのですが, それにしても情報は落ちているので, 様子を把握するには慣れが必要そうな感じがします。

 そこで, 本記事では, いくつかの複素関数に対し, 「x-y平面上で図形をアニメーションさせて, それがu-v平面上でどのように動くか」を見ることで, 複素関数に対する理解を深めていきたいと思います。

やること

 複素関数のプロットとアニメーションということで, いつも通りPythonでも良いのですが, せっかくなので Julia という言語を初めて使ってみることにします。

f:id:Optie_f:20180328233006p:plain:w300
(出展: https://github.com/JuliaLang/julia)

https://julialang.org/

 Julia とは, 数値計算・科学技術計算用に作られた言語で, 計算が非常に高速であることで知られています。

 また, 言語仕様がシンプルで, Python以上にすっきりしたコードを書くことができたり, 非常に強力な標準ライブラリやマクロ機能を備えていたり, Pythonなどを簡単に呼び出せるために他の資産の利用が容易などといった特徴があるようです。もちろん, Jupyter notebookで利用することもできます。

 言語としてはまだ若く, 利用者もPythonに比べてまだまだ多くはないのですが, 今後非常に将来性がありそうな感じがしています。

 ということで, 今回はタスク自体は簡単そうということもあり, この機に習得してみたいと思います。

2. 下準備

 Julia 本体(および Python, jupyter notebook )はインストール済みであるものとします。 Juliaのインストールについては以下などを参考。

Julia v0.6.1のインストールとJupyter Notebookで使うまで - Qiita

パッケージのインストール

Pythonなどと同様, ターミナル上などでjuliaを起動するとjuliaのREPLが開くので, そこでPkg.add() で利用したいパッケージ名を入力します。

その後, Pkg.update() を行い, アップデートなどを確認するのが一般的なパッケージ導入の流れです。

今回は以下のパッケージを利用しました。

Pkg.add("Ijulia")
Pkg.add("Plots")
Pkg.add("GR")
Pkg.add("LaTeXStrings")
Pkg.add("ProgressMeter")
Pkg.update()

それぞれ以下に解説します。

IJulia

Jupyterでjuliaを利用するためのカーネルPkg.add("Ijulia") とするだけで使えるようになる。

Plots

Juliaで標準的な, 描画用ライブラリのフロントエンド(インターフェース)。 pyplotやGRなどといったバックエンドと組み合わせて使う。

GR

JuliaでしばしばPlotsのバックエンドとして用いられる描画用パッケージ。C/C++で書かれているらしく, かなり速い。

LaTeXStrings

L"z=x^2" など書くことで, Plotsなどに渡した文字列が LaTeX文字列と認識され, $z=x^2$ というように表示されるようになる。(しかしあまり綺麗に出なかったので今回使わなかった)

ProgressMeter

ループなどを回すときに, 進捗メーターを表示させるためのパッケージ。

3. プロット & アニメーション

importとヘルパー関数など

Juliaでのモジュールのインポートは以下のようにします。

using Plots
import ProgressMeter
using LaTeXStrings

importusing がありますが, 以下のような違いがあるようです。

 例えば, import Plots とすると, Plots モジュール内の plot() 関数を, Plots.plot() と書くことで呼び出せるようになります。また, import Plots.plot などと書けば,Plots. を省略して plot() と書くことができます。これはPythonのimportと同じですね。

 using Plots とする場合, モジュール自体に加えて Plots モジュールの名前空間ごと Main に持ってくるため, Plots.plot()関数であればPlots. を省略して plot() と書くことができます。ただし, 関数名などが衝突する恐れがあると警告が出るようです。

以下は指定ディレクトリにあるgifファイルを読み込んで jupyter 上で表示するための関数です。 先にgifを生成してから, この関数を使用して表示します。

# open ... do ... endはdoの中でエラーが起きても必ずファイルを閉じる
displayfile(mime, file) = open(file) do stream
    base64 = base64encode(stream)
    display("text/html", """<img src="data:$mime;base64,$base64">""")
end

Plotsをimportして, 以下のように書けばGRバックエンドが立ち上がり, 後は plot() などの関数を実行することでプロットすることができます。

gr()

しかし, 最新版の GR を使うと, なぜか出力の画像サイズが指定の6倍になってしまうという仕様があります。そこでENV["PLOTS_TEST"] = "true" と書くことでこれを無効にします。

ENV["PLOTS_TEST"] = "true"

例1. $z=w^2$

まずは練習も兼ねて, $z=w^2$ をプロットします。

gr(
    legend=false,
    titlefont=Plots.font("sans-serif", 12),
    legendfont=Plots.font("sans-serif", 8),
    guidefont=Plots.font("sans-serif", 10),
    tickfont=Plots.font("sans-serif", 8),
    markersize=4, markerstrokewidth = 1, 
    markercolor= "#FF445B",markerstrokecolor="#4A4A4A",
    linewidth=2
)
n = 100  # フレーム数
prog = ProgressMeter.Progress(n,1)


anim = @animate for i in linspace(-10,10,n)
    k = i

    x1 = linspace(-50,50,201)
    y1 = k
    z1 = x1 + y1*im

    x2 = k
    y2 = linspace(-50,50,201)
    z2 = x2 + y2*im
    
    z3 =  k + k*im

    w1 = z1.^2
    w2 = z2.^2
    w3 = z3.^2

    p1 = plot(real(z1), imag(z1), linecolor="#F5A623", lims=(-20,20), title="x-y")
    plot!(p1, real(z2), imag(z2), linecolor="#4A90E2")
    vline!(p1,[0],line=(:black, 1)) # y軸
    hline!(p1,[0],line=(:black, 1)) # x軸
    scatter!(p1, (real(z3),imag(z3)))


    p2 = plot(real(w1),imag(w1), linecolor="#F5A623",  lims=(-200,200), title="u-v")
    plot!(p2, real(w2),imag(w2), linecolor="#4A90E2")
    vline!(p2,[0],line=(:black, 1))
    hline!(p2,[0],line=(:black, 1))
    scatter!(p2, (real(w3),imag(w3)))
    
    plot(p1,p2,size=(800,350), plot_title=L"w=z^2")
    
    ProgressMeter.next!(prog)
end


dir = "/"
fps = 30
filename = "01_fps$fps.gif"

gif(anim, filename, fps = fps)

displayfile(dir,filename)

f:id:Optie_f:20180328225511g:plain

これは有名な例ですが, $z=w^2$ は,

$ \begin{cases} z=x+yi \\ w=u+vi \end{cases} $

とおくと, $w = z^2 = x^2-y^2+2xyi$ より,

$ \begin{cases} u=x^2-y^2 \\ v=2xy \end{cases} $

となります。$x,y$ をそれぞれ固定して考えると, u-v平面では横倒しの二次関数の形となるため, プロットのように直線が二次関数に写されます。

また, x-y 平面で直交する二直線は, u-v 平面の二次関数に写されても互いに直交することが確認できます。

関数化

これらの処理を関数化してみます。任意の複素関数と, 定数kに応じて変化する図形群を与えられるようにします。

# 複素関数 f, 定数kを受け取り複素数列(点群)を返す shape の配列 shapes
function makeanim(f, shapes; n=100, kmin=-10, kmax=10, zsize=10, wsize=10, unitcircle=true)
    prog = ProgressMeter.Progress(n,1)
    
    anim = @animate for i in linspace(kmin,kmax,n)
        
        p1 = plot(lims=(-zsize, zsize), title="x-y")
        p2 = plot(lims=(-wsize, wsize), title="u-v")
        vline!(p1,[0],line=(:black, 1))
        hline!(p1,[0],line=(:black, 1))
        vline!(p2,[0],line=(:black, 1))
        hline!(p2,[0],line=(:black, 1))    
        
        if unitcircle
            θ = linspace(0,2π,60)
            plot!(p1, sin.(θ), cos.(θ), line=(:black, 1))
            plot!(p2, sin.(θ), cos.(θ), line=(:black, 1))
        end
            
        
        k = i
        
        for shape in shapes
            z = shape(k)
            w = f(z)
            
            if length(z)>1  # 点ならscatterを使う
                plot!(p1, real(z), imag(z))
            else
                scatter!(p1, (real(z),imag(z)))
            end
            
            if length(w)>1
                plot!(p2, real(w),imag(w))
            else
                scatter!(p2, (real(w),imag(w)))
            end
            
        end
        
        plot(p1,p2,size=(800,400))
        ProgressMeter.next!(prog)
        
    end
    
    return anim
    
end

図形は以下のようにします。

function yline(k)
    x = linspace(-50,50,1001)
    y = k
    z = x + y*im
end

function xline(k)
    x = k
    y = linspace(-50,50,1001)
    z = x + y*im
end

function xylinecrosspoint(k)
    x = k
    y = k
    z = x + y*im
end

function circle(k)
    r = k
    θ = linspace(-π,π,10000)
    z = r*(cos.(θ) + im*sin.(θ))
end

function circleline(k)
    r = linspace(0,100,10001)
    θ = 2k
    z = r*(cos.(θ) + im*sin.(θ))
end

function circlepoint(k)
    r = k
    z = r*(cos.(2k) + im*sin.(2k))
end

これらを以下のように使います。

function f1(z)
    z.^2
end
shapes=[yline, xline, xylinecrosspoint]
anim = makeanim(f1, shapes, zsize=20, wsize=200, unitcircle=false)

dir = "/"
fps = 30
filename = "01_fps$fps.gif"

gif(anim, filename, fps = fps)

displayfile(dir,filename)

出力は省きますが, 色を除いて上記と同じ結果となります。

例2. $w=\frac{1}{z}$

function f2(z)
    1./z
end
shapes1=[yline, xline, xylinecrosspoint]
shapes2=[circle, circleline, circlepoint]
dir = "/"
fps = 30
filename1 = "02_1_fps$fps.gif"
filename2 = "02_2_fps$fps.gif"
anim1 = makeanim(f2, shapes1, n=150, kmin=-2, kmax=2, zsize=2, wsize=2)
gif(anim1, filename1, fps = fps)
displayfile(dir,filename1)

f:id:Optie_f:20180328225652g:plain

anim2 = makeanim(f2, shapes2, n=250, kmin=2π, kmax=0, zsize=2π, wsize=2π)
gif(anim2, filename2, fps = fps)
displayfile(dir,filename2)

f:id:Optie_f:20180328225711g:plain

実軸対称 + 単位円に関して対称なやつです。直線の方は原点付近でリーマン球面をぐるりと回り込んでいますね。

例3. $w= \sin z$

function f3(z)
    sin.(z)
end
dir = "/"
fps = 30
filename1 = "03_1_fps$fps.gif"
filename2 = "03_2_fps$fps.gif"
anim1 = makeanim(f3, shapes1, n=150, kmin=-2, kmax=2, zsize=2, wsize=2)
gif(anim1, filename1, fps = fps)
displayfile(dir,filename1)

f:id:Optie_f:20180328225716g:plain

三角関数です。この辺は円を入れると異常な動きになります(よくわからない)。

anim2 = makeanim(f3, shapes2, n=500, kmin=π, kmax=0, zsize=2π, wsize=2π)
gif(anim2, filename2, fps = fps)
displayfile(dir,filename2)

f:id:Optie_f:20180328231606g:plain

例4. $w=\tan z$

function f4(z)
    tan.(z)
end
dir = "/"
fps = 30
filename1 = "04_1_fps$fps.gif"
filename2 = "04_2_fps$fps.gif"
anim1 = makeanim(f4, shapes1, n=150, kmin=-2, kmax=2, zsize=2, wsize=2)
gif(anim1, filename1, fps = fps)
displayfile(dir,filename1)

f:id:Optie_f:20180328231530g:plain

anim2 = makeanim(f4, shapes2, n=500, kmin=π, kmax=0, zsize=2π, wsize=2π)
gif(anim2, filename2, fps = fps)
displayfile(dir,filename2)

f:id:Optie_f:20180328231357g:plain

例5. $w=e^z$

function f5(z)
    e.^z
end
dir = "/"
fps = 30
filename1 = "05_1_fps$fps.gif"
filename2 = "05_2_fps$fps.gif"
anim1 = makeanim(f5, shapes1, n=150, kmin=-2, kmax=2, zsize=2, wsize=2)
gif(anim1, filename1, fps = fps)
displayfile(dir,filename1)

f:id:Optie_f:20180328231321g:plain

直線が円に写ります。x-y平面の紫の線は, $2\pi$ 上下にずらしてもu-v平面で同じ線に写ります。

anim2 = makeanim(f5, shapes2, n=500, kmin=π, kmax=0, zsize=2π, wsize=4π)
gif(anim2, filename2, fps = fps)
displayfile(dir,filename2)

f:id:Optie_f:20180328231145g:plain

例6. $w = \log z$

function f6(z)
    log.(z)
end
dir = "/"
fps = 30
filename1 = "06_1_fps$fps.gif"
filename2 = "06_2_fps$fps.gif"
anim1 = makeanim(f6, shapes1, n=150, kmin=-2, kmax=2, zsize=2, wsize=2)
gif(anim1, filename1, fps = fps)
displayfile(dir,filename1)

f:id:Optie_f:20180328231942g:plain

$w=e^z$の逆で, 円が直線に写ります。u-v平面では黄色い線は $2\pi$ おきに無限に繰り返されます。

例7. $w= \mathrm{sin{h}} z$

function f7(z)
    sinh.(z)
end
dir = "/"
fps = 30
filename1 = "07_1_fps$fps.gif"
filename2 = "07_2_fps$fps.gif"
anim1 = makeanim(f7, shapes1, n=150, kmin=-2, kmax=2, zsize=2, wsize=2)
gif(anim1, filename1, fps = fps)
displayfile(dir,filename1)

f:id:Optie_f:20180328231906g:plain

anim2 = makeanim(f7, shapes2, n=200, kmin=π, kmax=0, zsize=π, wsize=π)
gif(anim2, filename2, fps = fps)
displayfile(dir,filename2)

f:id:Optie_f:20180328231821g:plain

$ \mathrm{sin{h}} z$ の定義より, $w=\sin z$の$z$に$i$をかけて回転させたものになります。

例8. $w=z^\frac{1}{3}$

function f8(z)
    z.^(1/3)
end
dir = "/"
fps = 30
filename1 = "08_1_fps$fps.gif"
filename2 = "08_2_fps$fps.gif"
anim2 = makeanim(f8, shapes2, n=200, kmin=π, kmax=0, zsize=π, wsize=π)
gif(anim2, filename2, fps = fps)
displayfile(dir,filename2)

f:id:Optie_f:20180328232240g:plain

x-y平面の円全体がu-v平面の$\frac{1}{3}$の扇状の領域に対応します。

4. まとめ

 複素関数のふるまいについて理解が深まりました。また, Juliaの書き方とプロットのやり方がわかりました。

Python+matplotlibよりかなり高速で使いやすいので, 今後も使っていきたいです。

5. 参考文献