SymPy/Julia
https://www.math.csi.cuny.edu/~verzani/tmp/julia/symbolic.html
1. Julia で数式処理
コンピュータによる計算処理で，３ 通りくらいのアプローチがある。
- Mathenatica，Maple, Sage に代表される数式処理
- MATLAB, R による数値処理
- Python, Haskel による一般化された汎用処理
Julia は 数値計算処理のために，MATLAB や R の代替として存在する。Julia は PyCall を使うことによって，Python およびその SymPｙ パッケージへの連携を提供した。
Julia で SymPy パッケージを使えるようにする方法については，ここでは述べない。単に，Python を動くようにすることと，SymPy パッケージを使えるようにした後，Pkg.add("SymPy") するだけである。
2. 概略
数式処理をするために必用なのは，SymPy を使うことを宣言する（using SymPy）と，そこで使う変数を宣言すること（@syms x など）である。
using SymPy
@syms x
(x,)
これによって宣言（定義）される ｘ は，Julia における変数名ではなく，以下に引き続く処理の記述で使用される記号変数（Sym オブジェクト）である。
x の型(typeof) は Sym である。
x |> typeof |> println
Sym
Sym オブジェクト（以下では単に変数と呼ぶ）は，何通りかの方法で宣言（定義）することができるが，@syms マクロ（または @vars マクロ）を使うのがもっと便利であろう。h, y の 2 変数を定義するには以下のようにする。
変数名は空白あるいはカンマで区切って複数記述できる。
@syms h y
(h, y)
なお，実際の使用例を後で示すが，特定の条件をもつ変数として宣言することもできる。
@syms a::positive b::real c::integer
(a, b, c)
複数の条件を持つ場合は @vars を使い，条件をタプルで指定する。
@vars x y z real=true w nonnegative=true v (positive=true, real=true)
(x, y, z, w, v)
このようにして定義した変数を使用して，任意の数式を記述することができる。
たとえば x^2 -2x +2 のような数式はそのまま指定できる。
この記事を通じて，数式のあとに "|> typeof" あるいは "|> println" としているのは，左側の数式に対して typeof(数式)，さらに printiln(処理された数式） を実行し，最終的な結果を表示する（REPL の場合）ことを示しているだけである。
ユーザが実際に SymPy を使用する環境によってはそれらは不要，あるいはないほうがよい，更にはあっては困るということもあるであろうから，適宜省略・削除するようお願いしたい。（最後のいくつかの |> を省略したらどうなるか，事前に把握しておいて欲しい）。
以下では，定義した数式の型が何であるかを表示している。なれてくれば，このような確認は不用であろう。
x^2 - 2x + 2 |> typeof |> println
Sym
以下では，関数名 f で普通の Julia 関数を定義している。
f(x) = exp(-x^2/2)
f (generic function with 1 method)
普通通り，x = 3 のときの関数値は f(3) で得られる。
f(3) |> println
0.011108996538242306
しかし，f(x) は Sym 型を持つ数式である（typeof(f(x)) = Sym）。
f(x) |> println
exp(-x^2/2)
typeof(f(x)) |> println
Sym
f(x) が Sym オブジェクトであるということは，そのオブジェクト中の変数を特定の値として評価できるということである。
以下は，f(x) 中の x を 1 にするとどのような数式になるかを示すものである。
subs(f(x), x, 1) |> println
exp(-1/2)
f(x)(x => 1) |> println
exp(-1/2)
f(x) が表す式にある x を 1 で置き換えると exp(-1/2) になるということである。
subs() 関数によって式の一部が置き換えられた（代入された）ものも，前と同じ Sym オブジェクトなので，その式に別のものを代入することもできる。
以下の示すのは，f(x) において， x = 1 を代入 (subs) した結果の数値 (float) を表示する (println) という一連の処理を示している。
f(x) |> subs(x => 1) |> float |> println
0.6065306597126334
同じ処理をするのには何通りかの方法があることがあるが，今の場合は以下のほうがより簡潔かもしれなお。
f(x)(x => 1) |> float |> println
0.6065306597126334
多くの場合，最終的な（人間が見てわかる）表示形式にするのは一番最後になるであろう。
3. 方程式
SymPy は代数演算として，簡約化，因数分解，方程式の解など多くのタスクをこなすことができる。
しかし，数式処理をするために必用とされるテクニックはそんなに多くはない。
例えば，以下のように記述すれば，数式は自動的に簡約化される（最後の 「|> println」 は，処理された結果を画面表示するとそのようになるということであるが，環境によってはそれがないほうがわかりやすい数式表現がされることもあろう）。
(x+1) + (x+2) + (x+3) |> println
3*x + 6
書かれた数式をもう少しわかりやすい（人と場合によるのだが）形式で表現してほしいこともある。
例えば，以下の式は人によっては別の形式で表現してほしいかもしれない。現役数学者（中高生？）なら，我慢できないかもしれない。
sin(x)^2 + cos(x)^2 |> println
sin(x)^2 + cos(x)^2
このよう場合，simplyfi() 関数により有名すぎる公式「x が何であろうと（有効範囲はあるが），sin(x) と cos(x) の二乗の和和は 1」簡約化された解を得られる。
simplify(sin(x)^2 + cos(x)^2) |> println
1
sin(x)^2 + cos(x)^2 |> simplify |> println
1
factor() 関数により，数式の因数分解もできる。
factor(x^2 + 2x + 1) |> println
(x + 1)^2
因数分解された形式で表示された数式は，expand() 関数によって展開された数式にすることができる。
expand( (x+1)*(x+2)*(x+3) ) |> println
x^3 + 6*x^2 + 11*x + 6
4. 方程式を解く
多項式の実数解を数式解として求めるには real_roots() 関数を使う。
real_roots(x^2 + x - 1) |> println
Sym[-sqrt(5)/2 - 1/2, -1/2 + sqrt(5)/2]
もう少し一般的にいえば，多項式に限らず方程式の数式解を得るには solve() 関数を使う。
solve(x^2 + x -1, x) |> println
Sym[-1/2 + sqrt(5)/2]
得られるのは数式解なので，数値解を得るには float() 関数を適用すればよい。
solve(x^2 + x -1, x) |> float |> println
[0.6180339887498949]
三角関数方程式は無数の解を持つが，SymPy の sin では -90度 〜 90度の範囲の解を求める。
solve(sin(x)^2 - (1/2)^2) * (180 /pi) |> println
Sym[30.0000000000000, 150.000000000000, 210.000000000000]
式中に変数が 2 つ以上ある場合に，どの変数についての解を求めるかは第 2 引数で指定する。
fx = x^2 + y^2 - 1
out = solve(fx, y)
println(out)
Sym[-sqrt(1 - x^2), sqrt(1 - x^2)]
方程式は二次式なので，2 つの数式解を含むベクトルが返される。x = 1/2 のときの値は，各要素について x = 1/2 を代入して，float にする。
out[1](x => 1/2) |> float |> println
out[2](x => 1/2) |> float |> println
-0.8660254037844386
0.8660254037844386
以下のようにすればベクトルとしてまとめて表示できる。
println([o(x => 1/2) |> float for o in out])
[-0.8660254037844386, 0.8660254037844386]
println(map(arg -> arg(x => 1/2) |> float, out))
[-0.8660254037844386, 0.8660254037844386]
連立方程式も solve で解くことができる。 x + y - 2a = 0 x - y + b = 0 を x, y について解く。 連立方程式と，どの変数について解くかの指定をベクトルで指定する。
solve([x + y - 2a, x - y + b], [x, y]) |> println
Dict{Any, Any}(y => a + b/2, x => a - b/2)
別法として，Julia で数値解を得るときと同じように記述することもできる。
x + y = 2a x - y = -b
係数行列 [1 1; 1 -1] と定数ベクトル [2a, -b] と \ 演算子を使う。
Sym[1 1; 1 -1] \ Sym[2a, -b] |> println
Sym[a - b/2, a + b/2]
4.1. 課題
4.1.1. 問題 1.
(x-1)(x-2)(x-3)(x-4)(x-5) において，x^3 の係数はいくつか。
eq = expand((x-1)*(x-2)*(x-3)*(x-4)*(x-5))
eq |> println
x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120
目で見るだけで解は分かるが，係数を取り出す object.coeff(...) メソッドを使えばよい。
eq.coeff(x^3) |> println
85
4.1.2. 問題 2. Horner 法
((1x + 2)x + 3)x + 4 を展開せよ。
((2x + 3)x + 4)x + 5 |> expand |> println
2*x^3 + 3*x^2 + 4*x + 5
Horner 法は，多項式 2 x^3 + 3 x^2 + 4 x + 5 を記述する方法の一つで，((2 x + 3)x + 4)x + 5 のように ( ) を使って記述する。
((2 x + 3) x + 4) x + 5 では，加算が 3 回，乗算が 3 回であるのに対し
2 x^3 + 3 x^2 + 4 x + 5 では，加算は同じく 3 回であるが，乗算は（^ も加え） 6 回（3, 2, 1 回）である。 もっとも，最近のコンピュータの性能から見ればたいした違いではないだろうし，Julia のようなコンパイラ言語では数式の解釈時に Hoener 法を使って最適化するようなこともあるのかもしれない。
4.1.3. 問題 3.
x^3 - 6x^2 + 11x - 6 = 0 を解き，最も大きい解を求めよ。
以下で，N() 関数は，適切な数値表現への変換関数である。float() 関数は引数が整数でも実数化してしまう（この場合には，実害はないのだが）。
solve(x^3 - 6x^2 + 11x - 6) |> N |> println
[1, 2, 3]
目で見ても分かるが，最大値を求める関数 maximum() を使えばよい。
solve(x^3 - 6x^2 + 11x - 6) |> maximum |> println
3
4.1.4. 問題 4.
直角三角形において，直角を挟む二辺が 3cm と 4cm の場合に，斜辺の長さ x を求めよ。
この場合，x は非負の実数なので，@syms で変数を定義するときに「変数::変数の特性」の書き方を使い，nonnnegative であるという条件を付けておくと不要な解が含まれなくなる。
@syms x::nonnegative
solve(x^2 - (3^2 + 4^2)) |> N |> println
[5]
なお，Julia ではベクトルは要素が １ 個の場合でもベクトルなので，スカラーとして取り出すときは 「ベクトル[1]」 のように添え字を指定しないと「[要素]」となってしまうので注意が必要である。
要素数が 1 個の場合と，2 個以上の場合の両方に対応するには以下のようにする。
a = solve(x^2 - (3^2 + 4^2))
length(a) == 1 ? println(N(a[1])) : println(a)
5
@syms y
b = solve((y + 3) * (y - 5))
length(b) == 1 ? println(N(b[1])) : println(N(b))
[-3, 5]
4.1.5. 問題 5.
x^4 - 10x^3 + 32x^2 - 38x + 15 = 0 を解き，最も大きい解を求めよ。
問題 5 と同じである。
@syms x
solve(x^4 - 10x^3 + 32x^2 - 38x + 15) |> N |> println
[1, 3, 5]
solve(x^4 - 10x^3 + 32x^2 - 38x + 15) |> maximum |> println
5
いまの場合，4 次式なので解は 4 個あるはずだが，solve() では 3 個しか表示されなかった。
solve では重解は複数として数えないのである。factor を使えば詳細が分かる。
factor(x^4 - 10x^3 + 32x^2 - 38x + 15) |> println
(x - 5)*(x - 3)*(x - 1)^2
この場合は，x = 1 が二重解であることが分かる。
このような場合には，roots() 関数を使えば重解も複数として数える。解は辞書形式で返される。
roots(x^4 - 10x^3 + 32x^2 - 38x + 15) |> println
Dict{Any, Any}(3 => 1, 1 => 2, 5 => 1)
この場合は，x=3, x=5 は 1つずつ，x=1 は 2 個であることが分かる。
実数解だけを求める必要がある場合は，real_roots() 関数を使うことができる。
以下の p は 8 次式なので，実数解，虚数解を合わせて 8 個ある。
p = (x - 3)^2 * (x - 2) * (x-1) * x * (x + 1) * (x^2 + x + 1)
expand(p) |> println
x^8 - 7*x^7 + 13*x^6 + 2*x^5 - 11*x^4 - 13*x^3 - 3*x^2 + 18*x
real_roots() によれば，そのうち 6 個は実数解であることがわかる。
real_roots(p) |> N |> println
[-1, 0, 1, 2, 3, 3]
p のすべての解は roots() 関数でわかる。
虚数解は共役複素数 -1/2 ± sqrt(3)*I/2 （2 個）である。
roots(p) |> println
Dict{Any, Any}(-1 => 1, 3 => 2, 1 => 1, -1/2 - sqrt(3)*I/2 => 1, -1/2 + sqrt(3)*I/2 => 1, 0 => 1, 2 => 1)
4.1.6. 問題 6.
e^x = x^4 を解け。
Julia で，ℯ は REPL で \euler[tab] ，⩵ は \Equal[tab] として入力する。（ここからコピーペーストしても良い）
solve(ℯ^x ⩵ x^4) |> println
Sym[-4*LambertW(-1/4), -4*LambertW(1/4), -4*LambertW(-1/4, -1)]
⩵ が入力しづらければ Eq() 関数を使って以下のように書くこともできる。
solve(Eq(ℯ^x, x^4)) |> println
Sym[-4*LambertW(-1/4), -4*LambertW(1/4), -4*LambertW(-1/4, -1)]
LambertW() Lambert の W 関数というものだそうだ。float() を通せばその値がわかる。
solve(Eq(ℯ^x, x^4)) |> float |> println
[1.4296118247255556, -0.8155534188089607, 8.6131694564414]
N() では BigFloat として多くの桁数が表示される
solve(Eq(ℯ^x, x^4)) |> N |> println
BigFloat[1.42961182472555561227524441622361901326636222030404817451048179435868561038458, -0.8155534188089606577727273253085594805974099088406385398936250403272774562683745, 8.613169456441398596676396600371925502144825941277910783538046030885450338602536]
方程式の解として正しいか，evalf() を使って検算してみる。
3 番目の解では，左辺と右辺の数値は表示されている範囲では等しいが，全体としては誤差の範囲で一致しないということであろう。
for res in solve(Eq(ℯ^x, x^4))
left, right = ℯ^res.evalf(), (res^4).evalf()
println("left = $left right = $right right == left = $(right == left)")
end
left = 4.17707743900016 right = 4.17707743900016 right == left = true
left = 0.442394430203601 right = 0.442394430203601 right == left = true
left = 5503.66468907673 right = 5503.66468907672 right == left = false
Lambert の W 関数は Julia の LambertW パッケージに用意されている。
# import Pkg; Pkg.add("LambertW")
using LambertW
solve(Eq(ℯ^x, x^4)) |> println
Sym[-4*LambertW(-1/4), -4*LambertW(1/4), -4*LambertW(-1/4, -1)]
typeof(sympy.LambertW)
PyObject
PyCall.PyObject の LambertW() 関数を LambertW パッケージの lambertw() 関数で置き換える。
res = [-4*lambertw(-1/4), -4*lambertw(1/4), -4*lambertw(-1/4, -1)]
res |> float |> println
[1.4296118247255556, -0.8155534188089607, 8.6131694564414]
検算してみる。
res .^ 4 |> float |> println
[4.177077439000158, 0.44239443020360086, 5503.664689076726]
ℯ .^ res |> float |> println
[4.177077439000158, 0.44239443020360075, 5503.664689076728]
