Julia 1.0が出て、しばらくは言語仕様が大幅に変わることがなくなった。
ということで、物理分野で数値計算をする際に役に立ちそうなコードをサンプルとしてまとめておくことにする。Juliaは色々な書き方ができるので、ここにあるのはあくまで参考としてほしい。なるべく少ない労力で他の言語から移行できるような書き方を目標としてみた。
#対象とする人
- FortranやCを使って、簡単な数値計算をする人
- PythonでPythonのお作法を使わずにCやFortranと同じノリで書いたら実用的ではないなと感じた人
Juliaとは
- 書きやすい
- 明示的にコンパイルをする必要がない。
- 配列のインデックスが1始まりなのでFortranからの移行にやさしい
- 速い(Pythonでnumpyで頑張った場合とほぼ同等の速度が何も考えなくても出る)
- 行列の対角化、関数の積分、疎行列の取り扱い、などが便利
- CやFortranと同等の速度がでるPythonのような言語
- オブジェクト指向ではない（多重ディスパッチという概念を用いる）。しかし、物理で出てくるような簡単な数値計算をする際にはオブジェクト指向がなくても問題はなく、多重ディスパッチが便利。
- 日本語の文献が少ない（最近増えてきた）
参考：
「なぜ僕らはJuliaを作ったか」
https://www.geidai.ac.jp/~marui/julialang/why_we_created_julia/index.html
参考となりそうなサイト
Julia 1.0.0以降に対応してそうなものを以下に紹介しておく。
「Julia入門 - 入門者がまず読むべき参考サイトや本のまとめ【随時更新】」
http://st-hakky.hatenablog.com/entry/2018/02/01/100000
「Julia 入門①」
https://qiita.com/torininaritahi/items/6a3a431e1341e9057877
「Mac で Julia #1」
https://qiita.com/qoAop/items/c4df6794e44da21ee7e4
インストール
どのOSでも、
https://julialang.org
からバイナリファイルをダウンロードすれば、root権限なしですぐにJuliaを使用することができる。
使い方
基本的な計算
まず、基本的な計算がJuliaでどう書けるかをみてみる。
整数、実数、複素数
Juliaは型を自動的に推論してくれるので、整数、実数、複素数を定義して、普通に四則演算ができる。
コメントアウトは#を使う。全体をコメントアウトしたい場合には、#= と =#を使う。虚数はimである。
#整数
a = 1
b = 2
c = a + b
println(c)
d = a*b
println(d)
#=
以下は
実数の場合
=#
a = 1.0
b = 2.2
c = a + b
println(c)
d = a*3b
println(d)
#整数の複素数
a = 1+2im
b = 3+4im
c = a + b
println(c)
d = a*3b
println(d)
#実数の複素数
a = 1.0+2.0im
b = 3.0+4.0im
c = a + b
println(c)
d = a*3b
println(d)
#= 円周率πも普通に使えるし、ギリシャ文字も使用可能。
cosやsinもそのまま使用可能
=#
α = π
β = 3π
γ = cos(π+β)
println(γ)
#= 日本語だって使える
=#
りんご = 3
ばなな = 4
総数 = りんご +ばなな
println("りんごの数 $りんご ばななの数 $ばなな 総数は$総数")
あたらしいりんご = 3
総数 += あたらしいりんご #+=で追加できる。
println("総数は",総数,"となった")
ここに書いたように、数字と虚数2.0imや、数字と変数3bなど、掛け算の記号なしでくっつけて書くことができる。また、基本的な数学関数はPythonのようにnp.とかをつけなくてもよく、cosとかsinでよい。標準出力はprintlnとprintがあり、printlnは最後で改行し、printは改行しない。文章のなかで$aのようにすると、変数aの値を出力できる。最後のように、カンマで区切ってもよい。
なお、割り算は、
#結果は実数になる
a = 12
b = 4
c = a/b
println(c)
#結果は整数になる
a = 12
b = 4
c = div(a,b)
println(c)
#割り切れない場合
a = 12
b = 5
c = div(a,b)
println(c)
#あまりの計算
a = 12
b = 5
d = a %b
println(c)
となる。ここで、/を使うと実数になり、divを使うと、整数になる。割り切れない場合には、最大の商となる整数がでる。あまりは、%を使う。
累乗は
a = 3
b = a^2
println(b)
と^の記号を使う。
関数と多重ディスパッチ
次に、FortranやCから来ると便利に思える機能として、関数（function）がある。
関数は、
#足し算
function tasu(a,b)
c = a+b
return c
end
#足し算と掛け算
function tasukakeru(a,b)
c = a+b
d = a*3b
return c,d
end
a = 1.0
b = 2.2
c = tasu(a,b)
println(c)
c,d = tasukakeru(a,b)
println("$c,$d")
a = 1.0 +3im
b = 2.2
c = tasu(a,b)
println(c)
c,d = tasukakeru(a,b)
println("$c,$d")
#functionは一行で書くこともできる。
tasu(a,b) = a+b
a = 1.0
b = 2.2
c = tasu(a,b)
println(c)
となる。functionと書いてendで閉じる方法と、一行で書く方法がある。
面白い点としては、functionの定義において、変数の型（実数だとか複素数だとか）を指定していないことである。つまり、変数に何がきても、そのまま計算ができる。これは、実数と複素数をわける必要がないことを意味しており、物理で使う数値計算では非常に役に立つ。
一方、変数の引数が異なると異なる挙動をさせたい場合もあるだろう。たとえば、複素数のときだけ2πを足しておきたいとか。そういうときは、
#足し算
function tasu(a,b)
c = a+b
return c
end
#足し算
function tasu(a::Complex,b)
c = a+b+2π
return c
end
a = 1.0
b = 2.9
c = tasu(a,b)
println(c)
a = 1.0+0.0im
b = 2.9
c = tasu(a,b)
println(c)
などと型を明示する。ここで、functionの名前が同じなことを着目してほしい。実際にこの関数を使うときに、変数の型を意識せずに、自動的に該当する関数tasuを呼んでくれる。このような機能を「多重ディスパッチ」と呼ぶ。実は、四則演算の+、-、*や/も同じようになっており、型に応じて必要な関数を呼んでいるのである。
行列とその演算
物理系の数値計算で避けて通れないのは、行列である。Juliaでは、配列の添字は1から始まる。そして、配列は一番左側の添字が回ってからそのさらに左側の添字が回るような形でメモリに格納されているため、左側の添字が回るようにアクセスするとメモリアクセス的に効率的である。なお、このメモリ配置方法はFortranと同じ順番であり、cとは逆になっている。
#行列の定義
A = [1 2
3 4]
println(A)
println(typeof(A))
#別の定義
A = [1 2;3 4]
println(A)
println(typeof(A))
#要素がゼロの2x2行列
A = zeros(Int64,2,2)
A[1,1] = 1
A[1,2] = 2
A[2,1] = 3
A[2,2] = 4
println(A)
println(typeof(A))
#要素が1の2x2行列
A = ones(Int64,2,2)
println(A)
println(typeof(A))
#ベクトルの定義
A = [1,2] #これはベクトル。つまり、1次元配列
println(A)
println("ベクトル ",typeof(A))
A = [1 2] #これは上記の
#=
A = [1 2
3 4]
の3,4を削ったかたちなので、1x2行列 =#
println(A)
println("2次元配列 ",typeof(A))
A = [1;2] #これは2x1行列
println(A)
println(typeof(A))
ここで、typeofというのは、その変数の型を返す関数で、Array{Int64,1}だと整数の1次元配列、Array{Int64,2}だと整数の2次元配列である。
行列の積は、
A = [1 2
3 4]
B = [10 20
30 40]
C = A*B
println(C)
C = B*A
println(C)
と*を使って書くことができる。これは、先ほどの多重ディスパッチで*が行列の積として定義されているからである。配列のそれぞれの要素に演算をしたい場合には、任意の関数にドット.をつければよくて、たとえば、ある行列Aにたいしてそれぞれの要素にsinを適用したい場合には、
A = [1 2
3 4]
B = sin.(A)
println(B)
とすればよい。ドット.は自前の関数にも使えるので、
A = [1 2
3 4]
B = [10 20
30 40]
function tasu(a,b)
c = a+b
return c
end
C = tasu.(A,B)
println(C)
ともできる。
対角化
行列の対角化がしたいときには、
using LinearAlgebra
H = [1 3
3 1]
e,u = eigen(H)
println(e)
println(u)
とすれば、固有値eとユニタリー行列Uが出てくる。これがちゃんと行列を対角化しているか確認するには、
D = u'*H*u
println(D)
とすればよい。ここで、行列の転置は'でできる。なお、複素数の行列だと、エルミート共役をとってくれる。
配列の準備とコピー
配列のコピーは
#コピーする
B = copy(A)
#別の方法。配列次元がわかっているならこちらのがみやすいかもしれない。
B = A[:,:]
でできる。なお、
B = A
としてしまうと、BとAが同一のメモリをさし示すことになり、Aを変えるとBも変わってしまう。
例えば、
println("A: $A")
B = A
println("A: $A")
println("B: $B")
A[1,1] = 100
println("A: $A")
println("B: $B")
とすると、AもBも変更される。ここは注意。
行列の一部を取得、行列ベクトル積
二次元配列として行列を持っているとき、その一部を取り出したいことがある。たとえば、上で対角化した行列の固有ベクトルを取り出したいとき。
上で得られた固有値固有ベクトルがちゃんと固有値固有ベクトルになっているかを確認してみると
u1 = u[:,1]
u2 = u[:,2]
r1 = H*u1 - e[1].*u1
println(r1)
r2 = H*u2 - e[2].*u2
println(r2)
となる。ここで、ベクトルu1に、スカラーe[1]をドット.を使って要素積をしている。
特異値分解
特異値分解は
using LinearAlgebra
B = rand(Float64,3,3) #randという関数は乱数を発生させてくれる。
println(B)
U,S,V = svd(B) #特異値分解 B = U*diagonal(S)*V'
C = U*Diagonal(S)*V'
println(C)
でできる。ここで、一次元配列SはDiagonalという関数で行列にすることができる、ということを使った。また、rand(Float64,3,3)というのは、倍精度実数(Float64)の3,3行列を作る、という関数である。
逆行列
逆行列は
using LinearAlgebra
D = inv(B)
F = D*C
println(F)
で計算できる。なお、LinearAlgebraは線形代数関連のパッケージであり、一度usingしておけば何度も書く必要はない。
線型方程式を解く
$Ax=b$のような線型方程式を解きたい場合には、
using LinearAlgebra
using Random
Random.seed!(13526463) #乱数のシードを固定する。
A = rand(Float64,3,3)
b = rand(Float64,3)
x = A \ b #Ax=bを解く。
println(A*x-b)
とすればよい。ここで、using Randomを使って、乱数のシードをRandom.seed!(13526463)で固定した。これをすることにより、計算のたびに同じ乱数が使われるため、デバッグが楽になる。
行列の初期化
行列(あるいは多次元配列)に具体的な値をいれる前に初期化をしたいこともある。
n = 3
m = 4
A = Array{Int64}(undef, n, m) #要素が整数のn x m の二次元配列
A[3,4] = 3
println(A[3,4])
A = Array{Float64}(undef, n, m) #要素が倍精度実数のn x m の二次元配列
A[3,4] = 3 #ここで整数を入れても、定義で要素はFloat64なので、
println(A[3,4]) #実数で出力される。
A = Array{ComplexF64}(undef, n, m) #要素が倍精度複素数のn x m の二次元配列
A[3,4] = 3 #ここで整数を入れても、定義で要素は複素数なので、
println(A[3,4]) #複素数で出力される。
A = Array{Array{Int64,1}}(undef, n, m) #要素が1次元整数配列のn x m の二次元配列
A[3,4] = [2,4]
println(A[3,4]) #1次元配列が出力される。
A = Array{Array{Float64,2}}(undef, n, m) #要素が2次元字数配列のn x m の二次元配列
A[3,4] = [2 3
2 9]
println(A[3,4]) #2次元配列が出力される。
配列の要素の型としては配列を入れることもできる。
ループ
forループのやり方
FortranではDoループ、Python、CではForループと呼ばれるループのJulia版について。
Juliaでは、
n = 10
#一つ目の書き方。Fotranっぽく。
for i=1:n
println("i=$i")
end
#二つ目の書き方。Pythonっぽく。
for i in 1:n
println("i=$i")
end
と書ける。
##If文
JuliaでのIf文は、
n = 10
for i=1:n
println("i = $i")
if i % 2 == 0
println("$i, even")
elseif i == 5
println("$i")
else
println("(^_^)v")
end
end
このような感じで書く。
また、一文で書く方法としては、
function test()
n = 10
j= 0
for i=1:n
j += ifelse(i % 2 == 0,1,0)
println("j = $j")
end
end
test()
がある。このifelseは、条件がみたされたときに、1を足し、そうじゃないときには0を足す。trueとfalseを両方評価して代入する仕様になっているため、1とか0などの場合は普通のIf文とくらべてかなり速い。ただし、両方評価してから代入するために、
a = -2
b = ifelse(a > 0,sqrt(a),0)
はエラーでとまってしまう。なぜなら、sqrt(a)をa= -2のときに評価してしまうからである。このときは、if文で普通にやったほうがよい。
関数化
Juliaが1.0.0になってから気をつけるべき点がある。
forループなどをコードにべたがきするのではなく、ifelse文の例のように、functionでラップしてから、それを呼ぶ形にしたほうがよい。functionごとにJuliaがコードを最適化するので、より高速になる。べたがきでのfor文は全く最適化されないため、遅い。つねにfunctionを使う癖をつけておくと、問題が起きにくい。この記事かこの記事以降に、書きやすくわかりやすい例を載せる予定である。
数値積分
物理分野で非常に使われる数値計算といえば、数値積分である。
物理でよく出てくる波数空間の積分：
I = \frac{1}{2π} \int_{-\pi}^{\pi} dk f(k)
をやってみよう。
関数$f(k)$は試しに
f(k) = \sin (k) + k^2
としてみる。積分値は
I = \frac{1}{2\pi}\frac{2}{3}(\pi)^3 = 3.2898681336964524
となる。
台形公式による積分
まず、台形公式による積分をやってみよう。
積分を$N$分割すると、
\frac{1}{2π} \int_{-\pi}^{\pi} dk \sim \frac{1}{2π} \sum_i f(k_i) \frac{2\pi}{N} = \frac{1}{N} \sum_i f(k_i)
となる。ここで、台形公式の最初の点と最後の点が同一の点であることを利用し、ただの和にすることができることを用いた。
Juliaではfunctionを引数にすることができる。したがって、
function daikei(f,N)
dk = 2π/(N-1)
fsum = 0
for i=1:N
k = (i-1)*dk - π
fsum += f(k)
end
fsum /= N
return fsum
end
f(x) = sin(x) + x^2
N = 400
fsum = daikei(f,N)
exact = ((π)^3/3 -(-π)^3/3)/(2π)
println("daikei $fsum, exact $exact")
で台形公式による数値積分ができる。
パッケージを使った数値積分
次に、Juliaのパッケージを使ってみよう。
１次元の数値積分パッケージはQuadGK.jl
https://github.com/JuliaMath/QuadGK.jl
がある。使うには、
Juliaを立ち上げて、]キーを押してPkgモードにして、add QuadGKをすればよい。
quadgk(f,a,b)は、関数f(x)とaからbの積分区間で積分ができる。
using QuadGK
f(x) = sin(x) + x^2
fsum2 = quadgk(f,-π,π)./(2π)
println("quadgk $fsum2, exact $exact")
やってみると、すごく正確に出る。
[2019年5月11日、記述が間違っていたので修正]
ここで、./(2π)のドットは、quadgkが(I,E)のような形で
積分値をI、エラーをEでアウトプットするためにつけている。ここで、[]ではなく()はタプルというものである。ドットをつけると要素ごとに演算される。
データ点を用いた補間と積分
数値計算においては、ある関数$f(x)$の値が$x = x_1,x_2,\cdots,x_N$と有限個の点でのみ求められていることがある。このような場合に、持っている点以外での値が欲しいときには、補間が有効である。
補間の一つとしてスプライン補間というものがある。
これを扱うパッケージは、
https://github.com/kbarbary/Dierckx.jl
で、]を押してパッケージモードにしてから
(v1.0) pkg> add Dierckx
でインストールできる。
ある値$x$とその関数$y(x)$の組が与えられたとき、
using Dierckx
x = [0., 1., 2., 3., 4.]
y = [-1., 0., 7., 26., 63.] # x.^3 - 1.
spl = Spline1D(x, y)
とすれば、
spl([1.5, 2.5]) # result = [2.375, 14.625]
spl(1.5) # result = 2.375
とすることで任意の点での値を得ることができる。この例では$y(x) = x^3-1$である。
また、微分、積分、ゼロ点探索はそれぞれ
derivative(spl, 1.5) # x=1.5での微分; 値は 5.75
integrate(spl, 0., 4.) # x=0 から x=4までの積分; 値は 60.0
roots(spl) # 0点を求める。値は [1.0]
でできる。2次元も同様に、
x = [0.5, 2., 3., 4., 5.5, 8.]
y = [0.5, 2., 3., 4.]
z = [1. 2. 1. 2.; # サイズは (length(x), length(y))
1. 2. 1. 2.;
1. 2. 3. 2.;
1. 2. 2. 2.;
1. 2. 1. 2.;
1. 2. 3. 1.]
spline = Spline2D(x, y, z)
でできる。
特殊関数
Juliaでは特殊関数も簡単に扱うことができる。
Juliaを立ち上げて、]をおしてPkgモードにして、add SpecialFunctionsとすれば、特殊関数を扱うことができる。たとえば、Bessel関数であれば、
using SpecialFunctions
z = 0.2
a = besselj0(z) #ν=0のベッセル関数
println(a)
z = 0.1
for ν=0:10 #J_{ν}(z=0.1)の値
Jν = besselj(ν,z)
println("J_{$ν}(0.1) = $Jν")
end
として使うことができる。
Bessel関数のゼロ点が欲しいのであれば、add FunctionZerosとして、
using FunctionZeros
n=10
for ν=0:10 #J_{ν}(x)のn=10番目のゼロ点の位置
Jν_zero = besselj_zero(ν, n)
println("J_{$ν}の$n 番目のゼロ点の位置は $Jν_zero")
end
として得ることができる。
データの出力
データ出力の基本
計算したデータを出力する方法について。色々なやり方があるが、シンプルな方法を記す。
name1 = "test"
name2 = "_file.dat"
filename = name1*name2
#新しく書き込むファイルの用意
fp = open(filename,"w")
for i=1:10
println(fp,i)
end
#ファイルを閉じる
close(fp)
ここで、name1とname2は文字列型で、文字列型同士を結合するには、*を使う。+ではないのは、文字列型同士が可換ではないから、らしい。fpを開いて、閉じる。これで書き込みができる。Pythonっぽく書きたいのであれば、
open("r"*filename, "w" ) do fp
println(fp,atan(1)*4)
end
でもよい。
直接プロット
Juliaには、得られた結果をファイルに出力するのではなくそのままグラフとしてプロットする方法も色々ある。Plotsパッケージなどを使えば良い。
http://docs.juliaplots.org/latest/
可視化についてはさまざまなパッケージがあるので割愛する。FortranやCからやってきた場合にはgnuplotなどを使って可視化することが多いと思われるが、Julia自体にプロットさせるというのも面白いので考慮に値すると思う。
データの入力
数値データを読み込む方法についても記す。これもいろいろなやり方があるが、FortranやCからやってきた方にとってわかりやすそうな書き方の例をあげる。
例えば、読み込むファイルが
0 1.00 3.4 test0
1 -1.5 6.499 test1
2 3.01 4.5 test2
というものだとする。ここではスペースで区切っているが、tabでも構わない。
これを読み込むには、
function test(inputdata)
numofdata = countlines(inputdata) #ファイルの中のデータの数をカウント。ここでは3
data = readlines(inputdata) #データファイルの読み込み。配列として確保される。
data_int = [] #空の配列
data_float = [] #空の配列
data_string = [] #空の配列
for inum = 1:numofdata
u = split(data[inum]) #i番目のデータ。
println(u) #この状態ではstringになっている。
di = parse(Int64,u[1]) #整数の読み込み。parseは指定した型に変換する。
df = parse.(Float64,u[2:3]) #実数二つ読み込み。parse.と.をつけたことでそれぞれにparseを適用。
ds = u[end] #最後の文字列の読み込み。endは配列の最後を意味する。
push!(data_int,di) #データの追加
push!(data_float,df) #データの追加
push!(data_string,ds) #データの追加
end
println(data_int)
println(data_float)
println(data_string)
return data_int,data_float,data_string
end
inputdata = "test.dat"
test(inputdata)
などとする。openを使ったやり方もあるが、ここではreadlinesで一気に読み込んだ。
#マクロ
Juliaにはマクロというものがある。詳しいことは他の文献を参照ということにする。ここでは、@timeというマクロを紹介する。これは、時間を計測するマクロである。
@time fsum = daikei(f,N)
exact = ((π)^3/3 -(-π)^3/3)/(2π)
println("daikei $fsum, exact $exact")
これで、functionにかかった時間を測定できる。
#並列化
簡単な並列化についても記しておく。この部分は、Julia 0.6での記事
「Juliaでforループをプロセス並列化してみる」
https://qiita.com/cometscome_phys/items/44d09a6b67bd48a1fbfb
の一部をJulia 1.0でも動くようにしたものである。
並列化を行いたいときには、julia -p n test.jlのような形で実行する。ここでnは並列数である。
コードにおいては、
using Distributed
@everywhere using LinearAlgebra
とDistributedパッケージを読み込む。また、並列化したときの他のプロセスでもパッケージを使いたいので、usingするパッケージには@everywhereをつける。
たとえば、
function svdi(nm,i)
#Matrix{Float64}(I, nm, nm)はサイズnmの単位行列を作る。
u,s,v = svd((3*i).*Matrix{Float64}(I, nm, nm)+cos(i)*ones(Float64,nm,nm))
return s
end
nmax = 1000
nm = 100
A = [] #空の配列
@time for i=1:nmax
s = svdi(nm,i)
push!(A,s) #Aという配列にsを追加する。
end
println(A[6])
というようなforループを並列化する際には、
@everywhere function svdi(nm,i)
#Matrix{Float64}(I, nm, nm)はサイズnmの単位行列
u,s,v = svd((3*i).*Matrix{Float64}(I, nm, nm)+cos(i)*ones(Float64,nm,nm))
return s
end
@time A = pmap(i -> svdi(nm,i), 1:nmax)
println(A[6])
とすればよい。ここで、使っているfunctionがどのプロセスにも見えるように、@everywhereをつけておく。forループの並列化にはpmapを使う。ただし、pmapを使うときには、ループの中の処理が重くなければ、並列化のオーバーヘッドで並列化する前より遅くなってしまうことがある。このコードの場合には、特異値分解する行列が十分に大きければよい。
詳細は
https://docs.julialang.org/en/v1/manual/parallel-computing/index.html
を参照すること。
バッチジョブなどでコントロールされたクラスターマシン上では、
julia --machinefile=$PBS_NODEFILE /path/to/your/home/dir/test_julia.jl
で動くはずである。
[2018/10/25追記]：Julia 1.0では
julia --machine-file=$PBS_NODEFILE /path/to/your/home/dir/test_julia.jl
としなければならないようだ。
#微分方程式
##常微分方程式を解く
常微分方程式：
\frac{du}{dt} = f(u)
を解くには、DifferentialEquationsというパッケージを使用すればよい。
http://docs.juliadiffeq.org/latest/
以下はこのリンク先のドキュメントを参考にしている。
このパッケージを使うには、Juliaを起動し、]を押してパッケージモードにして、add DifferentialEquationsとすればよい。
例えば、
\frac{du}{dt} = 1.01 u
という微分方程式を初期値$u(t=0)=0.5$で$t=0$から$t=1$まで解きたい場合、
using DifferentialEquations
f(u,p,t) = 1.01*u
u0=1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan) #微分がfで初期値がu0の微分方程式du/dt = f(u)を解く。
@time sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8) #Tsit5は解き方を指定
nt = 50
t = range(0.0, stop=1.0, length=nt) #0.0から1.0までのnt点を生成する
for i=1:nt
println("t= $(t[i]), solution: $(sol(t[i])), exact solution $(0.5*exp(1.01t[i]))")
end
とすればよい。なお、厳密解は$u = 0.5 \exp(1.01 t)$なので、それとの比較を出力した。
微分方程式にパラメータがある場合やuがベクトルや行列の場合でも解くことができる。ドキュメントにあるようなx,y,zの三次元空間におけるローレンツ方程式
\frac{dx}{dt} =\sigma(y-x) \\
\frac{dy}{dt} =x(\rho-z)-y \\
\frac{dz}{dt} = xy-\beta z
は、
function parameterized_lorenz(du,u,p,t) #微分dy、関数u、パラメータp、変数t
du[1] = p[1]*(u[2]-u[1])
du[2] = u[1]*(p[2]-u[3]) - u[2]
du[3] = u[1]*u[2] - p[3]*u[3]
end
u0 = [1.0,0.0,0.0] #x=1,y=0,z=0を初期値とする。
tspan = (0.0,1.0) #時間は0から1まで。
p = [10.0,28.0,8/3] #パラメータ三つ
@time prob = ODEProblem(parameterized_lorenz,u0,tspan,p)
と解くことができる。
また、マクロ@ode_defを用いるともう少し数学っぽい書き方ができる。
using ParameterizedFunctions
g = @ode_def LorenzExample begin
dx = σ*(y-x) #dxとあるので、xが微分されるもの。
dy = x*(ρ-z) - y
dz = x*y - β*z
end σ ρ β #パラメータは三つある
u0 = [1.0;0.0;0.0] #x=1,y=0,z=0を初期値とする。
tspan = (0.0,1.0) #時間は0から1まで。
p = [10.0,28.0,8/3] #σ ρ β を設定
@time prob = ODEProblem(g,u0,tspan,p)
この書き方をした方が最適化がされて速いらしい。
[2019年5月１１日追記]マクロ@ode_defを使うためにはusing ParameterizedFunctionsをしなければならなくなったらしく、パッケージモードでadd ParameterizedFunctionsをしてParameterizedFunctions.jlをインストールする必要がある。
非線形関数の最小値
非線形関数$f(x)$を最小化する値$x$を求める問題もよくある問題である。
これは、Optimというパッケージを使えばよい。
http://julianlsolvers.github.io/Optim.jl/stable/#user/minimization/
Juliaを起動して、]でパッケージモードにし、add Optimとすれば入る。
上のリンク先のドキュメントを参考にして、以下に使い方を記す。
まず、
f(x,y) = (1.0 - x)^2 + 100(y - x^2)^2
を最小化するx,yの組を求めてみる。
using Optim
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
x0 = [0.0, 0.0] #初期値を0,0とした。
a1 = optimize(f, x0)
xsol = Optim.minimizer(a1) #関数fを最小化するxの値
println("xsol = $xsol")
fmin = Optim.minimum(a1) #関数fの最小値
println("fmin = $fmin")
minimizerというのが、最小化する変数を得るもので、minimumというのは関数の値を求めるものとなっている。
次に、
f(x) = 2x^2+3x+1
という関数が$x=-2$から$x=1$の範囲で最小となるxを求めてみる。
using Optim
f2(x) = 2x^2+3x+1
a2 = optimize(f2, -2.0, 1.0) #-2から1の間の最小値を探す。
xsol = Optim.minimizer(a2) #関数f2を最小化するxの値
println("xsol = $xsol")
fmin = Optim.minimum(a2) #関数f2の最小値
println("fmin = $fmin")
このように、非線形関数の最小値探索も簡単に計算することができる。