sympy で遊ぶ

目的

とりあえず動かす, というところまでの how to な記事が少なくいろいろ苦労したので, そのあたりのことを書いておきたい.

スクリプトから API 叩くというだけなら sympy 入れたら良いだけなのだが, インタラクティブな UI で遊ぶことを目指す.

概要

ipython notebook するとウェブサーバが上がってデフォルトブラウザでそいつを開くので, ここで sympy を使って微分方程式解いたり plot() したりする.

導入

Gentoo

# emerge numpy scipy sympy ipython jinja www-servers/tornado jsonschema

jinja やら tornado やらは ipython notebook に必要なのだが, 依存関係に入ってないので, 手動で入れる必要がある. なにか他にも必要かもしれない.

Homebrew

$ brew tap homebrew/science
$ brew tap homebrew/python
$ brew install numpy
$ brew install scipy
$ pip3 install sympy
$ pip3 install ipython
$ pip3 install notebook
$ pip3 install matplotlib

pip の場合 notebookipython と別になってるので, notebook に必要なものは一通り勝手に入るらしい.

起動

$ ipython3 notebook

とすると

f:id:crckyl:20151110210239p:plain

こういうのが起動するので, Notebook を新規作成する.

Pretty print

数式を表示するときに, 裏で latex を叩いて綺麗に表示してくれる機能がある. ぐぐる%load_ext ... なるコマンドを入れたら良いとの記述が見つかるが, この方法は古いらしい.

init_printing()

とすると綺麗になる.

f:id:crckyl:20151110210837p:plain

from sympy import *
init_printing()
t,m,k=symbols('t m k')
x=Function('x')
eq=Eq(m*x(t).diff(t,2),-k*x(t))
eq

おなじみのバネです.

微分方程式を解く

f:id:crckyl:20151110211311p:plain

classify_ode(eq)
x2=dsolve(eq)
x2

解けないときは classify_ode()() になる. 偏微分方程式classify_pde(), pdsolve(). Equality オブジェクトでなく式でも良い. その場合 =0 ということになるらしい.

プロットする

f:id:crckyl:20151110211629p:plain

x3=x2.subs(zip(symbols('C1 C2 k m'),(1,1,1,1)))
x3=x3.simplify()
x3
%matplotlib inline
plot(x3.rhs)

%matplotlib inline するとインラインに表示してくれる. pyglet かなにかが入ってると別ウィンドウで出してくる. 入ってないと何も出ない.

フーリエ級数展開

f:id:crckyl:20151111193256p:plain

from sympy.mpmath import *
r=[-2,2]
f=lambda x:1 if -1<x<1 else 0
def g(n):
    cs=fourier(f, r, n)
    return lambda x: fourierval(cs, r, x)
plot([f,g(1),g(3),g(10),g(50)],r)

どうやらこれは数値計算になるのかな.

f:id:crckyl:20151111200500p:plain

def fourier_series(f, arg, range_, N):
  T = range_[1] - range_[0]
  a, b = [], []
  for k in range(N):
    a_k = 2/T*integrate(f*cos(2*pi*(k+1)/T*arg), [arg] + range_)
    b_k = 2/T*integrate(f*sin(2*pi*(k+1)/T*arg), [arg] + range_)
    a.append(a_k)
    b.append(b_k)
    pass
  a_0=integrate(f, [arg] + range_)/T
  return a_0 + sum([a_k*cos(2*pi*(k+1)/T*arg) + b_k*sin(2*pi*(k+1)/T*arg)
                     for k, (a_k, b_k) in enumerate(zip(a, b))])

t = symbols('t')
f = Piecewise((0, t<-1), (0, t>1), (1, True))
g = lambda n: fourier_series(f, t, [-2, 2], n)
plot(*[g(n*2+1) for n in (1,2,5,20)],xlim=(-2,2))

というわけで手書き.

交流 RC 回路

f:id:crckyl:20151111224855p:plain

うまく解けない.