AdventCalendar
Julia
ode
JuliaDay 4

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

More than 1 year has passed since last update.

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

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

パッケージとは

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

パッケージの作りかた

JuliaではGitHubに自分で書いたパッケージを置くことで簡単に公開することができます。パッケージの作りかたと公開の方法はFinalizing Your Julia Package: Documentation, Testing, Coverage, and Publishingなどに書いてありますが、他人のパッケージのマネをするのが簡単です。わたくしはRungeKutta.jlのマネをしました。

Geometric integratorとは

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

GeometricIntegrators.jl

紹介

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

セットアップ方法

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

Pkg.clone("git://github.com/t-nissie/GeometricIntegrators.jl.git")

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

使用例

モジュールの中の例

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

$ cd src/
$ julia GeometricIntegrators.jl

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

1次元調和振動子

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

2体問題

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

twobody

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

テスト

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

Pkg.clone("git://github.com/t-nissie/GeometricIntegrators.jl.git"); Pkg.test("GeometricIntegrators"); Pkg.rm("GeometricIntegrators")

とすればOKです。

使用方法

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

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

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)と書けば動きました。

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