JuliaのパッケージGeometricIntegrators.jlを書いて公開してみた

Julia Advent Calendar 2016の４日目として書いてみます。
著者は、ふだんはFortranを使って物性のシミュレーションなどをしている、Julia初心者です。
Julia練習帳にメモを取りながら練習しています。
アドバイスなどをいただけないかと思い、Julia Advent Calendar 2016に挑戦してみます。

新型コロナウイルスが流行って外出が難しくなった2020年4月にGeometricIntegrators.jlをJulia 1.4に対応させたので、本文章も更新してみました。

パッケージとは

Juliaは簡単かつ強力に機能を拡張できるパッケージ・マネージャPkgを持っています。例えばusing Pkg; Pkg.add("Winston")などで$HOME/.julia/以下にそのパッケージがインストールされます。

パッケージの作りかた

JuliaではGitHubに自分で書いたパッケージを置くことで簡単に公開することができます。パッケージの作りかたと公開の方法はどこかに書いてありますが、他人のパッケージのマネをするのが簡単です。

Geometric integratorとは

Geometric integratorとは常微分方程式の解法（数値積分法）の一つです。調和振動子や振り子や惑星の運動のシミュレーションをするのに便利です。たいてい系のエネルギーが原理的によく保存されます。

GeometricIntegrators.jl

紹介

GeometricIntegrators.jlはJuliaでそのGeometric integrationをするモジュールを提供します。Julia-1.4で使えるはずです。ホームページは https://github.com/t-nissie/GeometricIntegrators.jl です。

セットアップ方法

依存するパッケージなどはないので

using Pkg; Pkg.add(PackageSpec(url="https://github.com/t-nissie/GeometricIntegrators.jl", rev="master"));

とするだけでインストールできます。

使用例

モジュールの中の例

モジュールsrc/GeometricIntegrators.jlの最後のほうに２つの独立な調和振動子の例が書いてあります。Julia-1.4とWinstonをお使いなら、

$ cd src/
$ julia GeometricIntegrators.jl

とするとGeometricIntegrators.epsというプロットが得られます。

１次元調和振動子

test/harmonic_oscillator.jlは１次元調和振動子の例です。harmonic_oscillator.jlを実行するとharmonic_oscillator_??.epsなどのいくつかのプロットが得られます。

２体問題

test/twobody.jlは２体問題の例です。twobody.jlを実行すると、図１のようなtwobody_??.epsのプロットが得られます。

図１２体問題。質量M=2とm=1の質点の軌跡。

テスト

GeometricIntegrators.jlパッケージをクローンして、テストして、消去するだけなら

using Pkg; Pkg.add(PackageSpec(url="https://github.com/t-nissie/GeometricIntegrators.jl", rev="master")); Pkg.test("GeometricIntegrators"); Pkg.rm("GeometricIntegrators")

とすればOKです。

使用方法

座標と運動量の微分をqdotとpdotとして用意してください。両方とも関数の配列です。ハミルトニアンを作って、その正準方程式を考えるのがよいでしょう。あとは初期値とともに例にならって使ってみてください。

Travis CI で Continuous Integration

travis-ci.orgにサインインして.travis.ymlというファイルを適宜用意しておけばgit-pushのたびにテストを実行してくれるように設定できます。Cron Jobsとして設定しておけば定期的かつ自動的に最新のJuliaでパッケージが動くかどうかチェックしてもらえます。さらにcoveralls.ioにサインインして.travis.ymlのafter_success:にデータを送るように設定しておけばcode coverageをチェックすることもできます。travis-ci.orgとcoveralls.ioのどちらもGitHubのアカウントでサインインできます。.travis.ymlなどは次のドキュメントを参考に設定しました。

パッケージを作っていて困ったこと

Julia-0.5での変更にハマりました

Julia-0.4.6:

julia> VERSION
v"0.4.6"

julia> typeof(x -> x[1])
Function

Julia-0.5:

julia> VERSION
v"0.5.0"

julia> typeof(x -> x[1])
##1#2

これはJulia-0.5の新しい仕様みたいです→Julia言語の0.5の変更点

こんなかんじで困っています：

julia> VERSION
v"0.5.0"

julia> function g(f::Array{Function,1})
           x = [1,2,3]
           map(f_each -> f_each(x), f)
       end
g (generic function with 1 method)

julia> g([x -> 1x[1], x -> 2x[2], x -> 3x[3]])
3-element Array{Int64,1}:
  1
  4
  9

julia> g([x -> 1x[1], x -> 2x[2]])
2-element Array{Int64,1}:
  1
  4

julia> g([x -> 1x[1]])
ERROR: MethodError: no method matching g(::Array{##13#14,1})
Closest candidates are:
  g(::Array{Function,1}) at REPL[1]:2

function g(f::Array{Function,1})と引数の型を指定しないで、ただfunction g(f)と書けば動きました。

Julia-1だかの変更にハマりました

where を使って書くことになったみたいです。

function energy(q::Array{T,2},
                p::Array{T,2}, m, k, n::Integer) where {T<:AbstractFloat}

あと、Project.tomlを用意することになったみたいです。パッケージ・マネージャPkgにリンクのある文章を参考にしてください。

他の常微分方程式を解くパッケージ

RungeKutta.jl （ただし、Runge-Kutta法はGeometricではない、もしくは、Symplecticではないです）更新が止まってしまったみたいです。
Sundials.jl
- #MonthOfJulia Day 21: Differential Equations
DifferentialEquations.jl
- Introducing DifferentialEquations.jl
その他JuliaDiffEqにあるもの