15. SymPy in Julia: 数式処理
Photo by Jeswin Thomas on Unsplash
数式処理 超入門
Julia ではいくつかの数式処理システムを利用できる． その中でもっともよく使われているものが Python で書かれた SymPy を Julia でも利用できるようにした SymPy package だろう．
また，Julia 言語で書かれているうえに最新の設計を用いた，(ただし目的はどちらかというと数値計算の補助の様子な) Symbolics package というものがだいぶ進歩してきているので，これからはこちらを使う方が良いかもしれない．
今回は，SymPy を中心に解説し，最後に，Symbolics の場合の違い等を付加的に解説しよう．
さて，この SymPy package の使い方だが， SymPyCore Tutorial というチュートリアルがあり，入門用資料として結構よくできているので，本来はこれを一通り読むのが望ましいなあ．できれば読んでおこう．
はじめに
そして，ここではなるべく簡単に SymPy の使い方を紹介しよう． このパッケージは，VIORAS環境だと既にインストール済みだ．
SymPy パッケージのインストール (VIORAS 環境では不要！)
個人環境等で SymPy が未インストールの場合は下記のようにしてインストールしておこう．
1using Pkg
2Pkg.add("SymPy")
さて，このパッケージの利用宣言をしよう．
1using SymPy # 数式処理パッケージ
初めて using SymPy を行ったときに、足りないものがあればインストールされるので、その場合は少し待とう．
そして，SymPy のポイントを先に説明しよう． そのポイントとは，
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してこれを使うとなると話が変わってくる．
1# t という名前の変数を数式処理の対象に！
2@syms t
まず，この変数を使うと違いが発生することを確認しよう．次の入力をしてみるとすぐわかる．
1eq = p(t)
すると出力が $t^2 + 3t + 2$ という数式で出てきて，「おや，これまでとなにか違うぞ」と感じるわけだ．
何が起きているか少し解説すると，これは数式処理対象として指定した変数 t を使っているので，
(p(x) は Julia の関数だが ) p(t) は Julia の関数ではなく，数式として扱われる ことによる違いが見えているのだ．
まあもう少し書くと，ここでは，「関数は計算手順を書いたもの」であるが「数式は(特殊な)文字列」だ，と思っておけば良い．
まあこれで eq が数式であることが認識できたので，この数式に対して以下のように，実際に数式処理を施してみよう．
変数が複数の場合．とくに問題ない．
1@syms u v
2
3eq3 = p(v) * (u^2 - 1) * exp(u/v)
$(u^2-1)(v^2+3v+2)e^{\frac{u}{v}}$
因数分解
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$
微分
式が複雑なときは手計算より頼りになるかも．
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つ書く)は無限大を表す記号である．
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$
と結果が得られる．
上のサンプル中の行列 $A$ の中の変数 $t$ に式や数字を代入しようとするとだいぶ遠回りの手続きが必要になる． これは Julia と Python での行列の扱い方を SymPy が合わせようとしてないから，ということになるだろう． SymPy が異なる言語をつなぐ仕組みである以上，どうしてもこうした「齟齬」があれこれとある． こうした点が気になる人は SymPy ではなく，Julia で書かれた Symbolics package を利用するのが良いだろう．
まあ話を戻すと、このように，数式に単になにか変形を加えるというのであれば本当に簡単だ．
数式を関数に戻す/変換するには？
上の例のように，関数から数式を作るのは単に「特別な変数を代入」するだけでよいので簡単だった．
では逆はどうだろう？
せっかく数式処理して得た結果の数式を関数に直して計算に使いたい，というのは自然な要求だ．
それは簡単で，たとえば上の数式 eq2 を 変数 t をもらって t の4次式に基づいた計算値を返す関数に変換するには，
1eval(Meta.parse("p4(t) = $eq2"))
とすれば良い2．
これで数式 eq2 を変換した関数 p4(t) を定義したことになる(ただしこの p4(t) の中の t は SymPy と関係なく，単なる普通の変数扱い)．
これで得られた p4(t) は確かに上で求めた 4次多項式の値を計算する関数になっている．
複数の入力変数があるときでも同様で，たとえば上の数式 eq3 に対しては
1eval(Meta.parse("q(u,v) = $eq3"))
とすれば，関数 q(u,v) を定義したことになる．
ちなみに，関数に変換するたびに毎回 eval(Meta.parse... と書くのは大変だという人は，using SymPy の直後にいつも
1toFunc(str) = eval(Meta.parse(str))
という感じの定義をつけておけば良い．
こうしておけば，上の p4(t) を定義するにも
1toFunc("p4(t) = $eq2")
と書けばよい．少しだけ楽ができるし，構文を覚えておかなくても良いね．
さて，変換した結果は通常の関数だから，たとえば p4 を使って t = 3 の時のこの4次多項式の値は
1p4(3)
160
と普通に計算できる．
これが普通の関数であることは，次のように通常の関数と比較することでも確認できる． まず，この4次多項式を普通の関数で作る．
1p5(x) = x^4 + 3x^3 + x^2 -3x -2
そして，同じ入力値に対してのそれぞれの動作を @code_lowered マクロで比較してみれば良い．
1@code_lowered p4(3)
CodeInfo(
1 ─ %1 = Main.:-
│ %2 = Main.:-
│ %3 = Main.:+
│ %4 = Main.:^
│ %5 = builtin Core.apply_type(Base.Val, 4)
│ %6 = dynamic (%5)()
│ %7 = dynamic Base.literal_pow(%4, t, %6)
│ %8 = Main.:*
│ %9 = Main.:^
│ %10 = builtin Core.apply_type(Base.Val, 3)
│ %11 = dynamic (%10)()
│ %12 = dynamic Base.literal_pow(%9, t, %11)
│ %13 = dynamic (%8)(3, %12)
│ %14 = Main.:^
│ %15 = builtin Core.apply_type(Base.Val, 2)
│ %16 = dynamic (%15)()
│ %17 = dynamic Base.literal_pow(%14, t, %16)
│ %18 = dynamic (%3)(%7, %13, %17)
│ %19 = Main.:*
│ %20 = dynamic (%19)(3, t)
│ %21 = dynamic (%2)(%18, %20)
│ %22 = dynamic (%1)(%21, 2)
└── return %22
)
1@code_lowered p5(3)
CodeInfo(
1 ─ %1 = Main.:-
│ %2 = Main.:-
│ %3 = Main.:+
│ %4 = Main.:^
│ %5 = builtin Core.apply_type(Base.Val, 4)
│ %6 = dynamic (%5)()
│ %7 = dynamic Base.literal_pow(%4, x, %6)
│ %8 = Main.:*
│ %9 = Main.:^
│ %10 = builtin Core.apply_type(Base.Val, 3)
│ %11 = dynamic (%10)()
│ %12 = dynamic Base.literal_pow(%9, x, %11)
│ %13 = dynamic (%8)(3, %12)
│ %14 = Main.:^
│ %15 = builtin Core.apply_type(Base.Val, 2)
│ %16 = dynamic (%15)()
│ %17 = dynamic Base.literal_pow(%14, x, %16)
│ %18 = dynamic (%3)(%7, %13, %17)
│ %19 = Main.:*
│ %20 = dynamic (%19)(3, x)
│ %21 = dynamic (%2)(%18, %20)
│ %22 = dynamic (%1)(%21, 2)
└── return %22
)
この2つの結果を比べると，入力変数名が異なるほかは全く同じであることがわかるだろう．安心だね．
Plot
SymPy の対象である数式は，実はある程度そのまま plot できる．
文法は少しだけいつもと異なるので， SymPy package Plotting を眺めておくと良い． この資料に載っている以上に細かい plot をしたければ，上に書いてあるように，いったん普通の関数に変換してから plot すると良いだろう．
Symbolics package を利用する場合
インストール
センターが用意している環境 VIORAS にも入っていないので，使うならばインストールが必要だ． 以下のように一回だけやれば良い．
1using Pkg
2Pkg.add("Symbolics")
3# Symbolics package 本体のインストール
4Pkg.add("SymbolicIntegration")
5# 積分計算用サブパッケージのインストール. ちょっと時間がかかる．
利用するには
1using Symbolics
2using SymbolicIntegration # 積分計算も行う場合
とすれば良い．
なお、積分機能パッケージ SymbolicIntegration に対して
初めて using SymbolicIntegration を行うときに不足なものをダウンロードするが、これにはかなり時間がかかる(10分以上？)．しかもたぶんまだ開発中だ．なので、とくに積分機能を試したい等でなければ、これはスキップしておくのも手だろう．
さて、 あとは，おおよそ以下のような感じだ．
変数の宣言
@variables というマクロを使う．
1# 一つでも，複数でも OK.
2@variables x y t
式の定義
1# SymPy と変わらないね．
2eq = x^2 + 3x + 2
$\displaystyle x^2 + 3x + 2$
因数分解
SymbolicUtils などの(サブ)パッケージを使うことを想定しているようだが、まだあまりうまく動かないようだ(エラーにはならないが、結果が望ましくない)．
展開
sympify という，式を整理するコマンドがあるのでそれに展開オプションを付ければ良い．
1simplify( eq * (x^2-1), expand = true )
$\displaystyle x^4 + 3x^3 + x^2 - 3x -2$
代入
substitute コマンドを使う．
1substitute( eq, Dict([x => x^2]))
$\displaystyle x^4 + 3x^2 + 2$
Symbolics package は Julia で書かれているので，次のような線形代数の例でも普通に代入が出来る．
1using LinearAlgebra
2
3A = [ 1 2 3
44 5 6
57 8 t ]
としてからの次のような代入は
1substitute(A, Dict([t => 15]))
無事に
3×3 Matrix{Int64}:
1 2 3
4 5 6
7 8 15
という結果になる．
微分
作用素表現などを用いた，かなり高度なことが出来るようだ． 下記にあるのはいちばん簡単な使い方かな．
1Symbolics.derivative( eq, x )
$\displaystyle 2x + 3$
積分
SymbolicIntegration という積分計算用(サブ)パッケージがあり、これが使える… が、まだ開発途上の様子で、まだ安心して使えるレベルではなさそう．
SymPy を呼び出す方法や SymbolicNumericIntegraton package と組み合わせて必要な積分の結果の式を得る方法もあるようだ．
1integrate( eq, x )
$\displaystyle 2x + \frac{3}{2}x^2 + \frac{1}{3}x^3$
Taylor展開
taylor というコマンドがある．
1taylor( sin(x), x, 0:5)
$\displaystyle x - \frac{1}{6}x^3 + \frac{1}{120}x^5$
級数
series というコマンドで級数を定義できる．
1ss(num) = series( [(1/n)^2 for n in 1:num] , 1 )
ちょっとまわりくどいが，級数の機能を使って $\displaystyle \sum_{n=1}^{\mbox{num}} \left( \frac{1}{n^2} \right)$ を定義している． 次のように使える．
1ss(100)
1.6349839001848923
… 全体にみて，Symbolics package は旧来の数式処理と異なる開発方針のようだし、なにより開発途上だ． 毎年、1年経つとだいぶ機能が増えている．そういう意味で、未実装ないしは不得意，という部分があちこちにみられるため、しばらくは SymPy と両方使い分けていくと良さそうだ．
レポート 15
以下の課題について能う限り賢明な調査と考察を行い，
学籍番号-氏名-15.pdf
というファイルとしてレポートを作成し、
webフォーム
から教官宛に提出しよう.
なお，レポートを「紙媒体」で教官に直接手渡す提出もありだが、物質によるレポート提出は破損や紛失の可能性があるのであまりお勧めしないぞ．
注意
近年はセキュリティ上の懸念から，実行形式のプログラムなどをメールに添付するとそのメールそのものの受信を受信側サーバが拒絶したりする．
そういうことを避けるため，レポートをメールで提出するときは添付ファイルにそういった懸念のあるファイルが無いようにしよう．
まあ要するに，レポートは pdf ファイルで提出するのが良い と思っておこう．
一応、レポートやメールには，自分の 所属・学年・学籍番号・氏名 を必ず書いておこう．なにかと安全だ．
課題
- この授業のこれまでで手計算が必要だった箇所を探し出し，2例ほど上の SymPy を使う形で数式処理してみよう．
そして，その結果と自分の手計算の結果を比較しよう．
- 自分の持っている微分積分学，線形代数学の教科書を開き，手計算でこれまで計算してきた計算を 2例ほど上の SymPy を使う形で数式処理してみよう．
そして，その結果と自分の手計算の結果を比較しよう．