Sympyを試す
Sympyが面白かったので色々試しました。
このページはjupyter notebookを作ったものをmarkdownに変換したものです。もととなったnotebookはここにあります。
初期設定
# ここはとりあえずコピペすればよさそう
%matplotlib inline
from sympy import *
from sympy.plotting import plot
from IPython.display import display
init_printing()シンボルの定義
x = Symbol("x")
y = Symbol("y") # これでもいいけど...
x, y, z = symbols("x y z") # こっちのほうが楽なのでこっちでいいと思う
f, g = symbols('f g', cls=Function)
# displayはnotebook側の関数で、引数を良い感じに出力してくれる。
display(x, y, z)\(\displaystyle x\)
\(\displaystyle y\)
\(\displaystyle z\)
式の定義
expr = (x + y) ** 5
# display()を使うとnotebookを通じて出力できる
display(expr)
display(x * expr * y) # 変数の並び順は簡略化されるらしい\(\displaystyle \left(x + y\right)^{5}\)
\(\displaystyle x y \left(x + y\right)^{5}\)
式の展開/因数分解
expr = (x + y) ** 5
display(expr.expand()) # expand()で式を展開できる
display((x**3 + 3*x**2*y + 3*x*y**2 + y**3).factor()) # factor()はexpand()の逆\(\displaystyle x^{5} + 5 x^{4} y + 10 x^{3} y^{2} + 10 x^{2} y^{3} + 5 x y^{4} + y^{5}\)
\(\displaystyle \left(x + y\right)^{3}\)
式の比較
expr = (x + y) ** 2
expanded = expand(expr)
display(expr)
display(expanded)\(\displaystyle \left(x + y\right)^{2}\)
\(\displaystyle x^{2} + 2 x y + y^{2}\)
display(expr == expanded) # この比較は代数的な比較ではない
display(Eq(expr, expanded)) # Eq()を使えば代数的な比較を作れる
display(Eq(expr, expanded).simplify()) # さらにSimplify()で恒等的に成り立つかどうかをチェックできるFalse\(\displaystyle \left(x + y\right)^{2} = x^{2} + 2 x y + y^{2}\)
\(\displaystyle \mathrm{True}\)
expr.equals(expanded) # equals()はランダムな点でテストをして等しいかどうかを判定するらしいTrue代入
expr = sin(x) + 1
expr.subs(x, pi / 3)\(\displaystyle \frac{\sqrt{3}}{2} + 1\)
expr = x**2 + 3*x + 4
expr.subs(x, expr) # 式の代入も可能\(\displaystyle 3 x^{2} + 9 x + \left(x^{2} + 3 x + 4\right)^{2} + 16\)
expr = x**3 + 5*x*y + x*z + 5*y**2 + y*z**2 + z**2
display(expr)
display(expr.subs([(x, 1), [y, 3], [z, 2]]))
display(expr.subs([(x, x), [y, x], [z, x]]))\(\displaystyle x^{3} + 5 x y + x z + 5 y^{2} + y z^{2} + z^{2}\)
\(\displaystyle 79\)
\(\displaystyle 2 x^{3} + 12 x^{2}\)
無理数の計算
pi.evalf(30) # 1000ぐらい余裕\(\displaystyle 3.14159265358979323846264338328\)
expr = sin(1)**2 + cos(1)**2
display((expr - 1).evalf())
display((expr - 1).evalf(chop = True))\(\displaystyle -4.0 \cdot 10^{-124}\)
\(\displaystyle 0\)
numpy関数の利用
import numpy
a = numpy.arange(100)
expr = sin(x)**2 + cos(x)**3
display(expr)
f = lambdify(x, expr, "numpy") # numpyを使わないと遅い
f(a)\(\displaystyle \sin^{2}{\left (x \right )} + \cos^{3}{\left (x \right )}\)
array([ 1. , 0.86580202, 0.75475425, -0.95036208, 0.29348079,
0.94236043, 0.96327991, 0.86012577, 0.97574947, -0.58654075,
-0.29478182, 0.9999805 , 0.88881004, 0.92378616, 0.98385952,
-0.01556116, -0.79539238, 0.90345118, 0.85189192, 0.98895835,
0.90142733, 0.53566986, -0.99980414, 0.5648118 , 0.89639378,
0.99135691, 0.85223411, 0.88972231, -0.81857033, 0.02180437,
0.97987667, 0.92866036, 0.88463129, 0.99982139, -0.33110279,
-0.55466852, 0.98152993, 0.86256579, 0.95901972, 0.94785945,
0.25853532, -0.93733553, 0.7760188 , 0.86290782, 0.99984336,
0.86900634, 0.73250099, -0.96191188, 0.32789382, 0.93680443,
0.96737779, 0.85798037, 0.96910399, -0.61757906, -0.25811662,
0.99952124, 0.89314651, 0.91891613, 0.98748892, -0.0530218 ,
-0.77101524, 0.91618973, 0.85189896, 0.98629159, 0.90661841,
0.50571095, -0.99823778, 0.59310782, 0.89154121, 0.99347747,
0.85292682, 0.87499481, -0.84050791, 0.05901876, 0.97557656,
0.93351789, 0.88062925, 0.99901083, -0.36702009, -0.5220171 ,
0.98646955, 0.8652878 , 0.95461493, 0.95326866, 0.22310345,
-0.92285576, 0.79628569, 0.86033431, 0.99937407, 0.87250858,
0.70927031, -0.97196403, 0.3617301 , 0.93122345, 0.97129634,
0.85614111, 0.96157091, -0.64772986, -0.22116672, 0.99847744])極限/微分/積分
expr = sin(x)/x
# 極限をその場で計算しないならLimitを、するならlimitを使う(以下同様に2種類の関数がある)
display(Eq(Limit(expr, x, 0), limit(expr, x, 0)))
# 極限の方向は片方しか指定できないみたい
display(Eq(Limit(expr, x, 0, dir="-"), limit(expr, x, 0, dir="-"))) \(\displaystyle \lim_{x \to 0^+}\left(\frac{1}{x} \sin{\left (x \right )}\right) = 1\)
\(\displaystyle \lim_{x \to 0^-}\left(\frac{1}{x} \sin{\left (x \right )}\right) = 1\)
expr = exp(sin(x))
Eq(Derivative(expr), diff(expr)) \(\displaystyle \frac{d}{d x} e^{\sin{\left (x \right )}} = e^{\sin{\left (x \right )}} \cos{\left (x \right )}\)
expr = x * (exp(x)*sin(x) + exp(x)*cos(x))
Eq(Integral(expr), integrate(expr))\(\displaystyle \int x \left(e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )}\right)\, dx = x e^{x} \sin{\left (x \right )} - \frac{e^{x}}{2} \sin{\left (x \right )} + \frac{e^{x}}{2} \cos{\left (x \right )}\)
expr = exp(-x**2)
# 無限大は"oo"で指定
Eq(Integral(expr, (x, -oo, oo)), integrate(expr, (x, -oo, oo)))\(\displaystyle \int_{-\infty}^{\infty} e^{- x^{2}}\, dx = \sqrt{\pi}\)
級数
expr = exp(x)
Eq(expr, expr.series(x=x, x0=0, n=10))\(\displaystyle e^{x} = 1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \frac{x^{5}}{120} + \frac{x^{6}}{720} + \frac{x^{7}}{5040} + \frac{x^{8}}{40320} + \frac{x^{9}}{362880} + \mathcal{O}\left(x^{10}\right)\)
expr = exp(sin(x)**2)
Eq(expr, expr.series(x=x, x0=0, n=10))\(\displaystyle e^{\sin^{2}{\left (x \right )}} = 1 + x^{2} + \frac{x^{4}}{6} - \frac{11 x^{6}}{90} - \frac{71 x^{8}}{2520} + \mathcal{O}\left(x^{10}\right)\)
プロット
x = symbols('x')
plot(x*x)
plot(sin(x))<sympy.plotting.plot.Plot at 0x271b9e4d828>