15. SymPy in Julia：数式処理

数式処理では，式を数値へ近似する前に，因数分解・展開・代入・微分・積分などを「式のまま」計算する．この回では，JuliaからSymPyを利用する基本操作を学び，最後にJulia製のSymbolics.jlとの違いを整理する．

1. 数式処理とは

Julia ではいくつかの数式処理システムを利用できる．代表的な選択肢の一つが，Pythonで書かれた成熟した数式処理ライブラリSymPyをJuliaから利用する方法である．現在は，共通基盤のSymPyCoreを通して，PyCallを使うSymPy.jlと，PythonCallを使うSymPyPythonCall.jlのほぼ同じJulia向けインターフェースを選べる．

また，Juliaで書かれ，数値計算・モデル生成・高速なJulia関数の生成との連携を重視するSymbolics.jlも発展している．SymPyとSymbolicsは単純な新旧関係ではなく，得意分野が異なる．伝統的で幅広い数式処理にはSymPy，Juliaの数値計算やSciMLとの統合にはSymbolics，という見方が分かりやすい．

今回はSymPyを中心に解説し，最後にSymbolicsを使う場合の違いを紹介しよう．

さて，この SymPy package の使い方だが， SymPyCore Tutorial というチュートリアルがあり，入門用資料として結構よくできているので，本来はこれを一通り読むのが望ましいなあ．できれば読んでおこう．

2. SymPyで記号と式を作る

そして，ここではなるべく簡単に SymPy の使い方を紹介しよう．

パッケージのインストール

個人環境等で SymPy が未インストールの場合は下記のようにしてインストールしておこう．

1using Pkg
2Pkg.add("SymPy")

新しい個人環境でPythonCall系を使いたい場合は，代わりに次を選べる．以後の基本的なコードはほぼ同じである．

1using Pkg
2Pkg.add("SymPyPythonCall")

さて，このパッケージの利用宣言をしよう．

1using SymPy # 数式処理パッケージ

SymPyPythonCallを選んだ場合は，ここだけ次のように読み替える．

1using SymPyPythonCall

そして，SymPy のポイントを先に説明しよう． そのポイントは，数式処理の対象となる変数を先に陽に指定するという仕組みにある． 言い換えると，「指定すると特別な変数」，「指定していない変数は，いつもどおり Julia の普通の変数」であるのだ． そのあたりを例で示そう．

まず，数式処理と関係なく，いつものように普通に Julia で 関数を定義してみよう．

1p(x) = x^2 + 3x + 2

ここで出てくる変数 x は関数定義のために仮に使われたもので特別なものではないので，

1p(x)

と入力すると(それまでに x を使っていなければ) UndefVarError: x not defined なんてエラーが出たりする． また，これまでの知識通り，関数 p の入力を数字に置き換えたり，関数 p をプロットできたりする． 例えば

1p(3)

とすると計算結果として 20 が出てくるし，

1using Plots
2plot(p)

とすると

というプロットが普通に得られる．

しかし，次のように @syms という命令を使って数式処理の対象変数を指定¹してこれを使うとなると話が変わってくる．

1# t という名前の変数を数式処理の対象に！
2@syms t 

まず，この変数を使うと違いが発生することを確認しよう．次の入力をしてみるとすぐわかる．

1eq = p(t)

すると出力が $t^2 + 3t + 2$ という数式で出てきて，「おや，これまでとなにか違うぞ」と感じるわけだ．

何が起きているか少し解説すると，これは数式処理対象として指定した変数 t を使っているので， (p(x) は Julia の関数だが ) p(t) は Julia の関数ではなく，数式として扱われる ことによる違いが見えているのだ．

もう少し正確に言えば，「関数」は入力から出力を計算する手順であり，ここで得られた「数式」はSymPyの式木をJulia側で包んだsymbolic expression objectである．単なる文字列ではないため，式の構造を保ったまま変形できる．

まあこれで eq が数式であることが認識できたので，この数式に対して以下のように，実際に数式処理を施してみよう．

複数の変数

1@syms u v
2
3eq3 = p(v) * (u^2 - 1) * exp(u/v)

$(u^2-1)(v^2+3v+2)e^{u/v}$

3. SymPyによる基本的な式操作

因数分解

1factor( eq )

$(t+1)(t+2)$

展開

1eq2 = expand( eq * (t^2-1) )

$t^4 + 3t^3 + t^2 -3t -2$

代入（式）

代入は大事だ．

1subs(eq, t => t^2)

$t^4 + 3t^2 + 2$

代入（数値）

1subs(eq, t => 3)

$20$

方程式を解く

式を0とみなして，変数について解くには solve を使う．

1solve(eq, t)

$[-2,-1]$

式の簡単化

同じ値を表す式を見通しのよい形へ整理したいときは simplify を使う．

1simplify(sin(t)^2 + cos(t)^2)

$1$

微分

式が複雑なときは手計算より頼りになるかも．

1diff( eq, t )

$2t+3$

不定積分

積分も頼りになるかも．

1integrate( eq, t )

$\displaystyle \frac{t^3}{3}+\frac{3t^2}{2}+2t$

1integrate( 1/(1+t^2), t )

$\mbox{atan}(t)$

定積分

1integrate( eq, (t, -2, PI) )

$\displaystyle \frac{2}{3}+2\pi+\frac{\pi^3}{3}+\frac{3\pi^2}{2}$

SymPy package には特別な定数記号がいくつかある． PI (アルファベット大文字の P と I)は円周率を表す記号であり，oo (アルファベットの o を2つ書く)は無限大を表す記号である．

厳密な記号をJuliaの近似数値へ変換するときはNを使う．

1N(PI)

3.141592653589793

1integrate( 1/(t^2+4), (t, 0, oo))

$\displaystyle \frac{\pi}{4}$

極限

1limit(sin(t)/t, t => 0)

$1$

Taylor展開

1eq4 = exp(t)
2
3series( eq4, t, 0, 5 )

$1 + t + \frac{t^2}{2} + \frac{t^3}{6} + \frac{t^4}{24} + O(t^5)$

総和

たとえば

1@syms n
2eq5 = summation( 1/n^2, (n, 1, t))

として $\mbox{eq5} = \displaystyle \sum_{n=1}^t \left( \frac{1}{n^2} \right)$ を定義してから， 下記のように $t$ に 10 を代入すると，

1subs(eq5, t => 10)

$\frac{1968329}{1270080}$
というように結果が得られるし，また， 下記のように $t$ に $\infty$ を代入すると

1subs(eq5, t => oo)

$\displaystyle \frac{\pi^2}{6}$
というように級数 $\displaystyle \sum_{n=1}^{\infty} \left( \frac{1}{n^2} \right)$ の値が得られる(もちろん，eq5 の定義の際に $t$ としているところに最初から oo を書いても良い)．

線形代数

線形計算ももちろん OK で，

1using LinearAlgebra
2
3A = [ 1 2 3
44 5 6
57 8 t ]
6
7det(A)

$27 - 3t$
と結果が得られる．

Juliaの行列に入っている各要素へ代入したいときは，broadcastを使えばよい．Refで代入規則を一つの値として扱うのがポイントである．

1A15 = subs.(A, Ref((t, 15)))

これで，$t$ が15に置き換わったJuliaの行列が得られる．SymPyはPythonのライブラリをJuliaから利用する仕組みなので，配列全体に対する処理では，このようにJulia側のbroadcastと組み合わせる場面がある．

このように，数式に単になにか変形を加えるというのであれば本当に簡単だ．

4. 数式をJulia関数へ変換する

上の例のように，Julia関数へsymbolic variableを入力すれば数式を作れる．逆に，数式処理で得た式を数値計算へ戻したい場合は lambdify を使う．文字列を組み立てて eval する必要はない．

たとえば，数式 eq2 を変数 t のJulia関数へ変換するには，次のようにする．

1p4 = lambdify(eq2, (t,))

この p4 は通常の数値を入力できる関数なので，

1p4(3)

160

と計算できる．複数の入力変数がある式 eq3 も同様である．引数の順序を明示しておくと安全だ．

1q = lambdify(eq3, (u, v))
2q(2.0, 3.0)

5. 数式をplotする

SymPy の対象である数式は，実はある程度そのまま plot できる．

文法は少しだけいつもと異なるので， SymPyCoreのIntroduction を眺めておくと良い．ただし，Plots.jlの豊富な機能を使うならば，いったん lambdify で通常のJulia関数へ変換する方法が分かりやすい．

1using Plots
2p4 = lambdify(eq2, (t,))
3plot(p4, -3, 3; xlabel="t", ylabel="p4(t)", label="expanded polynomial")

6. Symbolics.jlを利用する場合

インストールと利用宣言

利用環境に入っていない場合は，一度だけインストールする．

1using Pkg
2Pkg.add("Symbolics")

1using Symbolics

変数と式の定義

@variables というマクロを使う．

1# 一つでも，複数でも OK.
2@variables x y t
3
4eq = x^2 + 3x + 2

$\displaystyle x^2 + 3x + 2$

簡単化・展開・代入

simplifyで式を整理できる．expand=trueを指定すると，最初に式を展開してから簡単化する．

1Symbolics.simplify(eq * (x^2-1), expand=true)

$\displaystyle x^4 + 3x^3 + x^2 - 3x -2$

代入にはsubstituteを使う．

1Symbolics.substitute(eq, Dict(x => x^2))

$\displaystyle x^4 + 3x^2 + 2$

配列の各要素へ代入する場合はbroadcastを使う．

1using LinearAlgebra
2
3A = [ 1 2 3
44 5 6
57 8 t ]
6
7A15 = Symbolics.substitute.(A, (Dict(t => 15),))
8Symbolics.value.(A15)

3×3 Matrix:
 1  2   3
 4  5   6
 7  8  15

なお，SymPyのfactorに相当する一般的な因数分解が，Symbolicsの中心的な用途というわけではない．伝統的な代数操作を幅広く試したい場合はSymPy，Juliaの式やモデルを変換・最適化したい場合はSymbolics，と使い分けるとよい．

微分

現在の基本的な書き方では，Differentialで微分作用素を作る．微分は遅延表現されるため，expand_derivativesで実際に展開する．

1D = Symbolics.Differential(x)
2Symbolics.expand_derivatives(D(eq))

$\displaystyle 2x + 3$

積分

以前はSymbolics周辺の記号積分は限られていたが，現在は積分そのものを表すIntegralに加え，SymbolicIntegration.jlやSymbolicNumericIntegration.jlとの連携が用意されている．たとえば，開発中の機能であることに注意しつつ，次のように厳密な原始関数を求められる．

1using Pkg
2Pkg.add("SymbolicIntegration")
3
4using SymbolicIntegration
5SymbolicIntegration.integrate(x^2, x)

$\displaystyle \frac{x^3}{3}$

積分可能な式の範囲や安定性は方法によって異なる．幅広く成熟した積分機能が必要なら，現在もSymPyのintegrateが有力な選択肢である．

Taylor展開とseries

taylorでTaylor展開を計算できる．

1Symbolics.taylor(sin(x), x, 0:5)

$\displaystyle x - \frac{1}{6}x^3 + \frac{1}{120}x^5$

現在のseriesは，係数から形式的な冪級数を組み立てるための機能であり，一般の総和を求めるsumやSymPyのsummationとは役割が異なる．たとえば，係数をsymbolic variableとして次のように表せる．

1@variables a[0:3]
2Symbolics.series(a, x, 0)

数値計算用のJulia関数を作る

Symbolicsではbuild_functionにより，symbolic expressionからコンパイル可能なJulia関数を生成できる．

1f = Symbolics.build_function(eq, x; expression=Val(false))
2f(3)

20

全体として，Symbolicsは旧来の数式処理だけを目標にしたものではなく，数式と高性能な数値計算をつなぐ設計を重視している．SymPyの方が成熟している操作もあれば，Symbolicsが得意とするJulia固有の機能もある．しばらくは両方を目的に応じて使い分けるのがよいだろう．

1eval(Meta.parse("p4(t) = $eq2"))
2p4(3)

160

1eval(Meta.parse("q(u,v) = $eq3"))

1toFunc(str) = eval(Meta.parse(str))
2toFunc("p4(t) = $eq2")

1p5(x) = x^4 + 3x^3 + x^2 - 3x - 2

1@syms x::positive
2simplify(sqrt(x^2))

レポート No.15

以下の課題について，自らの将来のスキルアップにつながるように調査と考察を行い，
     学籍番号-氏名-15.pdf
というファイルとしてレポートを作成し、 webフォーム から教官宛に提出しよう.

なお，レポートを $\TeX$ 等で作成したものを印刷した「紙媒体」を教官に直接手渡す形で提出してもよいが、物質によるレポート提出は常に破損や紛失の可能性があるのであまりお勧めはしないぞ．

課題

1. この授業のこれまでで手計算が必要だった箇所を探し出し，2例ほど上の SymPy を使う形で数式処理してみよう． そして，その結果と自分の手計算の結果を比較しよう．
   
   
2. 自分の持っている微分積分学，線形代数学の教科書を開き，手計算でこれまで計算してきた計算を 2例ほど上の SymPy を使う形で数式処理してみよう． そして，その結果と自分の手計算の結果を比較しよう．