3.2. Sympy : Python での代数計算¶
著者: Fabian Pedregosa
目標
任意精度での数式の評価。
代数表現の代数的な操作の実行。
- 基本的な微積分(極限, 微分
, 積分) を代数表現で行なう。
多項式や超越方程式の求解。
いくつかの微分方程式の求解。
SymPy とは? SymPy は Python の代数計算ライブラリです. コードをシンプルに保ち、拡張しやすいように保ちコードのシンプルに保ちつつ Mathematica や Maple のようなシステムの代替となることを目指しています。SymPy は全て Python で書かれていて外部ライブラリを必要としません.
Sympy のドキュメントとインストール用パッケージは http://www.sympy.org/ にあります
章の内容
3.2.1. SymPy での第一歩¶
3.2.1.1. Sympy を計算機として使う¶
Sympy は3つの数値型を持っています: Real, Rational そして Integer.
Rational クラスは分子と分母の2つの Integer の対として有理数を表現します, つまり Rational(1,2) は 1/2 を表し, Rational(5,2) は 5/2 を表します:
>>> from sympy import *
>>> a = Rational(1,2)
>>> a
1/2
>>> a*2
1
Sympy はバックグラウンドで mpmath を利用します, これによって任意精度数値演算を実行できます. そうすることで, いくつかの特殊な定数 e, pi, oo (無限大) を symbol として扱い, さらに任意精度で評価することができます:
>>> pi**2
2
pi
>>> pi.evalf()
3.14159265358979
>>> (pi + exp(1)).evalf()
5.85987448204884
上でみるように, evalf は浮動小数点数として, 数式を評価します.
oo と呼ばれる数学的な無限大を表すクラスもあります:
>>> oo > 99999
True
>>> oo + 1
oo
3.2.1.2. 練習問題¶
\(\sqrt{2}\) を100桁まで計算しましょう.
\(1/2 + 1/3\) を有理数として計算しましょう.
3.2.1.3. Symbols¶
他の計算代数システムと対照的に, Sympy では symbol として使う変数を明示的に宣言しなければいけません:
>>> from sympy import *
>>> x = Symbol('x')
>>> y = Symbol('y')
そしてこれらを算術操作できます:
>>> x + y + x - y
2*x
>>> (x + y)**2
2
(x + y)
Symbol はいくつかの python の演算子で操作することができます: +, -, , * (算術演算), &, |, ~ , >>, << (論理演算).
Printing
ここで printing のために以下の設定をしておきましょう
>>> sympy.init_printing(use_unicode=False, wrap_line=True)
3.2.2. 代数的操作¶
Sympy は強力な代数操作を実行することができます. 最も頻繁に目にするであろうものとして: expand と simplify があります.
3.2.2.1. 展開¶
代数表現を展開するために使うことができます. 積やベキは降冪になるように試みます:
>>> expand((x + y)**3)
3 2 2 3
x + 3*x *y + 3*x*y + y
>>> 3*x*y**2 + 3*y*x**2 + x**3 + y**3
3 2 2 3
x + 3*x *y + 3*x*y + y
さらなるオプションをキーワードとして与えることができます:
>>> expand(x + y, complex=True)
re(x) + re(y) + I*im(x) + I*im(y)
>>> I*im(x) + I*im(y) + re(x) + re(y)
re(x) + re(y) + I*im(x) + I*im(y)
>>> expand(cos(x + y), trig=True)
-sin(x)*sin(y) + cos(x)*cos(y)
>>> cos(x)*cos(y) - sin(x)*sin(y)
-sin(x)*sin(y) + cos(x)*cos(y)
3.2.2.2. 簡単化¶
数式をより簡単な形式に変換する場合に simplify を使うことができます:
>>> simplify((x + x*y) / x)
y + 1
簡単化とはいくぶん曖昧な用語です, そのためより目的を明確にした simplify の代替が存在します: powsimp (指数の簡単化), trigsimp (三角関数を含む数式), logcombine, radsimp, togeter.
練習問題
\((x+y)^6\) の展開形を計算しましょう。
三角関数を含む式 \(\sin(x) / \cos(x)\) を簡単化しましょう
3.2.3. 微積分¶
3.2.3.1. 極限¶
極限は SymPy で簡単に計算することができ limit(function, variable, point) という構文に従います, つまり \(f(x)\) の \(x\rightarrow 0\) の極限を計算するには limit(f, x, 0) とします:
>>> limit(sin(x)/x, x, 0)
1
無限大での極限も計算することができます:
>>> limit(x, x, oo)
oo
>>> limit(1/x, x, oo)
0
>>> limit(x**x, x, 0)
1
3.2.3.2. 微分¶
Sympy の任意の数式は diff(func, var) を使って微分できます. 例:
>>> diff(sin(x), x)
cos(x)
>>> diff(sin(2*x), x)
2*cos(2*x)
>>> diff(tan(x), x)
2
tan (x) + 1
正しいかどうかは以下のようにして確認できます:
>>> limit((tan(x+y) - tan(x))/y, y, 0)
2
tan (x) + 1
高階微分は diff(func, var, n) メソッドで計算できます:
>>> diff(sin(2*x), x, 1)
2*cos(2*x)
>>> diff(sin(2*x), x, 2)
-4*sin(2*x)
>>> diff(sin(2*x), x, 3)
-8*cos(2*x)
3.2.3.3. 級数展開¶
Sympy はある点での数式の Taylor 展開をどう計算するかも知っています。series(expr, var):
>>> series(cos(x), x)
2 4
x x / 6\
1 - -- + -- + O\x /
2 24
>>> series(1/cos(x), x)
2 4
x 5*x / 6\
1 + -- + ---- + O\x /
2 24
練習問題
\(\lim{x\rightarrow 0, \sin(x)/x}\) を計算しましょう
\(log(x)\) の \(x\) 微分を計算しましょう。
3.2.3.4. 積分¶
SymPy は初等関数, 特殊関数の有限無限区間での積分も integrate() でサポートしています, これは強力な Risch-Norman の拡張アルゴリズムといくつかの発見的方法とパターンマッチングを利用しています。初等関数は以下のように積分できます:
>>> integrate(6*x**5, x)
6
x
>>> integrate(sin(x), x)
-cos(x)
>>> integrate(log(x), x)
x*log(x) - x
>>> integrate(2*x + sinh(x), x)
2
x + cosh(x)
特殊関数も容易に扱うことができます:
>>> integrate(exp(-x**2)*erf(x), x)
____ 2
\/ pi *erf (x)
--------------
4
有限区間での積分も計算できます:
>>> integrate(x**3, (x, -1, 1))
0
>>> integrate(sin(x), (x, 0, pi/2))
1
>>> integrate(cos(x), (x, -pi/2, pi/2))
2
広義積分もサポートしています:
>>> integrate(exp(-x), (x, 0, oo))
1
>>> integrate(exp(-x**2), (x, -oo, oo))
____
\/ pi
3.2.4. 方程式の求解¶
SymPy は1つまたはいくつかの変数についての代数方程式を解くことができます:
In [7]: solve(x**4 - 1, x)
Out[7]: [I, 1, -1, -I]
ここで見るように, 第1引数の数式は 0 と等しいと前提されます. 多くの多項式方程式を解くことができ, 連立方程式も変数の組を第2引数としてタプルを与えることで解くことができます:
In [8]: solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])
Out[8]: {y: 1, x: -3}
超越方程式も(限定的ですが)サポートされています:
In [9]: solve(exp(x) + 1, x)
Out[9]: [pi*I]
多項式の場合には factor が代替手段として利用できます. factor は多項式を既約でない項に因数分解し, 多くの領域で因数分解を計算できます:
In [10]: f = x**4 - 3*x**2 + 1
In [11]: factor(f)
Out[11]: (1 + x - x**2)*(1 - x - x**2)
In [12]: factor(f, modulus=5)
Out[12]: (2 + x)**2*(2 - x)**2
SymPy は真偽の論理式を解くこともできます, つまりある論理式が満される, または満されないかどうかを決定できます. そのためには satisfiable を使うことができます:
In [13]: satisfiable(x & y)
Out[13]: {x: True, y: True}
これは, (x & y) が True であるには x と y の両方が True である必要があるということを示しています. もし式が True になることが無い場合, つまり引数を True にできない場合 False を返します:
In [14]: satisfiable(x & ~x)
Out[14]: False
3.2.4.1. 練習問題¶
方程式系 \(x + y = 2\), \(2\cdot x + y = 0\) を解きましょう
(~x | y) & (~y | x)が真となるような真偽値x,yは存在しますか?
3.2.5. 線形代数¶
3.2.5.1. 行列¶
行列は Matrix クラスのインスタンスとして生成できます:
>>> from sympy import Matrix
>>> Matrix([[1,0], [0,1]])
[1 0]
[ ]
[0 1]
NumPy の配列と異なり, symbol を含むことができます:
>>> x = Symbol('x')
>>> y = Symbol('y')
>>> A = Matrix([[1,x], [y,1]])
>>> A
[1 x]
[ ]
[y 1]
>>> A**2
[x*y + 1 2*x ]
[ ]
[ 2*y x*y + 1]
3.2.5.2. 微分方程式¶
SymPy は (いくつかの) 常微分方程式を解くことができます。微分方程式を解くには dsolve を使います。まず cls=Function を symbols 関数に渡し、定義前の関数を作成します。
>>> f, g = symbols('f g', cls=Function)
f と g はここでは未定義の関数です。 f(x) を呼び出すこともできますが、未知の関数として表示されます:
>>> f(x)
f(x)
>>> f(x).diff(x, x) + f(x)
2
d
f(x) + ---(f(x))
2
dx
>>> dsolve(f(x).diff(x, x) + f(x), f(x))
f(x) = C1*sin(x) + C2*cos(x)
キーワード引数をこの関数に与えることで, 最適な解を見つけるための手助けができます. 例えば方程式が可分離であることを知っていればキーワードとして hint=’separable’ 使って dsolve に可分離な方程式として扱わせることができます:
>>> dsolve(sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x), f(x), hint='separable')
/ _____________\ / _____________\
| / C1 | | / C1 |
[f(x) = - asin| / ------- + 1 | + pi, f(x) = asin| / ------- + 1 | + pi,
| / 2 | | / 2 |
\\/ cos (x) / \\/ cos (x) /
/ _____________\ / _____________\
| / C1 | | / C1 |
f(x) = -asin| / ------- + 1 |, f(x) = asin| / ------- + 1 |]
| / 2 | | / 2 |
\\/ cos (x) / \\/ cos (x) /
練習問題
いくつかの Bernoulli 微分方程式の求解。
\[x \frac{d f(x)}{x} + f(x) - f(x)^2=0\]
同じ方程式を hint=’Bernoulli’ を使って解いてみましょう。何が起きますか?