Julia+SymPyの勝手なまとめ
自分の気になったところだけメモしていく形です
基本的にはJuliaのSymPyパッケージ内に入っているtutorial.mdというファイルから抜粋する形です。
SymPy の導入
まずはjuliaを立ち上げてSymPyのパッケージをインストールしましょう。
PyCallに依存しているので入れてない方はPyCallもインストールしましょう。
$ julia _ _ _ _(_)_ | A fresh approach to technical co (_) | (_) (_) | Documentation: http://docs.julia _ _ _| |_ __ _ | Type "?help" for help. | | | | | | |/ _ | | | | |_| | | | (_| | | Version 0.4.3 (2016-01-12 21:37) _/ |\__'_|_|_|\__'_| | Official http://julialang.org/ r |__/ | x86_64-w64-mingw32 julia> Pkg.add("PyCall") ...........なんかいろいろ出る julia> Pkg.add("SymPy") ...........なんかいろいろ出る julia>
<*追記*>
依存がPyCallとSymPyだけじゃありませんでした。
JuliaのPyCallとSymPyパッケージを入れた環境に、Python(2でも3でもいい)をインストールし、かつPythonのSymPyライブラリをインストールする必要があります。
python.jp
こちらのPython日本語ホームページからダウンロードできます。バージョンは2でも3でもお好きな方で大丈夫です。
Pythonのインストールが終わったら、次はPythonのSymPyをインストールする必要があります。
Pythonは2.7.9以上のバージョン(3を含む)ならば、PIPというパッケージ管理ツールが同梱していますのでそのままお使いいただけます。
$ python -m pip install sympy
で基本的にはインストールできます。
- 上記でインストール出来ない場合(Windowsのみ記載)
上記の方法でインストール出来ない場合、次の非公式パッケージがダウンロードできるサイトで.whlというファイルをダウンロードしましょう。
Python非公式パッケージ
現時点(2016/1)ではsympy-0.7.6.1-py2.py3-none-any.whlが最新のようです。
この.whlファイルを、ターミナルから操作しやすいディレクトリに入れておきましょう。
今回の場合、C:\test\sympy-0.7.6.1-py2.py3-none-any.whlとして保存しておくことにします。
$ pwd # 現在地 C:\test $ ls sympy-0.7.6.1-py2.py3-none-any.whl $ python -m pip install sympy-0.7.6.1-py2.py3-none-any.whl (色々出る)
で完了です。
- 上記でインストール出来ないかつ不安だからpython2と3の両方インストールしちゃったよの場合(Windowsのみ記載)
pythonはデフォルトで2と3が共存できます(3の方が新しいですが2への後方互換性が無いです)。
そのような場合 pyコマンドで2と3を明示的に切り替えることが出来ます。
# Python2を使う場合 $ py -2 -m pip install sympy-0.7.6.1-py2.py3-none-any.whl # Python3を使う場合 $ py -3 -m pip install sympy-0.7.6.1-py2.py3-none-any.whl
シンボルの生成
SymPyのJuliaでのラッパは色々なシンボルの生成の仕方があります。
1つずつシンボルを生成
julia> x = Sym("x") # Sym関数での生成 x julia> @syms x # symsマクロでの生成 (x,) julia> x x
複数のシンボルを一気に生成
julia> x, y, z = symbols("x, y, z") # symbols関数での生成 (x,y,z) julia> a, b, c, d, e = symbols("a:e") # ":"というイテレータで一気に生成 (a,b,c,d,e) julia> @syms x y z # symsマクロで生成 (x,y,z)
マクロでの生成が一番手軽でよさそう
特殊な定数 PI, E, oo
PIは円周率
Eは自然定数の底(ネイピア数)
ooは無限大です
julia> PI,E,oo
(pi,E,oo)
juliaはデフォルトでpiという定数がありますがPIとは異なります。
julia> [asin(1), asin(Sym(1))] 2-element Array{SymPy.Sym,1} [1.5707963267949] [ ] [ pi ] [ -- ] [ 2 ]
代入
シンボルで定義した式に、数値またはシンボルを代入します。
Julia固有の色々な代入の仕方があります。
subs関数による代入
前準備
julia> @syms x y (x,y) julia> f = x^2+2x+1 2 x + 2*x + 1
シンボルxにシンボルyを代入
julia> subs(f, x, y) 2 y + 2*y + 1
シンボルxに数値1を代入
julia> subs(f, x, 1) 4
便利な代入
前準備
julia> @syms x y z
(x,y,z)
julia> f=x+y+z
x + y + z
便利な代入1
これは便利な気がする
julia> f(x=>1, y=>pi) z + 1 + pi
便利な代入2
代入の順序がよくわかんない
julia> f(1, pi) y + 1 + pi
シンボルから数値への変換
evalf関数とN関数の二つの方法があります。
NはBigFloatという型らしいです。よく知りません。
evalfはお馴染みです。
julia> evalf(PI) 3.14159265358979 julia> N(PI) 3.141592653589793 julia> evalf(PI, 30) # 桁数指定 30桁 3.14159265358979323846264338328 julia> N(PI, 30) # 桁数じゃなかった。何者だろうこの引数は 3.1415926535897932384626433832793
方程式の解
solve関数を使います。
前準備
julia> @syms x (x,) julia> f=x^2+x+1 2 x + x + 1
解を得る1 x^2+x+1=0
julia> solve(f) 2***element Array{SymPy.Sym,1} [ ___ ] [ 1 \/ 3 *I] [- - - -------] [ 2 2 ] [ ] [ ___ ] [ 1 \/ 3 *I] [- - + -------] [ 2 2 ]
解を得る2 x^2+x+1=1
julia> solve(Eq(f, 1)) 2-element Array{SymPy.Sym,1} [-1] [ ] [0 ]
解を得る3 x^2+x+1=0
julia> solve(f⩵0) 2-element Array{SymPy.Sym,1} [ ___ ] [ 1 \/ 3 *I] [- - - -------] [ 2 2 ] [ ] [ ___ ] [ 1 \/ 3 *I] [- - - -------] [ 2 2 ] # ※注意!! f==0と書いていますがこれは \Equalの後にTabを押して出る特殊記号です。
極限
極限はlimit関数で操作できます。
基本の計算方法
julia> @syms x h (x,h) julia> limit(sin(x)/x, x, 0) 1 julia> limit(sin(x)/x, x=>0) 1 julia> limit((1+1/x)^x, x=>oo) E julia> limit((sin(x+h)-sin(x))/h, h=>0) cos(x)
数値計算
やっぱりなんだかんだ数値計算が出来てなんぼさ。
ということで数値計算の主な方法です。
関数の定義
Juliaでの関数の定義はいろいろなやり方があります。
まずSymPyを使うとシンボルを定義します。
シンボルはSymPy.Sym型の変数ですのでJuliaの性質である多重ディスパッチが使えます。
多重ディスパッチとは、同名の関数、例えばtestという関数を
function test(x::Int) return x end function test(x::Float64) return 2x end
として同名で、引数の型の違いによって出力や内容を変える性質のことです。
多重ディスパッチを使うと、
@syms x f(x) = x^2+x+1
という風に定義すると
julia> f(x) 2 x + 2*x + 1 julia> f(1) 4 julia> f(1.0) 4.0 julia> f(1:10) ERROR: MethodError: `*` has no method matching *(::UnitRange{Int64}, ::UnitRange{Int64}) Closest candidates are: *(::Any, ::Any, ::Any, ::Any...) *{T}(::Bidiagonal{T}, ::AbstractArray{T,1}) *(::Number, ::AbstractArray{T,N}) ... in power_by_squaring at intfuncs.jl:78 in f at none:1
1:10を引数にするとエラーが出ました。
そこでArray型にも対応できるように
f(x) = x.^2+2x+1
とします。 x^2がx.^2になっただけです。
このようにすれば
julia> f(1:10) 10-element Array{Int64,1}: 4 9 16 25 36 49 64 81 100 121
となります。
行列も引数に与えられます。
SymPyのままでも
@syms x f = x.^2+2x+1
としても
julia> f(x) 2 x + 2*x + 1 julia> f(1) 4 julia> f(1.0) 4.00000000000000 julia> f(1:10) ERROR: MethodError: `convert` has no method matching convert(::Type{SymPy.Sym}, ::UnitRange{Int64}) This may have arisen from a call to the constructor SymPy.Sym(...), since type constructors fall back to convert methods. Closest candidates are: SymPy.Sym(::Any) SymPy.Sym(::Any...) call{T}(::Type{T}, ::Any) ... in subs at C:\Users\takeshi\.julia\v0.4\SymPy\src\subs.jl:46 in call at C:\Users\takeshi\.julia\v0.4\SymPy\src\subs.jl:115
とこのように関数のように使うことはできます。
しかし1:10などの配列や行列などは代入できません。
この理由としては
f = x.^2+2x+1
と定義するとSymPy.Sym型であり
f(x) = x.^2+2x+1
と定義するとこれは多重ディスパッチによって無名関数になるからです。
型を確認してみましょう。
区別がつきやすいように
@syms x f(x) = x.^2+2x+1 g = x.^2+2x+1
としておきます。
julia> typeof(f)
Function
julia> typeof(g)
SymPy.Sym
fがFunction型、gがSymPy.Sym型になっています。
またさらにfの場合
julia> typeof(f(x)) SymPy.Sym julia> typeof(f(1)) Int64 julia> typeof(f(1.0)) Float64 julia> typeof(f(1:10)) Array{Int64,1}
となります。
gの場合は
julia> typeof(g(x)) SymPy.Sym julia> typeof(g(1)) SymPy.Sym julia> typeof(g(1.0)) SymPy.Sym
となっていずれもSympy.Sym型になります。
なかなかこれが使えます。
微分
当たり前ですがPythonのSymPyと同様に記号微分が出来ます。
julia> @syms x (x,) julia> f(x) = x.^2+2x+1 f (generic function with 1 method) julia> diff(f(x)) 2*x + 2 julia> diff(f(x), x) 2*x + 2 julia> diff(f) 2*x + 2
diff関数は多重ディスパッチで
diff(f:SymPy.Sym)だけでなく
diff(f::Function)も定義されているので、無名関数として定義したfを入れても微分してくれます。
diff関数のオプション
diff関数は任意の回数だけ微分することができます。
以下に例を書きます。
julia> f(x) 2 x + 2*x + 1 julia> diff(f(x), x, 0) # 0階微分 2 x + 2*x + 1 julia> diff(f(x), x, 1) # 1階微分 2*x + 2 julia> diff(f(x), x, 2) # 2階微分 2
もちろんfのままでも渡せます。が、少し文法が異なります。
要するにxが消えます。無名関数のまま渡すので、引数の情報がいらないわけです。
julia> f(x) 2 x + 2*x + 1 julia> diff(f, 0) # 0階微分 2 x + 2*x + 1 julia> diff(f, 1) # 1階微分 2*x + 2 julia> diff(f, 2) # 2階微分 2
偏微分
2つのシンボルで定義された関数は偏微分できます。
さっそく例を出していきましょう。
julia> @syms x y (x,y) julia> f(x, y) = x^2 + y^2 + x*y f (generic function with 2 methods) julia> diff(f(x,y), x) # xについて偏微分 2*x + y julia> diff(f(x,y), y) # yについて偏微分 x + 2*y julia> diff(f(x,y), x, y) # x,yについて偏微分 1 julia> diff(f(x,y)) ERROR: PyError (:PyObject_Call) <type 'exceptions.ValueError'> ValueError('\nSince there is more than one variable in the expression, the\nvariable(s) of differentiation must be supplied to differentiate x**2\n+ x*y + y**2',) File "C:\Python27\lib\site-packages\sympy\core\function.py", line 1638, in diff return Derivative(f, *symbols, **kwargs) File "C:\Python27\lib\site-packages\sympy\core\function.py", line 987, in __new__ must be supplied to differentiate %s''' % expr)) [inlined code] from C:\Users\takeshi\.julia\v0.4\PyCall\src\exception.jl:81 in pycall at C:\Users\takeshi\.julia\v0.4\PyCall\src\PyCall.jl:361 in fn at C:\Users\takeshi\.julia\v0.4\PyCall\src\conversions.jl:194 in call_sympy_fun at C:\Users\takeshi\.julia\v0.4\SymPy\src\SymPy.jl:262 in sympy_meth at C:\Users\takeshi\.julia\v0.4\SymPy\src\SymPy.jl:268 in diff at C:\Users\takeshi\.julia\v0.4\SymPy\src\math.jl:190 julia> diff(f) # 無名関数を偏微分 2*x + 2
積分
不定積分
julia> f(x) 2 x + 2*x + 1 julia> integrate(f(x), x) 3 x 2 -- + x + x 3 julia> integrate(f(x)) 3 x 2 -- + x + x 3 julia> integrate(f) 3 x 2 -- + x + x 3
より一般的な書き方
julia> @syms x n (x,n) julia> integrate(x^n, x) /log(x) for n = -1 | | n + 1 <x |------ otherwise |n + 1 \
めちゃめちゃわかりにくいですがこれは
ということです。
第n部分和と無限級数
基本的な使い方
まずは基本からです。summation関数を使います。
julia> @syms k n (k,n) julia> summation(2^k, (k, 1, n)) # 第n部分和 n + 1 2 *** 2 julia> summation(1/k^2, (k, 1, oo)) # 無限級数(この例はゼータ関数) 2 pi --- 6