yukiのブログ

作ったものなど

Sympyを試す

 

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))

png

png

<sympy.plotting.plot.Plot at 0x271b9e4d828>