A1. SymPy in Julia: 数式処理

Photo by Jeswin Thomas on Unsplash

数式処理 超入門

Julia ではいくつかの数式処理システムを利用できる． その中でもっともよく使われているものが SymPy package だろう．

この package の使い方だが，幸い，package 自身にドキュメントの素になる資料がついてくる． この資料を julia でコンパイルして html にすると web browser で読めるようになるのだが，面倒だろうから， 読めるようにした

package に付属する SymPy document

を用意した． この資料は入門用資料として結構よくできているので，本来はこれを一通り読むのが望ましいなあ．できれば読んでおこう．

はじめに

そして，ここではなるべく簡単に SymPy の使い方を紹介しよう． まず，インストールしていない人は SymPy package をインストールしよう．ちょっと時間がかかるぞ. なお，既にインストールしてある人はもちろんこれをしなくてよい．

1
2

using Pkg
Pkg.add("SymPy")

インストールが終わったら，次にこのパッケージの利用宣言をしよう．

1

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

そして，SymPy のポイントを先に説明しよう． そのポイントとは，

SymPy では 「数式処理の対象となる変数を先に陽に指定する」 という仕組み

である，とでもいうものである． 言い換えると，「指定すると特別な変数」，「指定していない変数は，いつもどおり Julia の普通の変数」であるのだ． そのあたりを例で示そう．

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

1

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

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

1

p(x)

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

例えば

1

p(3)

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

1
2

using Plots
plot(p)

とすると

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

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

1
2

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

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

1

eq = p(t)

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

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

まあもう少し書くと，ここでは，「関数は計算手順を書いたもの」であるが「数式は(特殊な)文字列」だ，と思っておけば良い．

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

1
2

# 因数分解
factor( eq )

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

1
2

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

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

1
2

# 微分
diff( eq, t )

$2t+3$

1
2

# 不定積分
integrate( eq, t )

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

1
2

# 定積分
integrate( 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つ書く)は無限大を表す記号である．

1
2
3
4

# 変数が複数の関数も問題ない．
@vars u v

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

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

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

数式を関数に戻す/変換するには？

上の例のように，関数から数式を作るのは単に「特別な変数を代入」するだけでよいので簡単だった．

では逆はどうだろう？ せっかく数式処理して得た結果の数式を関数に直して計算に使いたい，というのは自然な要求だ． それは簡単で，たとえば上の数式 eq2 を 変数 t をもらって t の4次式に基づいた計算値を返す関数に変換するには，

1

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

とすれば良い¹． これで数式 eq2 を変換した関数 p4(t) を定義したことになる(ただしこの p4(t) の中の t は SymPy と関係なく，単なる普通の変数扱い)． これで得られた p4(t) は確かに上で求めた 4次多項式の値を計算する関数になっている．

複数の入力変数があるときでも同様で，たとえば上の数式 eq3 に対しては

1

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

とすれば，関数 q(u,v) を定義したことになる．

ちなみに，関数に変換するたびに毎回 eval(Meta.parse... と書くのは大変だという人は，using SymPy の直後にいつも

1

toFunc(str) = eval(Meta.parse(str))

という感じの定義をつけておけば良い． こうしておけば，上の p4(t) を定義するにも

1

toFunc("p4(t) = $eq2")

と書けばよい．少しだけ楽ができるし，構文を覚えておかなくても良いね．

さて，変換した結果は通常の関数だから，たとえば p4 を使って t = 3 の時のこの4次多項式の値は

1

p4(3)

160

と普通に計算できる．

これが普通の関数であることは，次のように通常の関数と比較することでも確認できる． まず，この4次多項式を普通の関数で作る．

1

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

そして，同じ入力値に対してのそれぞれの動作を @code_lowered マクロで比較してみれば良い．

1

@code_lowered p4(3)

CodeInfo(
1 ─ %1  = Core.apply_type(Base.Val, 4)
│   %2  = (%1)()
│   %3  = Base.literal_pow(Main.:^, t, %2)
│   %4  = Core.apply_type(Base.Val, 3)
│   %5  = (%4)()
│   %6  = Base.literal_pow(Main.:^, t, %5)
│   %7  = 3 * %6
│   %8  = Core.apply_type(Base.Val, 2)
│   %9  = (%8)()
│   %10 = Base.literal_pow(Main.:^, t, %9)
│   %11 = %3 + %7 + %10
│   %12 = 3 * t
│   %13 = %11 - %12
│   %14 = %13 - 2
└──       return %14
)

1

@code_lowered p5(3)

CodeInfo(
1 ─ %1  = Core.apply_type(Base.Val, 4)
│   %2  = (%1)()
│   %3  = Base.literal_pow(Main.:^, x, %2)
│   %4  = Core.apply_type(Base.Val, 3)
│   %5  = (%4)()
│   %6  = Base.literal_pow(Main.:^, x, %5)
│   %7  = 3 * %6
│   %8  = Core.apply_type(Base.Val, 2)
│   %9  = (%8)()
│   %10 = Base.literal_pow(Main.:^, x, %9)
│   %11 = %3 + %7 + %10
│   %12 = 3 * x
│   %13 = %11 - %12
│   %14 = %13 - 2
└──       return %14
)

この2つの結果を比べると，入力変数名が異なるほかは全く同じであることがわかるだろう．安心だね．

Plot

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

文法は少しだけいつもと異なるので， SymPy package Plotting を眺めておくと良い． この資料に載っている以上に細かい plot をしたければ，上に書いてあるように，いったん普通の関数に変換してから plot すると良いだろう．