SymPy による数式処理を使ってみよう
目次
- 数式処理システムとは
- モジュールのインポート
- 変数の定義
- 数式の定義
- 変数に値を代入して数式を評価
- 数式の展開
- 因数分解
- 数式の簡略化
- 素因数分解と素数判定
- 方程式を解く
- 微分
- 積分
- 指数関数
- 確率分布
- 任意の関数のグラフを描く
- 円や楕円を描く
- 参考サイト
数式処理システムとは
数式処理システムとは \( x \) や \( y \) などの変数を含んだ数式を可能な限り代数的に扱い処理をするシステムである．代表的な数式処理システムは Mathematica や Maple などである．これらの数式処理システムは非常に高機能であるが，高価でもある．
Anaconda を使える環境であれば，SymPy モジュールをインポートすることで，Mathematica や Maple と同じような数式処理システムを無料で使えるようになる．ここでは，SymPy を使って数式処理を体験してみよう．
モジュールのインポート
まずは SymPy モジュールをインポートする．
モジュールのインポートimport sympy
変数の定義
使用したい変数はあらかじめ定義しなければならない．これから使う予定の変数をここで定義しておく．
変数の定義a = sympy.Symbol('a')
b = sympy.Symbol('b')
c = sympy.Symbol('c')
x = sympy.Symbol('x')
y = sympy.Symbol('y')
数式の定義
変数の定義を行っておれば，数式の定義は通常の算術演算子などを使って可能である．
数式の定義expr = 3*x
expr
数式の定義expr = a*x**2 + b*x + c
expr
上の例のように，変数名だけを入力して評価すると数式が表示されるが，print() を使うと次のように表示される．
数式の定義expr = 3*x
print(expr)
数式の定義expr = a*x**2 + b*x + c
print(expr)
変数に値を代入して数式を評価
定義した数式に値を代入して評価するには，subs() を使う．再度，数式を定義してみよう．
数式の定義expr = 3*x
expr
定義した数式に \(x = 5\) を代入して評価してみよう．
値を代入して評価expr.subs(x,5)
数式に \(x = 5\) を代入しても，数式自身が上書きされているわけではない．
数式の定義を確認expr
別の関数を定義する．
数式の定義expr = a*x**2 + b*x + c
expr
\(a\)，\(b\)，\(c\) の3つに値を代入して数式を評価する．
代入して評価expr.subs([(a,2),(b,3),(c,-4)])
さらに，\(x\) にも代入して数式を評価する．
代入して評価expr.subs([(a,2),(b,3),(c,-4),(x,2)])
数式の展開
数式の展開は sympy.expand() を使う．まず，数式を定義する．
数式の定義expr = (a+b)**2
expr
上の数式を展開する．
数式の展開sympy.expand(expr)
また別の数式を定義する．
数式の定義expr = (a+b)**3
expr
上の式を展開する．
数式の展開sympy.expand(expr)
さらに数式を定義して展開する．
数式の定義expr = (a+b+c)**2 + (b+c-a)**2 + (c+a-b)**2 + (a+b-c)**2
expr
数式の展開sympy.expand(expr)
因数分解
式の因数分解は sympy.factor() を使う．まず，数式を定義する．
数式の定義expr = a**2 + 2*a*b + b**2
expr
上の数式を因数分解する．
因数分解sympy.factor(expr)
別の数式を定義する．
数式の定義expr = a**3 + 3*a**2*b + 3*a*b**2 + b**3
expr
上の数式を因数分解する．
因数分解sympy.factor(expr)
さらに別の数式を定義する．
数式の定義expr = 3*a**2 + 7*a*b + 2*b**2 - 5*a - 5*b + 2
expr
上の数式を因数分解する．
因数分解sympy.factor(expr)
数式の簡単化
数式を簡単な形式に変換するためには sympy.simplify() を利用すれば良い．まず，数式を定義する．
数式の定義expr = (a + b + 2*c)**3 - (b + 2*c - a)**3 - (2*c + a - b)**3 - (a + b - 2*c)**3
expr
上の数式を簡単化する．
式の簡単化sympy.simplify(expr)
素因数分解と素数判定
素因数分解も可能である．\(4725 = 3^3 \times 5^2 \times 7^1 \) を素因数分解してみよう．
素因数分解sympy.factorint(4725)
素数であるかどうかのチェックをしてみよう．まず，9941 は素数である．
素数判定sympy.isprime(9941)
9943 は素数ではなく合成数であるので，素因数分解もしてみよう．
素数判定sympy.isprime(9943)
素因数分解の結果 \(9943 = 61 \times 163\) であることが分かった．
素因数分解sympy.factorint(9943)
非常に桁の大きな数についても素数判定をしてみよう．よく知られた \( 2^{n} - 1 \) で得られる「メルセンヌ素数」の19番目は \( 2^{4253} - 1 \) であり，これは 1281 桁の素数である（50桁ごとに改行して表示している）．
19番目のメルセンヌ素数2**4253 - 1
19079700752443907380746804296952917366935699474994 01773947418826735289797870050537063680498355149002 44303495954950709725762186311224148828811920216904 54220696074466616936422119528953843684539025016866 39328388051920551371543909126665275330073092926875 39092257043362517857366624699975402375462954490293 25923330313733064353155653973992192620143860643902 00751747230290568382725050515719675946083500634044 95977660656269020823960825567012344189908927956646 01199805798854863010763738099351982658238978188813 57054086530452196558017580812511640805546090574680 28203308718724654081055323215860189611391296030471 10844314674567196776630892585854727150731156376517 10083182486471100976148903135628565417841548817431 46033909602737947385055355960331855614540900081456 37865906837031726769698000118775099549109035010841 70509179915621679722810701613059725180448720483313 06383715094854938415738549894606070722584737978176 68642213435452698944302835364403718737538539783825 95118331664161343236956603676768977222879187734209 68982326089026150031515424165462111337527431154890 66632737492144627683356451977679763387550354866509 39145564820314822488831270237770396677079765598573 33357013727342079099064400455741830654320379350833 23624581934882406478358569292488102197833297494990 6122664421376034687815350484991
このようなメルセンヌ素数の判定は短時間で可能です．
素数判定sympy.isprime(2**4253 - 1)
面白いことに 31 から始まり 331 や 3,331 などはしばらく素数です．なお，桁の大きな整数はアンダースコアで区切ることができます．
素数判定sympy.isprime(31)
素数判定sympy.isprime(331)
素数判定sympy.isprime(3_331)
素数判定sympy.isprime(33_331)
素数判定sympy.isprime(333_331)
素数判定sympy.isprime(3_333_331)
素数判定sympy.isprime(33_333_331)
しかし，333,333,331 は素数ではありません．
素数判定sympy.isprime(333_333_331)
333,333,331 を素因数分解すると \(333,333,331 = 17 \times 19,607,843\) であることが分かりました．
素因数分解sympy.factorint(333_333_331)
その後はしばらく合成数が続きますが，333,333,333,333,333,331 はまた素数です．
素数判定sympy.isprime(333_333_333_333_333_331)
方程式を解く
方程式
例えば，\(x\) の方程式 \(8x^2-14x+3 =0\) を解くときは
\begin{eqnarray}
& & 8x^2-14x+3 = 0 \\
& & (4x-1)(2x-3) = 0
\end{eqnarray}
と因数分解して，\(\displaystyle x = \frac{1}{4},~\frac{3}{2}\) という解を得る．これを SymPy で解いてみよう．まだ変数 x が定義されていない場合は，後の連立方程式で使う y と一緒に次のコマンドで定義をしておく．
変数の定義x = sympy.Symbol('x')
y = sympy.Symbol('y')
SymPy で方程式を解くためには，sympy.solve() を利用すれば良い．
方程式sympy.solve(8*x**2-14*x+3, x)
次に，\(x\) の方程式 \((3x-4)^2=5\) は
\begin{eqnarray}
(3x-4)^2 &=& 5
\\
3x -4 &=& \pm \sqrt{5}
\\
3x &=& 4 \pm \sqrt{5}
\\
x &=& \frac{4 \pm \sqrt{5}}{3}
\end{eqnarray}
である．これを SymPy で解くためには
\begin{eqnarray}
(3x-4)^2 - 5 &=& 0
\end{eqnarray}
のように 「左辺 = 0」の形に変形して sympy.solve() に与えるとよい．
方程式sympy.solve((3*x - 4)**2 - 5, x)
さらに，\(x\) の方程式 \( 5x^2 - 7x + 1 = 0 \) は解の公式を利用して \begin{eqnarray} x &=& \frac{7\pm \sqrt{29}}{10} \end{eqnarray} である．これも sympy.solve で解くことができる．
方程式sympy.solve(5*x**2 - 7*x + 1, x)
解の公式をそのものの表示も可能である．
解の公式sympy.solve(a*x**2 + b*x + c, x)
連立方程式
連立方程式 \begin{eqnarray} \left\{ \begin{array}{l} x + 3y = 2\\ -2x + 3y = 5 \end{array} \right. \end{eqnarray} も sympy.solve で解くことができる．
連立方程式sympy.solve([x + 3*y - 2, -2*x + 3*y - 5], [x,y])
次の連立方程式 \begin{eqnarray} \left\{ \begin{array}{l} x^2 - 3xy + 2y^2 = 0\\ x^2 + y^2 + x - y = 4 \end{array} \right. \end{eqnarray} も sympy.solve で解くことができる．
連立方程式sympy.solve([x**2 - 3*x*y + 2*y**2, x**2 + y**2 + x - y - 4], [x,y])
微分
SymPy では関数の微分を代数的に求めることも可能である．例えば，\( f(x) = ax^2 - bx+c \) の導関数は次のように求めることができる．
微分sympy.diff(a*x**2 - b*x - c, x)
積分
不定積分
SymPy で積分を計算するには，sympy.integrate を利用する．例えば，
\begin{eqnarray*}
\int \left( 8x^3 + x^2 - 4x + 2 \right)dx
&=& 2x^4 + \frac{x^3}{3} - 2x^2 + 2x + C
\end{eqnarray*}
を SymPy で求めてみよう．
不定積分sympy.integrate(8*x**3 + x**2 - 4*x + 2, x)
なお，上の不定積分の結果には積分定数 \( C \) が表示されていないことに注意が必要である．
次に \begin{eqnarray*} \int x^n dx \end{eqnarray*} を求めよう．
不定積分n = sympy.Symbol('n')
sympy.integrate(x**n, x)
定積分
定積分 \begin{eqnarray*} \int_{0}^{2} \left( x^3 - 3x^2 - 1 \right)dx &=& \left[ \frac{x^4}{4} - x^3 - x \right]^{2}_{0} \\ &=& (4-8-2) - 0 \\ &=& -6 \end{eqnarray*} を SymPy で求めてみよう．
定積分sympy.integrate(x**3 - 3*x**2 - 1, (x,0,2))
もちろん定数を含むような定積分 \begin{eqnarray*} \int_{a}^{b} \left( 3x^2 + 2x - 1 \right) dx \end{eqnarray*} も SymPy で求めることができる．
定積分sympy.integrate(3*x**2 + 2*x - 1, (x,a,b))
指数関数
SymPy での自然対数の底（ネイピア数）は sympy.E，指数関数は sympy.exp() である．
自然対数の底sympy.E
指数関数sympy.exp(x)
\( e^x \) は \( x \) で微分しても \( e^x \) である．
指数関数expr = sympy.exp(x)
expr
指数関数の微分sympy.diff(expr, x)
\( e^{ax} \) を \( x \) で微分すると \( ae^{ax} \) となる．
指数関数expr = sympy.exp(a * x)
expr
指数関数の微分sympy.diff(expr, x)
指数関数の積分も SymPy で確認しよう．
指数関数expr = - a * sympy.exp(- a * x)
expr
指数関数の積分sympy.integrate(expr, x)
確率分布
正規分布
正規分布 \(\mathcal{N}\left(\mu, \sigma^2\right)\) は統計学の基本的な分布である．なお，平均が \(\mu\)，分散が \(\sigma^2\) である．正規分布の密度関数
\begin{eqnarray*}
f(x)
= \frac{1}{\sqrt{2\pi} \sigma}
\exp \left( - \frac{(x-\mu)^2}{2\sigma^2}\right)
\end{eqnarray*}
について，定数項を除いた \(\exp(\cdot)\) を \(-\infty\) から \(\infty\) の範囲で積分すると，\(\sqrt{2\pi}\sigma\) になることが知られている．このことを確認してみよう．まずは，変数 \(x\)，\(\mu\)，\(\sigma\) を定義する．ただし，\(\sigma > 0\) の条件を positive=True によって指定しておく．また，TeXの表記を使うとギリシア文字も利用できることも注意する．
変数の定義x = sympy.Symbol('x')
mu = sympy.Symbol('\mu')
sigma = sympy.Symbol('\sigma', positive=True)
次に，積分したい関数を定義する．
関数の定義fn = sympy.exp(-(x-mu)**2 / (2 * sigma**2))
fn
上で定義した関数を \(-\infty\) から \(\infty\) の範囲で積分し，
\begin{eqnarray*}
\int_{-\infty}^{\infty}
\exp \left( - \frac{(x-\mu)^2}{2\sigma^2}\right) dx
= \sqrt{2\pi}\sigma
\end{eqnarray*}
となることを確認してみよう．このとき，\(\infty\) は SymPy では sympy.oo （アルファベット小文字の「オー」を2個）で入力することができる．
積分するexpr = sympy.integrate(fn, (x,- sympy.oo,sympy.oo))
sympy.simplify(expr)
つまり上の結果から，正規分布の密度関数にある \(\frac{1}{\sqrt{2\pi}\sigma}\) という項は積分した結果が1になるようにするための定数項であることがわかります．
さらに，正規分布の平均（期待値）が \(\mu\)，分散が \(\sigma^2\) になることも確認しておこう．期待値の定義は \begin{equation} E[X] = \int_{-\infty}^{\infty} x f(x) dx \end{equation} であり，分散の定義は \begin{equation} {\rm{Var}}[X] = \int_{-\infty}^{\infty} x^2 f(x) dx - \mu^2 \end{equation} である．
まず，密度関数を定義する．
密度関数の定義density = 1 / (sympy.sqrt(2 * sympy.pi) * sigma) * sympy.exp(-(x-mu)**2 / (2 * sigma**2))
density
期待値の定義通り積分すると \(\mu\) になることが確認できました．
正規分布の平均expr = sympy.integrate(x * density, (x,- sympy.oo,sympy.oo))
sympy.simplify(expr)
分散の定義通り計算するとやはり \(\sigma^2\) になることが確認できました．
正規分布の分散expr = sympy.integrate(x ** 2 * density, (x,- sympy.oo,sympy.oo)) - mu ** 2
sympy.simplify(expr)
任意の関数のグラフを描く
Matplotlib では x と y に対応する値の配列を準備することで任意の関数のグラフを作成することができました．しかし，Sympy では数学関数をそのままグラフとして描くことが可能です．ここでは \( y = - x^2 + 8 x - 4 \) と \( y = x^2 - 8 x + 4 \) のグラフ描くことにする．
グラフを描くために予め matplotlib ライブラリをインポートしておきます．
import matplotlib.pyplot as plt
# 以下は高解像度ディスプレイ用の設定
from IPython.display import set_matplotlib_formats
# from matplotlib_inline.backend_inline import set_matplotlib_formats # バージョンによってはこちらを有効に
set_matplotlib_formats('retina')
利用する変数 \(x\) を定義しておきます．
x = sympy.Symbol('x')
最初の関数 \( y = - x^2 + 8 x - 4 \) を sympy.plot() を使って描画します．なお，第2引数の (x, 0, 8) はグラフを描画する x の範囲をタプル形式で指定しています．
p1 = sympy.plot(-x**2 + 8*x - 4, (x, 0, 8), xlabel='x', ylabel='y')
もう一つの関数 \( y = x^2 - 8 x + 4 \) のグラフも描画します．
p2 = sympy.plot(x**2 - 8*x + 4, (x, 0, 8), xlabel='x', ylabel='y')
2つのグラフを重ねるには extend() メソッドを利用します．このとき，それぞれのグラフ作成時に show=False を指定することで，個別のグラフの描画をしないようにすることもできます．
p1 = sympy.plot(-x**2 + 8*x - 4, (x, 0, 8), xlabel='x', ylabel='y', show=False)
p2 = sympy.plot(x**2 - 8*x + 4, (x, 0, 8), xlabel='x', ylabel='y', show=False)
p1.extend(p2)
p1.show()
円や楕円を描く
半径 1 の円の方程式は次の通りです． \begin{equation} x^2 + y^2 = 1 \end{equation} この円を描くには次のような媒介変数表示にします． \begin{equation} x = \cos\theta, \quad y = \sin \theta \end{equation} グラフを描くために予め matplotlib ライブラリをインポートしておきます．
import matplotlib.pyplot as plt
# 以下は高解像度ディスプレイ用の設定
from IPython.display import set_matplotlib_formats
# from matplotlib_inline.backend_inline import set_matplotlib_formats # バージョンによってはこちらを有効に
set_matplotlib_formats('retina')
媒介変数として利用する変数 \(t\) を定義しておきます．また，\(x\), \(y\) の定義がまだであればこれらも定義しておきます．
t = sympy.Symbol('t')
x = sympy.Symbol('x')
y = sympy.Symbol('y')
媒介変数表示でグラフを描くには sympy.plot_parametric() を利用します．
plt.rcParams['figure.figsize'] = (6, 6)
sympy.plot_parametric(sympy.cos(t), sympy.sin(t), (t, 0, 2*sympy.pi));
また，媒介変数表示にすることなく円の方程式から直接描画することも可能です．このためには sympy.Eq() を使って方程式を定義します．
eqn = sympy.Eq(x**2 + y**2, 1)
eqn
次に sympy.plot_implicit() を使って描画します．
plt.rcParams['figure.figsize'] = (6, 6)
sympy.plot_implicit(eqn, (x, -1, 1), (y, -1, 1));
楕円の方程式を定義して描画します．
eqn = sympy.Eq(0.5* x**2 + 2 * y**2, 2)
eqn
plt.rcParams['figure.figsize'] = (6, 6)
sympy.plot_implicit(eqn, (x, -2, 2), (y, -2, 2));
次は \begin{equation} x^2 + \left( y - \sqrt{|x|} \right)^2 = 1 \end{equation} なる方程式を描いてみよう．
eqn = sympy.Eq(x**2 + (y - sympy.sqrt(sympy.Abs(x))) ** 2, 1)
eqn
plt.rcParams['figure.figsize'] = (6, 6)
sympy.plot_implicit(eqn, (x, -1, 1), (y, -1, 1.7));
ハート型を描くことができました．
さらに \begin{equation} \left( x^2 + y^2 \right)^2 - \left( 3x^2 - y^2 \right)y = 0 \end{equation} なる方程式を \( -1 \leq x \leq 1\) の範囲で描いてどのようなグラフになるか確認しよう．