pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

入力例で学ぶPython (SymPy) の使い方(入門)

本記事ではプログラミング言語Python』の数式処理ライブラリ『SymPy』の使い方を紹介します。

数式変形をすることそれ自体も楽しいものですが、単純作業という一面があることも否めません。効率よく勉強するには面倒な計算はコンピュータに任せてしまうのも一つの方法でしょう。 

数式を処理するソフトウェアとしてよく知られているのはMathematicaでしょうか?利用者も多いので、わからないことがあったときなどは、その充実したヘルプ機能に加えて、ネット検索によっても解決策が得られやすいという長所があります。しかしいかんせん有料(結構高い)なので個人利用で導入するのは気乗りしません。

ここではMathematicaの代替として、Pythonで使える数式処理ライブラリ『SymPy』の使い方を紹介します。といっても、いち利用者として具体例を扱いながら、習うより慣れろで使い方を学んだだけなので、ちょっと高度で便利な関数電卓の(自己流の)使い方を紹介します、という感覚の記事です。



SymPy の使い方入門

以下では具体的に Python(SymPy)の入力例を示します。

動作は JupyterNotebook で確認していますが、(もし問題があれば)適宜ご自分の動作環境のものに読み替えてください。

コピー&ペーストして計算を実行させたり、入力をいじったりして遊んでみてください。

SymPy ライブラリの読み込み

なにはともあれまずは SymPy ライブラリを読み込んでおきます:

import sympy

加減乗除

3+2 # ハッシュ記号の右側はコメント(評価の際に無視される文字列)になります

\quad\displaystyle{
3+2=5
}

3-2

\quad\displaystyle{
3-2=1
}

3*2

\quad\displaystyle{
3\times2=6
}

sympy.Rational(3,2)

\quad\displaystyle{
3\div2=\frac{3}{2}
}
また

3/2

と入力すると  1.5 と小数値を返します。

モジュロ演算(剰余演算)

11 % 2   #  m % n で、 m を n で割った余りを出力します

他方、商の整数部分は

11 // 2

で与えられます。

数値の型

整数型(int)
type(3)
浮動小数型(float)
type(3.)

冪乗・階乗・二重階乗・ルート

3**2

\quad\displaystyle{
3^2=9
}

10の冪乗については次のようにしても入力できます:

1.2e2

\quad\displaystyle{
1.2\times 10^2=120
}
ただしこのように入力した場合にはfloat型になります:

type(1.2e2)

また小数については

3e-3

\quad\displaystyle{
3\times10^{-3}=0.003
}

とします。

sympy.factorial(10) 

\quad\displaystyle{
10!=3628800
}

sympy.factorial2(10)

\quad\displaystyle{
10!!=3840
}

sympy.factorial2(9)

\quad\displaystyle{
9!!=945
}

sympy.sqrt(50)

\quad\displaystyle{
\sqrt{50}=5\sqrt{2}
}

複素数

虚数単位
sympy.I**2   # 虚数単位は大文字"アイ"

\quad\displaystyle{
i^2=-1
}

複素共役
sympy.conjugate(2+3*sympy.I)

\quad\displaystyle{
\overline{2+3i}=2-3i
}

実部・虚部
sympy.re(2+sympy.I*3)

\quad\displaystyle{
{\rm Re}(2+3i)=2
}

sympy.im(2+sympy.I*3)

\quad\displaystyle{
{\rm Im}(2+3i)=3
}

絶対値
sympy.Abs(2+sympy.I*3)

\quad\displaystyle{
\left|2+3i\right|=\sqrt{13}
}

偏角
sympy.arg(2+sympy.I*3)

\quad\displaystyle{
{\rm Arg}(2+3i)=\arctan\frac{3}{2}
}

代数的操作

素因数分解する
sympy.factorint(sympy.factorial(10))

\quad\displaystyle{
10!=2^8\cdot3^4\cdot5^2\cdot7^1
}

変数を宣言する
sympy.var('x') # 変数の宣言 別入力法として  x = sympy.symbols('x')  があります 
因数分解する
sympy.factor(x**2-4*x+3)

\quad\displaystyle{
x^2-4x+3=(x-3)(x-1)
}

部分分数分解する
sympy.apart(1/(x**2-4*x+3))

\quad\displaystyle{
\frac{1}{x^2-4x+3}=-\frac{1}{2(x-1)}+\frac{1}{2(x-3)}
}

sympy.apart(1/(x**6-1))

\quad\displaystyle{
\frac{1}{x^6-1}=\frac{x-2}{6(x^2-x+1)}-\frac{x+2}{6(x^2+x+1)}-\frac{1}{6(x+1)}+\frac{1}{6(x-1)}
}

多項式を展開する
sympy.var('x, y')
sympy.expand((x+y)**6)

\quad\displaystyle{
(x+y)^6=x^6+6x^5y+15x^4y^2+20x^3y^3+15x^2y^4+6xy^5+y^6
}

式を簡単化する・通分する
sympy.simplify((2*x-2)/((x-1)**2*(x-2)))

\quad\displaystyle{
\frac{(2x-2)}{(x-1)^2(x-2)}=\frac{2}{(x-2)(x-1)}
}

sympy.simplify(1/(x-1)**2+1/((x-1)*(x-2)))

\quad\displaystyle{
\frac{1}{(x-1)^2}+\frac{1}{(x-1)(x-2)}=\frac{2x-3}{(x-2)(x-1)^2}
}

式を簡単化・通分したうえで展開する
sympy.cancel((2*x-2)/((x-1)**2*(x-2)))

\quad\displaystyle{
\frac{(2x-2)}{(x-1)^2(x-2)}=\frac{2}{x^2-3x+2}
}

sympy.cancel(1/(x-1)**2+1/((x-1)*(x-2)))

\quad\displaystyle{
\frac{1}{(x-1)^2}+\frac{1}{(x-1)(x-2)}=\frac{2x-3}{x^3-4x^2+5x-2}
}

式の真偽を判定する
0==sympy.sqrt(2)-2**(sympy.Rational(1,2))
0==sympy.sqrt(2)-2**(sympy.Rational(1,3))
sympy.factorial(50)==sympy.factorial2(50)*sympy.factorial2(49)

ただし数式を適当に処理しないとうまく判定しない場合も多いです: 

x**2-2*x+1==(x-1)**2
x**2-2*x+1==sympy.expand((x-1)**2)
多項式から特定の項の係数を取り出す
sympy.expand((x+y)**6).coeff(x, 3) # expr.coeff(var, n) で数式 exp の変数 var の n 次の係数を表示
値を代入する
sympy.var('x, a, b, c')
(a*x**2+b*x+c).subs(x, 3)
sympy.expand((a*x**2+b*x+c).subs(x, (-b + sympy.sqrt(-4*a*c + b**2))/(2*a)))
総和
k=sympy.symbols('k',integer=True)
sympy.summation(k, (k, 1, 10) )

\quad\displaystyle{
\sum_{k=1}^{10}k=55
}

k, N=sympy.symbols('k N',integer=True)
sympy.factor(sympy.summation(k, (k, 1, N) ))

\quad\displaystyle{
\sum_{k=1}^{N}k=\frac{N(N+1)}{2}
}

sympy.factor(sympy.summation(k**3, (k, 1, N) ))

\quad\displaystyle{
\sum_{k=1}^{N}k^3=\frac{N^2(N+1)^2}{4}
}

総乗
k=sympy.symbols('k',integer=True)
sympy.product(k, (k, 1, 10) )

\quad\displaystyle{
\prod_{k=1}^{10}k=3628800
}

sympy.factorial(10)==sympy.product(k, (k, 1, 10) )

方程式

代数方程式の根を求める
sympy.var('x, a, b, c')
sympy.solve (a*x**2+b*x+c, x) # sympy.solve(eq, var) で方程式 eq=0 を変数 varについて解いた根を表示
連立方程式を解く
sympy.var('x, y')
sympy.solve ([3*x+4*y-5, 2*x+3*y-3], [x, y]) 
# sympy.solve([eq1, eq2, ...], [var1, var2, ...]) で連立方程式 eq1=0, eq2=0, ... の解を表示

より詳しくは
pianofisica.hatenablog.com
で解説しています。

関数の定義・変数への代入

文字に式を割り当てる、値を代入する

\quad\displaystyle{
u(x,y)=(x+y)
\qquad
v(x,y)=(x+y)^3
}

sympy.var('x, y')
u=x+y
v=(x+y)**3
u.subs(x,5)

\quad\displaystyle{
u(5,y)=y+5
}

v.subs(x,5).subs(y,3)

\quad\displaystyle{
v(5,3)=512
}

割り当てた文字を使う・四則演算
u+v
sympy.expand(u+v)
sympy.factor(u+v)
sympy.expand(u*v)
sympy.simplify(u/v)
割り当てた文字を使う・微分する
sympy.diff(v, x)   # v の x についての1階偏微分係数

\quad\displaystyle{
\frac{\partial v}{\partial x}=3(x+y)^2
}

sympy.diff(v, y, 2)   # u の y についての2階偏微分係数

\quad\displaystyle{
\frac{\partial^2 v}{\partial y^2}=6(x+y)
}

0==sympy.diff(v, x, 1)-sympy.diff(v, x)
割り当てた文字を使う・積分する
sympy.integrate(u, x)   # u を x について不定積分

\quad\displaystyle{
\int dx\, u=\frac{x^2}{2}+xy
}

sympy.integrate(v, (x, 0, 6))   # v を x について 0 から 6 までの区間で定積分

\quad\displaystyle{
\int_0^6 dx\, v=6y^3+54y^2+216y+324
}

関数を定義する
def f(x):
    return -sympy.log(1-x)/x

\quad\displaystyle{
f(x)=-\frac{\log(1-x)}{x}
}

f(2)

\quad\displaystyle{
f(2)=-\frac{1}{2}\log(-1)=-\frac{i\pi}{2}
}

関数を微分する
sympy.diff(f(x), x, 1)

\quad\displaystyle{
\frac{df}{dx}=\frac{1}{x(1-x)}+\frac{\log(1-x)}{x^2}
}

関数を積分する
sympy.integrate(f(x), x)

\quad\displaystyle{
\int dx\, f(x)={\rm Li}_2(x)
}

ここで  {\rm Li}_2 は二重対数関数(dilogarithm)と呼ばれる関数です。

あまり馴染みがないかもしれませんが、このすぐ後で簡単に紹介します。

Python(SymPy)で多重対数関数(polylogarithm)  {\rm Li}_s

sympy.polylog(s,x)

と入力します。

関数をTaylor展開する
sympy.series(f(x), x, 0, 7)  # f(x) を x について x=0 の周りで7次未満までTaylor展開

\quad\displaystyle{
f(x)=-\frac{\log(1-x)}{x}=1+\frac{x}{2}+\frac{x^2}{3}+\frac{x^3}{4}+\frac{x^4}{5}+\frac{x^5}{6}+\frac{x^6}{7}+\mathcal{O}(x^7)
}

ちなみに、ここまでの計算から

\quad\displaystyle{
{\rm Li}_2(x)= \int_0^x dx'\, f(x')=x+\frac{x^2}{2^2}+\frac{x^3}{3^2}+\frac{x^4}{4^2}+\frac{x^5}{5^2}+\frac{x^6}{6^2}+\frac{x^7}{7^2}+\mathcal{O}(x^7)
}

がわかります。一般に、多重対数関数は逐次的に

\quad\displaystyle{
{\rm Li}_{s+1}(x)= \int_0^x dx'\, \frac{{\rm Li}_{s}(x')}{x'}
}

で定義され、その  x=0 のまわりでのTaylor展開は

\quad\displaystyle{
{\rm Li}_s(x)=x+\frac{x^2}{2^s}+\frac{x^3}{3^s}+\frac{x^4}{4^s}+\frac{x^5}{5^s}+\frac{x^6}{6^s}+\frac{x^7}{7^s}+\mathcal{O}(x^7)
}

で与えられます。たとえば

sympy.series(sympy.polylog(3,x), x, 0, 4)

などです。ちなみに、変数を  x=1 と代入して  s の関数とみなすと

\quad\displaystyle{
{\rm Li}_s(1)=\sum_{n=1}^\infty \frac{1}{n^s}=\zeta(s)
}

ですが、この関数  \zeta(s) はRiemannのゼータ関数と呼ばれる整数論の重要な研究対象です。

等式の右辺・左辺、式の分母・分子を取り出す
F=x+1
G=x**2+x-2
E=sympy.Eq(x+1, F/G) # sympy.Eq(exp1, exp2) で exp1=exp2 という等式を表す
sympy.solve(E, x)
E.rhs  # 等式 E の右辺
E.lhs  # 等式 E の左辺
sympy.denom(E.rhs)  # 等式 E の右辺の分母
sympy.numer(E.rhs)  # 等式 E の右辺の分子

三角関数双曲線関数)の諸操作

引数の和を展開する
sympy.var('a, b')
sympy.expand(sympy.sin(a+b), trig=True)
sympy.expand(sympy.cos(a+b), trig=True)
sympy.expand(sympy.tan(a+b), trig=True)
簡単化する
sympy.trigsimp((sympy.sin(a))**2+(sympy.cos(a))**2)
sympy.trigsimp((sympy.sinh(a))**2+(sympy.cosh(a))**2)

微分方程式

常微分方程式を解く
sympy.var('a, x')
y = sympy.Function('y')(x)
eq = sympy.Eq( sympy.diff(y, x),  -a*y )
sympy.dsolve(eq)

より詳しくは
pianofisica.hatenablog.com
で解説しています。

リスト

リストを作る・要素を参照する
S = list([1, 10, 100])  # リストの定義
S[0]  # 要素の参照
S[0]+S[1]*S[2]
S[3] # 定義されていないものはエラー
リストの要素の和
sum(S)
リストの要素の最大値・最小値
max(S)
min(S)
リストの要素の個数
len(S)
リストを生成する(FOR文)
N_list = []  # 空のリストを準備する
def n(k): return k**2
for i in range(10):
    N_list.append(n(i))  # リストの各要素に関数 n(i) を加える
print(N_list)

条件分岐

Heaviside関数(階段関数、ステップ関数)

Heaviside関数は実軸上

\qquad\displaystyle{H(x)=\left\{\begin{array}{ll} \ 0 & (x<0) \\ \ 1 & \text {(otherwise)}\end{array}\right.}

と定義されます。これをPythonで自分で定義してみます:

def H(x):
   if x < 0: return 0
   else: return 1

より詳しい条件分岐の例については、以下の記事を参照してください:

pianofisica.hatenablog.com

グラフを描く

x = sympy.symbols('x')
sympy.plot(x**2)

変数の範囲を設定するには

x = sympy.symbols('x')
sympy.plot(sympy.tanh(x), (x, -5, 5))

とし、複数の対象を同時に描きたい場合には

x = sympy.symbols('x')
sympy.plot(sympy.log(x), sympy.sqrt(x), x, (x, -3, 3))

とします。より詳細は次の記事を読んでみてください:

pianofisica.hatenablog.com

TeX形式で計算結果を出力する

E2=sympy.Eq(a*x**2+b*x+c, 0)
sol=sympy.solve(E2, x)
sympy.init_printing()
display(sol[1]) # TeXでコンパイルして結果の数式を出力
print(sympy.latex(sol[0])) # TeXのソース形式で出力


さて、以上見てきたように SymPy で関数を扱うには、

いろいろな箇所で sympy. という文字列をおまじないのようにつけなければいけません。

これが面倒なようであれば、最初のライブラリの読み込みの段階で

from sympy import *

としておくという処方があるようです。

ただし、このとき、大文字のアイにはすでに虚数単位が割り当てられています:

I**2

このように、思いがけず使えない文字があったりするので、

SymPy を上記のようにして読み込んだ場合には注意が必要になりそうです。

または

import sympy as sp

として sympy と入力していた箇所を sp と省略することができます。


今回学んだ事項の例題として、Bernoulli数(ベルヌーイ数)を具体的に求める問題を扱っているのが

pianofisica.hatenablog.com

です。その(数学的にもPyhton的にも)発展編として

Bernoulli多項式(ベルヌーイ多項式)をPythonの計算を交えて解説しているのが

pianofisica.hatenablog.com

です。さらに、NumPy、Matplotlibライブラリを組み合わせた使い方について解説しているのが

pianofisica.hatenablog.com

です、本記事の続編、具体的な応用例として、是非あわせて読んでみてください。





キーワードPython、SymPy、JupyterNotbook