読者です 読者をやめる 読者になる 読者になる

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

うまく解けない.