Thread by @genkuroki

NASAでどのようにJuliaを使っているか
宇宙船の分離のJuliaによるシミュレーション

多分 #Julia言語 ユーザーはNASAでどのようなパッケージを使っているのか気になると思います。動画のスクショを見て下さい。

https://twitter.com/JuliaComputing/status/1370358602788048901

#Julia言語 NASAでの宇宙船の分離シミュレーションで使っているパッケージ続き

添付が画像③によればMATLABでやっていた場合よりも15000倍高速化されたらしい。

DifferentialEquations.jl で微分方程式を解いて、Optim.jl で最適化を実行しているらしい。

Pluto.jlを使っていることもわかる。

#Julia言語 NASA では Pluto.jl 上で宇宙船の分離のシミュレーションを行っているらしいのですが、スライダーを設置したり、MonteCarloMeasurements.jl で摂動をかけた場合の結果の違いをシミュレーションのコードを変更せずに表示させているようです。

#Julia言語 では、摂動をかけた場合の結果の違いを後で見たくなることに気を使わずに書かれたシミュレーションのコードを一切変更せずに、MonteCarloMeasurements.jl を組み合わせることによって摂動をかけた場合の結果を表示させることができます！

初めて経験したときにはまるで魔法に感じられる。

#Julia言語 Juliaは確かに速いのですが、その点だけを強調するのはひどくミスリーディングで、1つ前に書いたように、自分が書いたコードを変更せずに、他人が書いたパッケージの機能を適用することによって、コードを複雑にすることなく、複雑な仕事を簡単にできてしまう点も強調されるべきです。

#Julia言語 似たようなシンプルな例を Van der Pol 方程式について作って来ました。

NASAでも使っている Pluto.jl のノートブックのソースコードを

gist.github.com/genkuroki/9013…

で公開しておきました。添付画像のようなことをできます。

#Julia言語 ソースコードも添付画像の形式で貼り付けておきます。

添付画像1に微分方程式を数値的に解くためのシンプレクティック・オイラー法のコードがあります。

その部分には、初期条件を摂動したときの効果を計算するためのコードが一切含まれていないことを確認して下さい！

#Julia言語 実は私は初期条件を摂動したときの結果の違いをプロットするためのコードを何も書いていません。

単に、using MonteCarloMeasurements して、初期条件を

u0 = u00 ± ε

と与えているだけです。± ε の部分を削ると、初期条件を摂動した場合の結果の違いは表示されなくなります。

#Julia言語 微分方程式を解くコードも自分で書かずに、NASAと同じように DifferentialEquations.jl にまかせようとかと最初思ったのですが、それだと自分で書くコードが「パラメーターを与えて結果を表示するだけ」になってしまうので、シンプレクティック・オイラー法の部分は自分で書きました(笑)

#Julia言語 さらに微分方程式を数式で表示するためのコードも追加したのに、添付画像のようなことをするために自分が書いたコードは60行程度でしかありません。

Van der Pol方程式のリミットサイクルがしっかり見えています。

Pluto notebook のソースコード↓
gist.github.com/genkuroki/9013…

#Julia言語 NASAでもJuliaを使っているのだから、以下のリンク先のような仕事でもJuliaを使うと便利なのではないかと思いました。

https://twitter.com/madnoda/status/1366886071288492034

#Julia言語 上での Van der Pol 方程式の例では使わなかったのですが、DifferentialEquations.jl には多数の優れた数値解法がすでに実装されており、ものすごく沢山の機能を持っており、DifferentialEquations.jl 以外のパッケージとも組み合わせて使えるように書かれています。

#Julia言語 誰か、太陽系内での宇宙船の軌道のシミュレーションを DifferentialEquations.jl でやって見せて欲しいです。

各種定数の取得には

github.com/cadojo/Unitful…

を使えます。

#数楽 添付画像では最大εの大きさで初期条件に摂動がかけられた場合の様子もプロットされているのですが、リミットサイクルに引き込まれて行っているので、初期条件の違いの影響が時間とともに減少している様子も見えます。

#Julia言語 常微分方程式の数値解法にも種類がものすごく沢山あって、精度と計算速度に間にはトレードオフがあります。自分の目的に合わせて適切な数値解法を選択する必要がある。

数値解法を自分の手でハードコーディングしてしまうとそういう試行錯誤をすることが恐ろしく大変になります。

#Julia言語 解法の選択でも試行錯誤したい場合には、DifferentialEquations.jl のような超強力超巨大パッケージを利用した方が得だと思う。

新しめの数値解法を自分で実装して試すことはかなりの手間が必要な大変な作業になりそうです。

#Julia言語 添付画像は私が書いたシンプレクティック・オイラー法のコードです。

このコードのどこにもFloat64のような浮動小数点数の型にあたる情報が何も書かれていないことに注目！

このような書き方でもJuliaは高速に計算してくれます(今回はpush!を使っているので計算速度に気を使っていない)。

#Julia言語 仮に添付画像のように Float64 (倍精度)のような型の指定をしないで済ますのではなく、Float64を指定してしまっていると、さまざまな害が起こります。

* Float32やFloat128では使えなくなる。

* 摂動した結果を表示するためのパッケージと組み合わせて使えなくなる。

#Julia言語

function f(x)
～
end

と一切型名を書かずに函数を定義していれば、多彩な型の引数で使える函数になるのに、

function f(x::Float64)::Float64
～
end

のように書いて、損しているのをよく見る。条件を摂動した結果の違いを見るためのコードも自分で書く羽目に陥いる。

#Julia言語

function f(x::Float64)::Float64
～
end

のようなコードを書いた途端に、NASAの人が

Solves the dependent input parameter problem
(no more "exploding complexity")

と強調しているところの、

no more "exploding complexity"

に至る道が閉ざされてしまうのです。

#Julia言語 上のVan der Pol方程式のPlutoノートブックを使う方法

⓪Juliaの公式バイナリをインストール

julialang.org/downloads/

v1.6.0-rc2またはnightly buildがお勧め。

julialang.org/downloads/nigh…

①juliaを起動して、]と押して

pkg> add Pluto PlutoUI Plots MonteCarloMeasurements

#Julia言語

②バックスペースを押して

julia> using Pluto; Pluto.run()

すると、ブラウザにPlutoノートブックのタブが開く。

③以下のリンク先のURLを貼り付ける。

gist.github.com/genkuroki/9013…

④OKボタンを押す。

#Julia言語

⑤すると、Plutoノートブックがダウンロードされて実行される！しばらく待つ。

⑥⑦⑧添付画像を参照。

Juliaを使ったことがない人であってもトラブらなければ10分もあればここまでたどり付けると思います！

#Julia言語 3.0 ± 0.1 のような「不定性を持った数」を実装しておき、Plots.jl用のレシピも実装することによって「魔法」が実現されています。

添付画像

① sin(3x)のプロット

② sin((3±0.1)x)のプロット

①と②の違いはパッケージの読み込みと a = 3.0 と a = 3.0 ± 0.1 の違いだけ！

https://twitter.com/sarkov28/status/1371355063440969728

#Julia言語

①質点の質量m=1、バネ定数k=3の調和振動子をDifferentialEquations.jlでYoshida6()法を使って数値的に解いている。

②MonteCarloMeasurements.jlを使ってバネ定数を k = 3.0 ± 0.03 とした場合。

この2つのコードの違いはほんの少ししかない❗️

gist.github.com/genkuroki/5230…

#Julia言語 このように、微分方程式をシンプレクティック法の1つであるYoshida6()法で数値的に解いて、パラメータを変化させたときの挙動を視覚化することを、空行を除けばたったの9行で実行できている！

こういう便利なものをNASAの宇宙開発で使っているわけです。

gist.github.com/genkuroki/5230…

#Julia言語

しかも、パラメータを変化させたときの挙動を視覚化を行うためには、本質的に、3.0に設定していたパラメータの値を 3.0 ± 0.03 に書き換えるだけで良い！

NASAの人もその点を強調しています。

この辺はこのスレッドのトップに貼り付けたYouTubeの動画を見ても認識し難い点だと思った。