sympy
これまで数値計算はすべてmathematicaを用いてやってきましたが、webシステム開発者としてpythonに移行しようと思います。
まずは最初の基本として、
に記載された問題を実際に手を動かしてやって行こうと思います。
github
- jupyter notebook形式のファイルはこちら
google colaboratory
- google colaboratory で実行する場合はこちら
筆者の環境
筆者のOSはmacOSです。LinuxやUnixのコマンドとはオプションが異なります。
!sw_vers
ProductName: Mac OS X
ProductVersion: 10.14.6
BuildVersion: 18G6020
!python -V
Python 3.7.3
基本的なライブラリをインポートしそのバージョンを確認しておきます。
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib
import matplotlib.pyplot as plt
import scipy
import numpy as np
print('matplotlib version :', matplotlib.__version__)
print('scipy version :', scipy.__version__)
print('numpy version :', np.__version__)
matplotlib version : 3.0.3
scipy version : 1.4.1
numpy version : 1.16.2
sympyの3つの数値型
sympyは
- Real
- Rational
- Integer
の三つの型を持つようです。
分数の扱い
from sympy import *
a = Rational(1,2)
a
1/2
円周率の扱い
pi.evalf()
3.14159265358979
pi * 2
2*pi
pi ** 2
pi**2
自然対数の扱い
exp(1) ** 2
exp(2)
(exp(1) ** 2).evalf()
7.38905609893065
無限大
oo > 999
True
oo + 1
oo
evalfで桁数の指定
引数に桁数を入れるようです。
pi.evalf(1000)
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420199
対数
log(10)
log(10)
log(2,2**4).evalf()
0.250000000000000
E.evalf()
2.71828182845905
シンボルの定義
代数的な操作を可能にするために、$xや$yなどの変数を宣言します。
x = Symbol('x')
y = Symbol('y')
x + y ** 2 + x
2*x + y**2
(x + y) ** 3
(x + y)**3
代数操作
expandやsimplifyはmathematicaと同じのようで助かります。
展開
expand((x + y)**3)
x**3 + 3*x**2*y + 3*x*y**2 + y**3
簡易化
simplify((x + x*y) / x)
y + 1
因数分解
factor(x**3 + 3*x**2*y + 3*x*y**2 + y**3)
(x + y)**3
微積分
極限
limit(sin(x)/ x, x, 0)
1
limit(x, x, oo)
oo
limit(1 / x, x, oo)
0
limit(x**x, x, 0)
1
微分
diff(sin(x), x)
cos(x)
diff(sin(x) ** 2, x)
2*sin(x)*cos(x)
diff(sin(2 * x) , x)
2*cos(2*x)
diff(tan(x), x)
tan(x)**2 + 1
微分が正しいかの確認
極限をとり、微分が正しいかどうかチェックできます。微分の定義に従って計算させるだけです。
limit((tan(x+y) - tan(x))/y, y, 0)
tan(x)**2 + 1
高階微分も可能です。
diff(sin(x), x, 1)
cos(x)
2階微分で元の値にマイナスをかけた値になっている事がわかります。
diff(sin(x), x, 2)
-sin(x)
diff(sin(x), x, 3)
-cos(x)
diff(sin(x), x, 4)
sin(x)
級数展開
Taylor展開も可能です。
series(exp(x), x)
1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)
series(cos(x), x)
1 - x**2/2 + x**4/24 + O(x**6)
series(sin(x), x)
x - x**3/6 + x**5/120 + O(x**6)
第三引数で展開する次数を指定出来るかと思いやってみました。
series(exp(x), x, 6)
exp(6) + (x - 6)*exp(6) + (x - 6)**2*exp(6)/2 + (x - 6)**3*exp(6)/6 + (x - 6)**4*exp(6)/24 + (x - 6)**5*exp(6)/120 + O((x - 6)**6, (x, 6))
違うみたいで、第三引数に数値計算する際の中心の数で、第四引数で展開させる次数のようです。
series(exp(x), x, 0, 6)
1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)
積分
微分が出来れば、積分ももちろん対応しています。
integrate(x**3,x)
x**4/4
integrate(-sin(x),x)
cos(x)
integrate(exp(x),x)
exp(x)
integrate(log(x),x)
x*log(x) - x
integrate(exp(-x**2),x)
sqrt(pi)*erf(x)/2
積分区間も指定することが出来ます
integrate(x**2, (x, -1, 1))
2/3
integrate(sin(x), (x, 0, pi/2))
1
範囲が無限大の広義積分も可能です。
integrate(exp(-x), (x, 0, oo))
1
ガウス積分をやってみます。
integrate(exp(-x**2), (x, -oo, oo))
sqrt(pi)
integrate(exp(-x**2 / 3), (x, -oo, oo))
sqrt(3)*sqrt(pi)
方程式を解く
代数方程式を解くことも可能です。ほぼmathematicaと同じインターフェースです。
solve(x**3 - 1, x)
[1, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2]
solve(x**4 - 1, x)
[-1, 1, -I, I]
多変数の連立方程式も対応しています。
solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])
{x: -3, y: 1}
オイラーの式も計算してくれるようです。
solve(exp(x) + 1, x)
[I*pi]
行列演算
行列演算はnumpyでずっとやっているのでおそらく行列に関しては利用していませんが、一応ここでは手を動かして実行してみます。
from sympy import Matrix
a = Matrix([[1,0], [0,1]])
b = Matrix([[1,1], [1,1]])
a
1/2
b
Matrix([
[1, 1],
[1, 1]])
a + b
Matrix([
[2, 1],
[1, 2]])
a * b
Matrix([
[1, 1],
[1, 1]])
微分方程式
常微分方程式も解く事が可能です。dsolveを利用します。
f, g = symbols('f g', cls=Function)
f(x)
f(x)
f(x).diff(x,x) + f(x)
f(x) + Derivative(f(x), (x, 2))
dsolve(f(x).diff(x,x) + f(x), f(x))
Eq(f(x), C1*sin(x) + C2*cos(x))
実際に数値計算に利用されているのがよくわかる高機能ぶりでした。
もちろん有料であるmathematicaよりも劣る部分もあると思いますが、積極的にmathematicaからの移行を進めていきます。すごいなぁ〜
参考ページ
こちらのページを参考にしました。