Julia入門

国立研究開発法人日本原子力研究開発機構
理化学研究所
1

システム計算科学センター

革新知能統合研究センター (AIP)

永井佑紀

2

アウトライン
イントロダクション
数値計算をする、とはどういうことか
Juliaとはどんなプログラミング言語かJuliaでのライブコーディング：天体の運動
まとめ

その場でゼロからコーディング

3

イントロダクション

4

自己紹介
永井佑紀
専門：物性理論等

日本原子力研究開発機構システム計算科学センター
理化学研究所 革新知能統合センター(AIP)
スパコンを使って、数千-数万CPUコアMPI並列計算を実行している
最近の研究分野：機械学習を物理学に使う
量子多体系へとか格子量子色力学へとか、分子動力学へとか機械学習を「高速化」に使う話が好き（人間の代わりに頑張れ）
使えるプログラミング言語(得意な順)Julia, Fortran2008以降, Python, c++Juliaを知ったきっかけ
2016-2017(MIT客員研究員時代@米国ボストン)
同僚にJuliaを勧められたのでやってみる
ちょうどMITにいるしFortranと同じ配列が1始まりだし便利なパッケージがたくさんあるし速いし -> Juliaの数値計算で論文を書く5

今日の講義について目的
1. Juliaというプログラミング言語を知る
2. Juliaで数値計算するのが便利だということを知る
後で試したくなる？

1週間で学べる！Julia数値計算プログラミング

講談社

この本を読んでみたくなる？？

6

数値計算をする、とはどういうことか

7

数値計算とは？
数値計算 人間ができないor面倒or大変なことを計算機にやらせる
物理系の学科の学部で習うこと
力学、熱力学、統計力学、電磁気学、量子力学etc...
紙と鉛筆で学習：講義に数値計算は出てこない

何故？
容易に手で解ける問題しか扱っていないから手で解ける領域
残りは解けない8

物理で現れる計算
2つの物体同士に働く万有引力運動方程式を立てると⃗ma ⃗= FGMm
F ⃗ = − 3 r⃗r

2本の微分方程式となる
r:2つの物体間の距離dv⃗
GM= − 3 r⃗
dtr

dr ⃗= v⃗dt

これを解けば2つの物体の運動がわかる
2つなら手で色々やることで性質がわかる
N個の物体がある場合は？

9

物理で現れる計算
N個の物体同士に働く万有引力

Fi = − Gmimj∑ | r ⃗ − r ⃗ |2ji

jri: i番目の物体の位置

運動方程式を立てると⃗ma ⃗ = Fi

F i⃗ = − Gi

2N本の微分方程式となる
r i⃗ − r j⃗
d v i⃗=−Gmj∑ | r ⃗ − r ⃗ |2dt
ji

j

∑
jmimj

r i⃗ − r j⃗
| r i⃗ − r j⃗ |d r i⃗
= v i⃗
dt

これを解けばN個の物体の運動がわかる->手では解けない-> 自然界の振る舞いを調べるには、紙と鉛筆だけでは足りない

2

10

物理で現れる計算
1粒子の時間依存しないシュレーディンガー方程式

̂Hψn(x, y, z) = Enψn(x, y, z)ℏ∂∂∂Ĥ = −+++
V(x,y,z)2222m ( ∂x∂y∂z )222

微分方程式を満たすようなエネルギーEnと波動関数ψnを求める問題-> V(x,y,z)=0なら、手で（フーリエ変換で）解くことができる-> V(x,y,z)が特別な形をしている時のみ手で解ける
多くの場合は解けない
まだ電子1つの問題：現実には沢山電子がいる-> 自然界の振る舞いを調べるには、紙と鉛筆だけでは足りない

𝒞

11

物理で現れる計算

統計力学

2次元イジング模型

例：磁性J

σi: i番目のスピンの値,+1か-11
H( )⟨A⟩ =exp −A( ) 物理量の期待値∑( kBT )Z[]：あるスピン配置全ての可能なスピン配置に対して和を取る
電子はスピンを持つ

ある温度以下になると、突然磁石になる（相転移）->小さな磁石𝒞

向きが全部揃うと磁石になる𝒞𝒞-> 自然界の振る舞いを調べるには、紙と鉛筆だけでは足りない12

物理で現れる計算

統計力学-> 自然界の振る舞いを調べるには、紙と鉛筆だけでは足りない
例：磁性J

2次元イジング模型
ただし、数学を駆使すれば解けることもある
2次元イジング模型は手で解けている
手で解けるに越したことはない

電子はスピンを持つ->小さな磁石
向きが全部揃うと磁石になる

2x +2x +1 = 0の解を数値的に探すよりも2(x+1) =0, x=-1とわかっていた方が良い-> 手元の有益な情報が多いほど計算コストが下がる13

数値計算をする、とはどういうことか手で解ける領域

手で解けない面白い問題が沢山ある
計算機は道具である

残りは解けない
物理屋は数学と計算機を道具とする
数学科ほど数学を深めない
情報系学科ほど計算機を深めない
どちらも例外あり

現実世界に存在する面白いものを対象とする
どんな方法であれ解ければ良い
物理学者が数学者を刺激する場合もあり14(多くの)物理屋にとって数値計算とは物理の問題に注力したい
物理以外の問題に手間をかけ過ぎたくない
自動車を運転して遠くに行きたいのであって、そのために一から車の部品を組み立てたいわけじゃない

現実世界に存在する面白いものを対象とする
どんな方法であれ解ければ良い

一方、目的地に到達するためにはカスタマイズした移動手段が必要な場合もある
カスタマイズが楽しくなってくる場合もある
計算機科学者との違い：見据えているのは目的地

自由度が高い、書きやすい、そして速いプログラミング言語がいい
Julia言語??

15

数値計算の種類
1. 理論はすでにあり、適応する対象が無数にあるケース

Fortran,C++

例：固体物理における「第一原理計算」
2.どのような理論を使うか試行錯誤が必要なケース
いっぱいある

Fortran,C++,Python

グリーン関数による摂動論
ハミルトニアンの数値的対角化

3. 機械学習が絡むケース

Python

最先端の機械学習の手法は個人で実装するには大変すぎる -> Pythonライブラリの使用

4. 上記の様々な組み合わせ

Juliaは2,3,4で有用

16

数値計算の種類
1. 理論はすでにあり、適応する対象が無数にあるケース
密度汎関数理論で多体電子系を扱う
固体物理における「第一原理計算」
原子位置の情報がわかれば（原理的に）全てがわかる
物質をデザインしたり、実験結果と比較したりもちろん、手法を開発する研究もある

数値計算として重要な点:可能な限り高速に動くクラスターマシンやスパコンなどで大規模に動かせる-> FortranやC++で書かれていることが多い
この領域はすでに良いコードが複数あるので、そのまま使った方が良さげ
使い方さえわかれば、内部のコーディングがどうなっているか気にする必要がない

17

数値計算のためのプログラミング言語
数値計算をする場合に、どのプログラミング言語を選べば良いのか
原理的には、どんな言語を選んでも数値計算は可能

例えば
FORTRAN: (formula translation)
有名なFORTRAN

1957年にコンパイラがリリース

FORTRAN77 1978年
1991-92年
Fortran90

45年前！
31年前！

速い！ 数値計算でいまだに使われているプログラミング言語

Python

機械学習分野での基本言語の地位を確立（研究における英語みたいな）
Pythonは「知っていないとまずい」言語に他の言語で書かれた高速なコードを呼ぶ、という形で使われることが多い18

数値計算のためのプログラミング言語
数値計算をする場合に、どのプログラミング言語を選べば良いのか
原理的には、どんな言語を選んでも数値計算は可能

物理の問題に注力したい

どんなプログラミング言語が良いか

数式で表現されたものを数値計算したい
Mathematica等
その言語特有のお作法に悩まされたくない
Python: forループが遅いためにnumpyの書き方を習得する必要がある
c++: 行列やベクトルや複素数を扱うためには言語に慣れる必要がある
なるべく「楽に」数値計算がしたい19

数値計算の例
1粒子の時間依存しないシュレーディンガー方程式

̂Hψn(x, y, z) = Enψn(x, y, z)ℏ∂∂∂Ĥ = −+++
V(x,y,z)2222m ( ∂x∂y∂z )2221. 計算機では連続的な量を扱えないので、微分を差分に直す
2∂ ψ(x)
1h≡x−x∼(ψ(x)−2ψ(x)+ψ(x))
i+1ii+1ii−1∂x 22h 2-> シュレーディンガー方程式は行列の固有値問題になる
2. NxNyNz x NxNyNzの行列Hを作成
2と3を計算機にやらせる
Pythonは（素朴にやると）2が遅い
3. 固有値問題Hψ=Eψを解く
慣れれば難しくないが…… <-C++は3のためにEigen等外部ライブラリをインストールする必要がある20

Juliaはどんなプログラミング言語か21

数値計算のためのプログラミング言語
数値計算をする場合に、どのプログラミング言語を選べば良いのか
原理的には、どんな言語を選んでも数値計算は可能

物理の問題に注力したい

どんなプログラミング言語が良いか

物理の数値計算では線形代数を多用する
特殊関数も使う
積分も微分もする-> 数値計算のできる「電卓」が欲しい！

使うのが簡単がいい
インストールとかコンパイルとかに時間を溶かしたくない

バグとりが楽ならなお良い

Julia言語が最適22

Juliaはどんな言語か
1次元シュレーディンガー方程式

̂Hψn(x) = Enψn(x)2ℏ∂Ĥ = −+V(x,y,z)
2m ∂x 2差分化して解いて、V=0の解析解と比較してみる

本当にこれだけのコード
他に”おまじない”はいらない23

Juliaはどんな言語か
1次元シュレーディンガー方程式

̂Hψn(x) = Enψn(x)2ℏ∂Ĥ = −+V(x,y,z)
2m ∂x 2差分化して解いて、V=0の解析解と比較してみる

julia qm.jlと実行すると

本当にこれだけのコード
他に”おまじない”はいらないグラフができる

24

Juliaはどんな言語か
トポロジカル超伝導体のハミルトニアン
HBdG =HN(kx, ky, kz)

ΔΔ†−HN(−kx, − ky, − kz)*)(ϵ(kx, ky, kz) + M(kx, ky, kz)
HN =0

P3(kx, ky, kz)

P1(kx, ky, kz) + iP2(kx, ky, kz)

ギリシャ文字（日本語も）使える
行列が直感的に書ける

0

P3(kx, ky, kz)

P1(kx, ky, kz) − iP2(kx, ky, kz)ϵ(kx, ky, kz) − M(kx, ky, kz)

0ϵ(kx, ky, kz) + M(kx, ky, kz)

P1(kx, ky, kz) + iP2(kx, ky, kz)−P3(kx, ky, kz)

0

P1(kx, ky, kz) − iP2(kx, ky, kz)−P3(kx, ky, kz)ϵ(kx, ky, kz) − M(kx, ky, kz)

25Juliaはどんな言語か
例：“トポロジカル超伝導体の励起スペクトルの波数空間依存性を詳細に見たい”
適当な範囲で3次元k空間での固有値を求めて3次元スプライン補間を行なってより細かな波数空間のメッシュを切ってファイルに書き出す

行列の対角化は1行で書ける

26

Juliaはどんな言語か“量子力学を機械学習で解く最新論文の結果(の一部)を再現したい”̂ n(x) = Enψn(x)
Hψ

2

ℏ
∂
Ĥ = −+V(x,y,z)
2m ∂x 2En =̂∫ dxψn(x) Hψn(x)††∫ dxψn(x) ψn(x)エネルギー固有値Enが最小となるような波動関数を機械学習で見つけるψn(xi) = f(xi)f(x):ニューラルネットワーク

A simple method for multi-body wave function of groundand low-lying excited states using deep neural network
Tomoya Naito, Hisashi Naito, Koji Hashimoto, arXiv:2302.08965

やること

1. ハミルトニアンの微分演算子を差分化
2. エネルギーが最小となるような波動関数を探索
論文では、Pythonの機械学習フレームワークTensorFlowを使用しているらしい27

Juliaはどんな言語か“量子力学を機械学習で解く最新論文の結果(の一部)を再現したい”
ハミルトニアンの微分演算子を差分化

using Flux
n1 = 4model = Chain(x ->[x],Dense(1,n1,Flux.softplus),Dense(n1,1,Flux.softplus),x -> sum(x)) |> Flux.f64

ニューラルネットワークを構築
隠れ層1層、活性化関数softplusのニューラルネット

あとは訓練するだけ

調和振動子ポテンシャル項

運動エネルギー項

function energy(model,xs)
ψ = model.(xs)c = ψ'*ψE = ψ'*H*ψ/creturn E
end

En =̂ n(x)∫ dxψn(x) Hψ†∫ dxψn(x)†ψn(x)エネルギーの定義

PyTorchやTensorFlowと違い、テンソルにする必要がない

28

Juliaはどんな言語か“量子力学を機械学習で解く最新論文の結果(の一部)を再現したい”
Adamというオプティマイザーを使用
ハミルトニアンの微分演算子を差分化

opt = Flux.setup(Adam(), model)
for i=1:50000
energyを最小化するようにmodelを訓練
Flux.train!(energy, model,(xs,),
opt)
ニューラルネットワークを構築
e = energy(model,xs) modelを使ったenergyの計算
if i % 1000 == 0
調和振動子ポテンシャル項
println("i= $i energy = $e")function energy(model,xs)
end
̂ (x)
∫ dxψ (x) Hψψ = model.(xs)E =∫ dxψ (x) ψ (x)c = ψ'*ψ
end
これだけで最新の物理学with機械学習
E = ψ'*H*ψ/cn

訓練はこれだけend
運動エネルギー項
return E†nnn†

nの結果を再現できるエネルギーの定義

Juliaはどんな言語か29

初学者がとっつきやすいシンプルな文法 and FortranやC言語並に高速
重要な点：複雑なコードでも高速！

格子量子色力学（格子QCD）のJuliaコード、LatticeQCD.jl :富谷氏の個別講義Fortranコードに匹敵する速度
注：Fortranコードはガチチューニング済みコードではない30

Julia言語2018年にバージョン1がリリースされた比較的新しい言語
特徴：数式に近い形で計算ができる
無料：インストールも簡単

電卓のような手軽さ and 高速！

https://julialang.org/downloads/

自分のOSに合わせてダウンロードする31

Julia言語の使い方
大きく分けて三種類の方法がある
1. アプリを起動する（ダブルクリックで開く）

手軽：電卓のように使える
2. Jupyter notebookを使うブラウザで使える。Pythonのnotebookと似た操作感。Mathematicaに似ている

3. 端末で実行する

重たい計算をしたい時やクラスターやスパコンで使いたい時複数の計算を投げたい時

Julia言語の使い方
大きく分けて三種類の方法がある
1. アプリを起動する（ダブルクリックで開く）

REPL(Read Eval Print Loop)と呼ばれているファイルに書いたコードも呼べる

32

33

Julia言語の使い方
大きく分けて三種類の方法がある
2. Jupyter notebookを使うブラウザで使える。Pythonのnotebookと似た操作感。Mathematicaに似ている34

Julia言語の使い方
大きく分けて三種類の方法がある
2. Jupyter notebookを使うブラウザで使える。Pythonのnotebookと似た操作感。Mathematicaに似ている

入力をしてshift + enterで結果が下に表示される数式も扱えるので、一つの計算ノートが出来上がる35

Julia言語の使い方
大きく分けて三種類の方法がある
3. 端末で実行する

重たい計算をしたい時やクラスターやスパコンで使いたい時
julia findnodes.jl > output.txt

FortranやCのようなコンパイル作業がいらない
どの方法でやっても同じ。
“ある（誰かが作った）コードを走らせたいのに、コンパイルが通らなくて動かせない” ということがない36

Juliaの文法
詳しい事は電卓的に使える

虚数はim
2*imじゃなくて2imで良い三角関数や指数関数はそのまま使える数値計算に特化したJuliaの解説書
量子力学や統計力学等の問題を解く

円周率はπをそのまま使える37

Juliaの文法
詳しい事は
行列の定義

文字列は非可換なので+じゃなくて*ベクトルの定義

数値計算に特化したJuliaの解説書
量子力学や統計力学等の問題を解く

aの転置複素共役はaʼ行列とベクトルの積Juliaの文法
詳しい事は

38

物理屋にとって良い点：ギリシャ文字が使える
シンプルな関数の定義
functionとして関数を定義してもよい

線形代数も特殊関数もシンプルに使える

数値計算に特化したJuliaの解説書
量子力学や統計力学等の問題を解く

行列Aのエルミート共役はAʼで計算できる

39

Julia言語の特徴：速さ
Juliaは、”コンパイル”という作業をユーザーで行わない Pythonのようなスクリプト言語みたい書いてすぐ実行という手軽さが実現
Juliaは型を指定しない

しかし、FortranやC言語に匹敵する速さ
実行時コンパイル(JITコンパイル)という技術を採用しているため初めて関数が呼ばれた時、動的にコンパイルが実行されている型を指定する煩雑さから解放されている

1. 型の推論
2. 高度な最適化
3. コンパイル

1.整数が入り整数が返ることが特定されている
2. a=1+2はこの時点で3が代入されているし、aという変数は代入で消去されている

40

Julia言語の特徴：速さ
初めて関数が呼ばれた時、動的にコンパイルが実行されているtest(2)1. 型の推論
2. 高度な最適化
3. コンパイル

test(2.1)

引数が2.1と整数ではないので、倍精度実数用のコードが作られる。倍精度用の和と積が用意されている
引数にどんな型が入っているかで、その都度最適なコードがコンパイルして生成される

41

Juliaの特徴：パッケージ管理システム
他の人が作ってくれているライブラリを簡単に導入することができる
FotranやC:自分の環境でコンパイルできるかどうかがわからない
それを動かすために別のライブラリが必要かもしれない
そもそもライブラリが置かれている共通スペースがない

WindowsでもMacでもLinuxでも同じようにインストールできる
REPL上で]キーを押す

add Fluxなどとする

依存パッケージも含めて自動的にインストールされる

数値計算を行う上で、他の人が開発したパッケージを使えるというのは極めて重要「車輪の再発明」を避ける
他の人が動作検証している（はず）なので安心して使える42

Juliaの特徴：パッケージ管理システム
例：SparseArrays.jl
疎行列（要素がほとんどゼロの行列）を簡単に扱える

SpecialFunctions.jl

TightBinding.jl

様々な特殊関数が扱える自作。公式ライブラリ化

QuadGK.jl
1次元数値積分ができる

ベッセル関数の直交性の確認43

Juliaの特徴：パッケージ管理システム“Juliaって新しい言語だからパッケージがPythonより少ないんでしょ…”
大丈夫！Pythonのコードも簡単に呼べるよPythonのライブラリ、TensorFlowのKerasを使ってニューラルネットワークを定義してみる
using PyCall
using Plots
const tensorflow = pyimport("tensorflow")
const keras =tensorflow.keras
const layers =keras.layers
const optimizers =keras.optimizers
function build_model(d_input,d_middle)
inputs = keras.Input(shape=(d_input,))x = layers.Dense(d_middle, activation="relu")(inputs)
y = layers.Dense(1)(x)adam = optimizers.Adam()model = keras.Model(inputs=inputs, outputs=y)model.compile(optimizer=adam,
loss="mean_squared_error")
return model
end

Juliaで作った配列をそのまま扱える
numpy配列としてインプットされる

c言語も呼べるしFortranも呼べる44

Julia言語の特徴:多重ディスパッチ
オブジェクト指向言語ではないが、ほとんど等価なことが可能
多重ディスパッチという機能で実装されている

FORTRAN,C

オブジェクト指向

犬が鳴く(10回)犬.鳴く(10回)鳴く(犬,10回)

蛙が鳴く(3回)

蛙.鳴く(3回)鳴く(蛙,3回)異なる関数を呼ぶ必要がある

多重ディスパッチ

C++, Python等

一つ目の引数を特別扱い

一つ目と二つ目に何が入るかで挙動が変わる

ユーザーは、犬か蛙かを気にせずに、「鳴く」を選べば良い45

Julia言語の特徴：多重ディスパッチ
オブジェクト指向“動物”クラスに「鳴く」を定義
動物クラスを継承した

多重ディスパッチクラスの代わりに「型」を使う“動物”という抽象型を定義“動物”を引数とする「鳴く」を定義
鳴く(a::動物,回数)“犬”クラスで「鳴く」をオーバーライド:動物型の下部タイプでワンワン！

ある「犬」型を定義

鳴く(a::犬,回数)
ワンワン！ を定義

動物クラスを継承した“蛙”クラスで「鳴く」

動物型の下部タイプでをオーバーライド:

ある「蛙」型を定義

ケロケロ！

を定義

犬と蛙が同時に鳴く場合のみゲゲゲワンと鳴く
シチュエーションの場合は？

鳴く(a::蛙,回数)
ケロケロ！ を定義

鳴く(a::犬,b::蛙,回数)
を定義

46

Juliaでのライブコーディング

47

万有引力に従う天体の運動
2つの物体同士に働く万有引力⃗⃗
ma = F

無次元化しておくGMm⃗
F = − 3 r⃗r

時間を離散化して、各時刻での天体の位置をシミュレートする

2本の微分方程式となる

dv⃗
GM= − 3 r⃗
dtr

dr ⃗= v⃗dt

簡単のため、z方向には動かないとする
原点に天体がいて、その周りを回っている

これを解けば2つの物体の運動がわかる
時間に関する微分方程式を逐次的に解くことにする

これをJuliaでコーディングしてみるエネルギーが保存するような手法を使ってみる-> リープフロッグ法を採用

48

万有引力に従う天体の運動時間に関する微分方程式を逐次的に解くことにする

1 2
12
E = (u + v ) −2r
リープフロッグ法：

エネルギーが保存するような手法を使ってみる-> リープフロッグ法を採用
速度と位置を交互に更新する

時刻tでの運動量は最初に必要なt=-Δt/2の時のuはu(−Δt/2) = u(0) − Φ(x)Δt/2

Juliaでのライブコーディング

49

万有引力に従う天体の運動
nprint:何回に1回書き出すか

u(t): uu(t+Δt/2): u̲pu(t-Δt/2): u̲mu2 + v2
E=−2v(t): v
v(t+Δt/2): v̲p

アニメーションのためのパーツv(t-Δt/2): v̲musing Plots
ENV["GKSwstype"] = "nul"anim = Animation()Δt = tmax/nn:離散点の数
tmax:最終時刻

1x2 + y2u(−Δt/2) = u(0) − Φ(x)Δt/2

plt = scatter([0,x],[0,y],label="t =$t",xlims=(-2,2),ylims=(-2,2),aspect_ratio = 1)
frame(anim,plt)

gif(anim, filename, fps = 30)

E0: 最初のエネルギー
E: 最初のエネルギー

これらの材料でライブコーディング

x = 1.0
y = 0.0
u = 0.0
初期パラメータv = 1.0n = 10000
tmax = 20
Δt = tmax/nt = 0.0
filename = "anim_v$v.gif"
nprint = 50

50

3体問題：太陽、地球、月それぞれが万有引力に従う

Fi = − Gmimj-> r1,r2,r3,v1,v2,v3みたいにすると大変∑ | r ⃗ − r ⃗ |2ji

Juliaっぽく書くには？

j例えばstruct Body
r::Vector{Float64}
v::Vector{Float64}
f::Vector{Float64}
mass::Float64
function Body(r,v,mass)f = zero(r)return new(r,v,f,mass)end
end

mutable struct SolarSystembodies::Vector{Body}
nbody::Int64
function SolarSystem()bodies = Body[]
nbody = 0
return new(bodies,nbody)end
end

天体の情報を格納する型

太陽系全体の情報を格納する型

これらを使って微分方程式を解く

51

3体問題：太陽、地球、月それぞれが万有引力に従う

Fi = − Gmimj-> r1,r2,r3,v1,v2,v3みたいにすると大変∑ | r ⃗ − r ⃗ |2jij

標準のpush!をSolarSystemに対応させる
function Base.push!(s::SolarSystem,body::Body)
push!(s.bodies,body)s.nbody = length(s.bodies)endfunction Base.getindex(s::SolarSystem,i)return s.bodies[i]ends[1]やs[2]で天体の情報を取り出せるようにする
多重ディスパッチの利用の例

Juliaっぽく書くには？
function calc_force!(s::SolarSystem)
nbody = s.nbody
for i=1:nbodys[i].f .= 0endfor i=1:nbody
r_i = s[i].r
mass_i = s[i].mass
for j=i+1:nbody
r_j = s[j].r
mass_j = s[j].mass
rij = r_i - r_js[i].f .+= -G*mass_i*mass_j*rij/norm(rij)^3s[j].f .+= -s[i].fend
end
end

万有引力の計算

52

3体問題：太陽、地球、月それぞれが万有引力に従う

Fi = − Gmimj-> r1,r2,r3,v1,v2,v3みたいにすると大変∑ | r ⃗ − r ⃗ |2jij

function update!(s::SolarSystem,Δt)
nbody = s.nbody
for i=1:nbody
mass = s[i].mass
s[i].v .+= s[i].f*Δt/(2*mass)s[i].r .+= s[i].v*Δtendcalc_force!(s)for i=1:nbody
mass = s[i].mass
s[i].v .+= s[i].f*Δt/(2*mass)endend

速度ベレル法で時間発展Juliaっぽく書くには？
function calc_force!(s::SolarSystem)
nbody = s.nbody
for i=1:nbodys[i].f .= 0
end
for i=1:nbody
r_i = s[i].r
mass_i = s[i].mass
for j=i+1:nbody
r_j = s[j].r
mass_j = s[j].mass
rij = r_i - r_js[i].f .+= -G*mass_i*mass_j*rij/norm(rij)^3s[j].f .+= -s[i].fend
end
end

万有引力の計算

53

3体問題：太陽、地球、月それぞれが万有引力に従う

Fi = − Gmimj-> r1,r2,r3,v1,v2,v3みたいにすると大変∑ | r ⃗ − r ⃗ |2jij

solarsystem = SolarSystem()sun = Body([0,0],[0,0],mass_Sun )push!(solarsystem,sun)r_Earth = [radius_Earth,0]v_Earth = [0,velocity_Earth]earth = Body(r_Earth,v_Earth ,mass_Earth )push!(solarsystem,earth)r_Moon = [radius_Earth+radius_Moon ,0]
v_Moon = [0,velocity_Earth+velocity_Moon]
moon = Body(r_Moon,v_Moon ,mass_Moon )push!(solarsystem,moon)

初期速度や位置を入れる
push!で太陽系に天体を追加しているJuliaっぽく書くには？
for i = 2:n
update!(solarsystem,Δt)t += Δtend
本体コード（可視化部分は省略）
この方法だと、pushで別の天体をいくらでも追加できる天体を追加しても本体コードの変更はない

54

3体問題：太陽、地球、月それぞれが万有引力に従う

Fi = − Gmimj-> r1,r2,r3,v1,v2,v3みたいにすると大変∑ | r ⃗ − r ⃗ |2
ji

Juliaっぽく書くには？

j

質量とか速度をちゃんと入れると、365日で一周するのを確認できる

55

まとめ

56

Julia言語とは
数値計算に最適な言語
Pythonのようなシンプルさ
C,Fortranに匹敵する速さ
matlabのような線形代数操作
Mathematicaのように積分したりできる
新しい言語なので、数値計算における”痒いところに手が届く”ような言語設計になっている

この春の学校中、Juliaに関する質問はいつでも受け付けます

57

58

Julia言語とは数値計算に最適な言語
Pythonのようなシンプルさ
C,Fortranに匹敵する速さ

数値計算しようぜmatlabのような線形代数操作
Mathematicaのように積分したりできる
新しい言語なので、数値計算における”痒いところに手が届く”ような言語設計になっている

この春の学校中、Juliaに関する質問はいつでも受け付けます