読者です 読者をやめる 読者になる 読者になる

たけし備忘録

自分の好奇心の赴くままに勉強メモ LL系が大好き Python bash Julia C

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
便利な代入3

パイプライン(|>という演算子)による代入です

julia> f |> subs(x, 1)
y + z + 1


便利な代入1の=>による代入が一番わかりやすい気がする。

シンボルから数値への変換

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)
右極限と左極限

sign関数を使って例を書きます。
sign関数とは

という関数です。つまり左極限と右極限が異なる不連続関数です。

julia> limit(sign(x), x=>0, dir="-") # 左極限
-1

julia> limit(sign(x), x=>0, dir="+") # 右極限
1

数値計算

やっぱりなんだかんだ数値計算が出来てなんぼさ。
ということで数値計算の主な方法です。

関数の定義

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
\

めちゃめちゃわかりにくいですがこれは

ということです。

積分
julia> integrate(f, 0, 3)
21

julia> integrate(f, 0, x)
 3
x     2
-- + x  + x
3

julia> integrate(f, x, 3)
   3
  x     2
- -- - x  - x + 21
  3

julia> integrate(f, x, 2x)
   3
7*x       2
---- + 3*x  + x
 3

積分範囲もシンボルが使えます。

積分

同様に重積分も行えます。

julia> @syms x y
(x,y)

julia> integrate(x*y, (y, 0, 1), (x, 0, 1))
1/4

julia> integrate(x^2*y, (y, 0, 1-x^2), (x, -1, 1))
8/105

シンボルを含んだ積分をした後に数値のみの積分も出来ます。

テイラー展開

テイラー展開もできちゃいます。
めんどくさいので略。

第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
limit関数を使った無限級数

limitを使えば部分和を無限級数にできちゃいます。

julia> Sn = summation(1/k^2, (k, 1, n))
harmonic(n, 2)

julia> limit(Sn, n=>oo)
  2
pi
---
 6

harmonic(n, 2)はnをシンボルとした2次の調和級数です。要するにζ(2)です。