Python入門トップページ

SymPy による数式処理を使ってみよう

目次

  1. 数式処理システムとは
  2. モジュールのインポート
  3. 変数の定義
  4. 数式の定義
  5. 変数に値を代入して数式を評価
  6. 数式の展開
  7. 因数分解
  8. 数式の簡略化
  9. 素因数分解と素数判定
  10. 方程式を解く
    1. 方程式
    2. 連立方程式
  11. 微分
  12. 積分
    1. 不定積分
    2. 定積分
  13. 指数関数
  14. 確率分布
    1. 正規分布
  15. 任意の関数のグラフを描く
  16. 円や楕円を描く
  17. 参考サイト

目次に戻る

数式処理システムとは

数式処理システムとは \( 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
\( 3x \) 
数式の定義
expr = a*x**2 + b*x + c
expr
\( ax^2 + bx + c \)

上の例のように，変数名だけを入力して評価すると数式が表示されるが，print() を使うと次のように表示されます． 

数式の定義
expr = 3*x
print(expr)
3*x 
数式の定義
expr = a*x**2 + b*x + c
print(expr)
a*x**2 + b*x + c

目次に戻る

変数に値を代入して数式を評価

定義した数式に値を代入して評価するには，subs() を使います．再度，数式を定義してみよう． 

数式の定義
expr = 3*x
expr
\( 3x \)

定義した数式に \(x = 5\) を代入して評価してみよう． 

値を代入して評価
expr.subs(x,5)
\( 15 \)

数式に \(x = 5\) を代入しても，数式自身が上書きされているわけではない． 

数式の定義を確認
expr
\( 3x \)

別の関数を定義します． 

数式の定義
expr = a*x**2 + b*x + c
expr
\( ax^2 + bx + c \)

\(a\)，\(b\)，\(c\) の3つに値を代入して数式を評価します． 

代入して評価
expr.subs([(a,2),(b,3),(c,-4)])
\( 2x^2 + 3x - 4 \)

さらに，\(x\) にも代入して数式を評価します． 

代入して評価
expr.subs([(a,2),(b,3),(c,-4),(x,2)])
\( 10 \)

目次に戻る

数式の展開

数式の展開は sympy.expand() を使います．まず，数式を定義します． 

数式の定義
expr = (a+b)**2
expr
\( (a + b)^2 \)

上の数式を展開します．

数式の展開
sympy.expand(expr)
\(a^2 + 2ab + b^2\)

また別の数式を定義します． 

数式の定義
expr = (a+b)**3
expr
\( (a + b)^3 \)

上の式を展開します． 

数式の展開
sympy.expand(expr)
\(a^3 + 3a^2b + 3ab^2 + b^3\)

さらに数式を定義して展開します． 

数式の定義
expr = (a+b+c)**2 + (b+c-a)**2 + (c+a-b)**2 + (a+b-c)**2
expr
\( (-a + b + c)^2 + (a - b + c)^2 + (a + b -c)^2 + (a + b + c)^2 \) 
数式の展開
sympy.expand(expr)
\(4a^2 + 4b^2 + 4c^2\)

目次に戻る

因数分解

式の因数分解は sympy.factor() を使います．まず，数式を定義します． 

数式の定義
expr = a**2 + 2*a*b + b**2
expr
\( a^2 + 2ab + b^2 \)

上の数式を因数分解します． 

因数分解
sympy.factor(expr)
\(\left( a+b \right)^2\)

別の数式を定義します． 

数式の定義
expr = a**3 + 3*a**2*b + 3*a*b**2 + b**3
expr
\( a^3 + 3a^2b + 3ab^2 + b^3 \)

上の数式を因数分解します． 

因数分解
sympy.factor(expr)
\(\left( a+b \right)^3\)

さらに別の数式を定義します． 

数式の定義
expr = 3*a**2 + 7*a*b + 2*b**2 - 5*a - 5*b + 2
expr
\( 3a^2 + 7ab - 5a + 2b^2 - 5b + 2 \)

上の数式を因数分解します． 

因数分解
sympy.factor(expr)
\( \left( a+2b-1 \right)\left( 3b+b-2 \right) \)

目次に戻る

数式の簡単化

数式を簡単な形式に変換するためには sympy.simplify() を利用すれば良いでしょう．まず，数式を定義します． 

数式の定義
expr = (a + b + 2*c)**3 - (b + 2*c - a)**3 - (2*c + a - b)**3 - (a + b - 2*c)**3
expr
\( -(-a+b+2c)^3 - (a-b+2c)^3 - (a+b-2c)^3 + (a+b+2c)^3 \)

上の数式を簡単化します． 

式の簡単化
sympy.simplify(expr)
\( 48abc \)

目次に戻る

素因数分解と素数判定

素因数分解も可能です．\(4725 = 3^3 \times 5^2 \times 7^1 \) を素因数分解してみよう． 

素因数分解
sympy.factorint(4725)
{3: 3, 5: 2, 7: 1}

素数であるかどうかのチェックをしてみよう．まず，9941 は素数です． 

素数判定
sympy.isprime(9941)
True

9943 は素数ではなく合成数であるので，素因数分解もしてみよう． 

素数判定
sympy.isprime(9943)
False

素因数分解の結果 \(9943 = 61 \times 163\) であることが分かった． 

素因数分解
sympy.factorint(9943)
{61: 1, 163: 1}

12345678910987654321 も素数です． 

素数判定
sympy.isprime(12345678910987654321)
True

非常に桁の大きな数についても素数判定をしてみよう．よく知られた \( 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)
True

面白いことに 31 から始まり 331 や 3,331 などはしばらく素数です．なお，桁の大きな整数はアンダースコアで区切ることができます． 

素数判定
sympy.isprime(31)
True 
素数判定
sympy.isprime(331)
True 
素数判定
sympy.isprime(3_331)
True 
素数判定
sympy.isprime(33_331)
True 
素数判定
sympy.isprime(333_331)
True 
素数判定
sympy.isprime(3_333_331)
True 
素数判定
sympy.isprime(33_333_331)
True

しかし，333,333,331 は素数ではありません． 

素数判定
sympy.isprime(333_333_331)
False

333,333,331 を素因数分解すると \(333,333,331 = 17 \times 19,607,843\) であることが分かりました． 

素因数分解
sympy.factorint(333_333_331)
{17: 1, 19607843: 1}

その後はしばらく合成数が続きますが，333,333,333,333,333,331 はまた素数です． 

素数判定
sympy.isprime(333_333_333_333_333_331)
True

目次に戻る

方程式を解く

方程式

例えば，\(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)
[1/4, 3/2]

次に，\(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)
[4/3 - sqrt(5)/3, sqrt(5)/3 + 4/3]

さらに，\(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)
[7/10 - sqrt(29)/10, sqrt(29)/10 + 7/10]

解の公式をそのものの表示も可能です． 

解の公式
sympy.solve(a*x**2 + b*x + c, x)
[(-b + sqrt(-4*a*c + b**2))/(2*a), -(b + sqrt(-4*a*c + b**2))/(2*a)]

連立方程式

連立方程式 \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])
{x: -1, y: 1}

次の連立方程式 \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])
[(-2, -1), (8/5, 4/5), (-sqrt(2), -sqrt(2)), (sqrt(2), sqrt(2))]

目次に戻る

微分

SymPy では関数の微分を代数的に求めることも可能です．例えば，\( f(x) = ax^2 - bx+c \) の導関数は次のように求めることができます． 

微分
sympy.diff(a*x**2 - b*x - c, x)
\( 2ax - b \)

目次に戻る

積分

不定積分

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)
\( \displaystyle 2x^4 + \frac{x^3}{3} - 2x^2 + 2x \)

なお，上の不定積分の結果には積分定数 \( C \) が表示されていないことに注意が必要です．

次に \begin{eqnarray*} \int x^n dx \end{eqnarray*} を求めよう． 

不定積分
n = sympy.Symbol('n')
sympy.integrate(x**n, x)
\begin{equation} \left\{ \begin{array}{ll} \displaystyle \frac{x^{n+1}}{n+1} & \mbox{for } n \neq -1 \\ \log (x) & \mbox{otherwise} \end{array} \right. \end{equation}

定積分

定積分 \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))
\( -6 \)

もちろん定数を含むような定積分 \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))
\( \displaystyle -a^3 - a^2 + a + b^3 + b^2 - b \)

目次に戻る

指数関数

SymPy での自然対数の底（ネイピア数）は sympy.E，指数関数は sympy.exp() です． 

自然対数の底
sympy.E
\( e \) 
指数関数
sympy.exp(x)
\( e^x \)

\( e^x \) は \( x \) で微分しても \( e^x \) です． 

指数関数
expr = sympy.exp(x)
expr
\( e^x \) 
指数関数の微分
sympy.diff(expr, x)
\( e^x \)

\( e^{ax} \) を \( x \) で微分すると \( ae^{ax} \) となります． 

指数関数
expr = sympy.exp(a * x)
expr
\( e^{ax} \) 
指数関数の微分
sympy.diff(expr, x)
\( ae^{ax} \)

指数関数の積分も SymPy で確認しよう．

\( \displaystyle \int -ae^{-ax}dx = e^{-ax} \) 
指数関数
expr = - a * sympy.exp(- a * x)
expr
\( -ae^{-ax} \) 
指数関数の積分
sympy.integrate(expr, x)
\( e^{-ax} \)

目次に戻る

確率分布

正規分布

正規分布 \(\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
\( e^{-\frac{(-\mu + x)^2}{2\sigma^2}} \)

上で定義した関数を \(-\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)
\( \sqrt{2}\sqrt{\pi}\sigma \)

つまり上の結果から，正規分布の密度関数にある \(\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
\( \frac{\sqrt{2} e^{-\frac{(-\mu + x)^2}{2\sigma^2}}}{2 \sqrt{\pi} \sigma} \)

期待値の定義通り積分すると \(\mu\) になることが確認できました． 

正規分布の平均
expr = sympy.integrate(x * density, (x,- sympy.oo,sympy.oo))
sympy.simplify(expr)
\( \mu \)

分散の定義通り計算するとやはり \(\sigma^2\) になることが確認できました． 

正規分布の分散
expr = sympy.integrate(x ** 2 * density, (x,- sympy.oo,sympy.oo)) - mu ** 2
sympy.simplify(expr)
\( \sigma^2 \)

目次に戻る

任意の関数のグラフを描く

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')
sympy-plot-1

もう一つの関数 \( y = x^2 - 8 x + 4 \) のグラフも描画します． 


p2 = sympy.plot(x**2 - 8*x + 4, (x, 0, 8), xlabel='x', ylabel='y')
sympy-plot-2

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()
sympy-plot-3

目次に戻る

円や楕円を描く

半径 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));
circle2022-1

また，媒介変数表示にすることなく円の方程式から直接描画することも可能です．このためには sympy.Eq() を使って方程式を定義します． 


eqn = sympy.Eq(x**2 + y**2, 1)
eqn
\( x^2 + y^2 = 1 \)

次に sympy.plot_implicit() を使って描画します． 


plt.rcParams['figure.figsize'] = (6, 6)
sympy.plot_implicit(eqn, (x, -1, 1), (y, -1, 1));
circle2022-2

楕円の方程式を定義して描画します． 


eqn = sympy.Eq(0.5* x**2 + 2 * y**2, 2)
eqn
\( 0.5x^2 + 2y^2 = 2 \) 

plt.rcParams['figure.figsize'] = (6, 6)
sympy.plot_implicit(eqn, (x, -2, 2), (y, -2, 2));
circle2022-3

次は \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
\( x^2 + \left( y - \sqrt{|x|} \right)^2 = 1 \) 

plt.rcParams['figure.figsize'] = (6, 6)
sympy.plot_implicit(eqn, (x, -1, 1), (y, -1, 1.7));
circle2022-4

ハート型を描くことができました．

さらに \begin{equation} \left( x^2 + y^2 \right)^2 - \left( 3x^2 - y^2 \right)y = 0 \end{equation} なる方程式を \( -1 \leq x \leq 1\) の範囲で描いてどのようなグラフになるか確認しよう．

目次に戻る

参考サイト

SymPy のオリジナルサイトはこちら

目次に戻る