pianofisica

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

トップ > Python > 入力例で学ぶPython (SymPy) の使い方（入門）

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

Python Python-SymPy Python-SymPy-基本操作

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

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

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

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



SymPy の使い方入門

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

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

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

SymPyライブラリの使い方の前に、Pythonの基本的な構文などを確認しておきたい方は以下の記事を参照してください。

pianofisica.hatenablog.com



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.pi
自然対数の底（ネイピア数）
sympy.E
オイラーの公式
sympy.E**(sympy.pi*sympy.I)

\quad\displaystyle{ e^{i\pi}=-1 }

無限大
sympy.oo

数値の参照

数学定数や、特定の点での関数の値などを数値的にみたいときには

sympy.pi.evalf(10)

などとします。evalfの引数は表示する桁数を指定します。

代数的操作

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

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

変数を宣言する
sympy.var('x') # 変数 x の宣言

または

x = sympy.symbols('x')

としても 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 次の係数を表示
変数に値を代入する

２次式

\quad\displaystyle{ ax^2+bx+c }

で、 x=3 を代入してみます。

sympy.var('x, a, b, c')
(a*x**2+b*x+c).subs(x, 3)

代入できるものは数字に限らず、文字式でも問題ありません。２次式

\quad\displaystyle{ ax^2+bx+c }

に、解の公式から得られる根の１つ

\quad\displaystyle{ x=\frac{-b+\sqrt{b^2-4ac}}{2a} }

を代入して整理してみましょう。

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

変数の性質を規定する

整数であると宣言する

すでに上の総和を求めるさいの添字として、変数が整数であると宣言する書式を使っています。

k = sympy.symbols('k', integer = True)
実数であると宣言する
sympy.var('a, b', real = True)
sympy.conjugate(a+b*sympy.I)

とくに何も宣言がなければ、文字は複素数として扱われてしまいます：

sympy.var('z')
sympy.conjugate(z)
非負であると宣言する
sympy.var('r', nonnegative = True)
sympy.sqrt(r**2)

とくに何も宣言がなければ、絶対値などでケアが必要な場合に処理が進みません：

sympy.var('s')
sympy.sqrt(s**2)
正であると宣言する
sympy.var('t', positive = True)
sympy.Abs(t)
変数依存性を宣言する
from sympy.core.function import Function
R = Function('R')(r, t)

したがって

sympy.diff(R, r)

によって

\quad\displaystyle{ \frac{\partial}{\partial r}R(r,t) }

などという計算ができます。





初等関数

三角関数

正弦関数（サイン）

sympy.var('x')
sympy.sin(x)

余弦関数（コサイン）

sympy.cos(sympy.pi)

正接関数（タンジェント）

sympy.tan(sympy.pi/2)

正割関数（セカント）

sympy.sec(x)

余割関数（コセカント）

sympy.csc(x)

余接関数（コタンジェント）

sympy.cot(x)
逆三角関数
sympy.asin(1) 

sympy.acos(1) 

sympy.atan(1)
双曲線関数
sympy.sinh(x) 

sympy.cosh(x) 

sympy.tanh(x)
指数関数
sympy.exp(1)

指数関数と三角関数の関係は次のオイラーの公式の一般形

\quad\displaystyle{ e^{i\theta}=\cos\theta +i\sin\theta }

で与えられます。ただし \theta は実数です。これをみてみます。左辺の実部は

sympy.var('theta', real = True)
sympy.re(sympy.exp(sympy.I*theta))

より \cos\theta であることが、虚部は

sympy.im(sympy.exp(sympy.I*theta))

より \sin\theta であることが出力されます。

対数関数
sympy.log(sympy.E)

指数関数と対数関数は互いに逆関数の関係です。実際、正の引数 \theta に対して

sympy.var('theta', real = True)
sympy.exp(sympy.log(theta)) 

sympy.log(sympy.exp(theta))

より

\quad\displaystyle{ \theta=\exp\left(\log(\theta)\right)=\log\left(\exp(\theta)\right) }

が成り立ちます。しかし２つめの等号については対数関数の多価性をケアして、変数の性質を規定しない場合には勝手に処理されることはありません。

sympy.var('z')
sympy.log(sympy.exp(z))

三角関数（双曲線関数）の諸操作

引数の和を展開する（加法公式）
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)

方程式

代数方程式の根を求める

２次方程式

\quad\displaystyle{ ax^2+bx+c=0 }

の解はSymPyの組み込み関数 solve で求めることができます。

sympy.var('x, a, b, c')
sympy.solve (a*x**2+b*x+c, x) # sympy.solve(eq, var) で方程式 eq=0 を変数 varについて解いた根を表示
連立方程式を解く

組み込み関数 solve は連立方程式の場合にも使えます。

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^2-y^2 \qquad v(x,y)=2xy }

sympy.var('x, y')
u = x**2-y**2
v = 2*x*y
変数に値を代入する
u.subs(x,5)

\quad\displaystyle{ u(5,y)=25-y^2 }

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

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

割り当てた文字を使う・四則演算
u+v

sympy.expand(u*v)

sympy.factor(u*v)

sympy.simplify(u/v)
割り当てた文字を使う・微分する
sympy.diff(v, x)   # v の x についての１階偏微分係数

\quad\displaystyle{ \frac{\partial v}{\partial x}=2y }

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

\quad\displaystyle{ \frac{\partial^2 v}{\partial y^2}=0 }

0==sympy.diff(v, x, 1)-sympy.diff(v, x)

いま扱っている関数の具体形に対して

sympy.diff(u, x, 1)==sympy.diff(v, y, 1)

sympy.diff(u, y, 1)==-sympy.diff(v, x, 1)

より

\quad\displaystyle{ \frac{\partial u}{\partial x}=\frac{\partial v}{\partial y}\qquad\quad \frac{\partial u}{\partial y}=-\frac{\partial v}{\partial x} }

という関係式が成り立ちます。これはCauchy-Riemann（コーシー・リーマン）の関係式と呼ばれ、複素関数論で基本となる関係式です。詳細は述べませんが、この関係式は以下のような考察から導かれます。いま２次元平面上の点 (x,y) を複素数 z=x+iy とみなし、関数値の組 (u,v) を同様に複素数 w=u+iv とみなしたとき、これを複素数 z から複素数 w へ移す関数 w=f(z) とみなすことができます。ここで変数変換 (x,y)\mapsto(z,\bar{z}) を考えたとき、ここで定めた関数 f が \bar{z} によらない（正則である）ための必要十分条件としてCauchy-Riemann関係式が得られます。

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

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

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

\quad\displaystyle{ \int_0^6 dx\, v=36y }

いま経路に依存する線積分

\quad\displaystyle{ I(C)=\int_C \left(udx-vdy\right) }

を考えてみます。はじめに、経路 C_1 として、原点 (0,0) から x-軸上を [0,1] の区間移動し、そのあと、 y-軸方向に [0,1] の区間に沿って移動する経路をとってみましょう。つまり

\quad\displaystyle{ C_1\:\ (0,0)\quad\overset{(x,0)}{\longrightarrow}\quad(1,0)\quad\overset{(1,y)}{\longrightarrow}\quad(1,1) }

に対して上記の積分を考えてみます。

sympy.integrate(u.subs(y, 0), (x, 0, 1))-sympy.integrate(v.subs(x, 1), (y, 0, 1))

より

\quad\displaystyle{ I(C_1)=-\frac{2}{3} }

次に、２つめの経路 C_2 としてを原点 (0,0) から y=x の直線上を (1,1) の点まで移動してみます。つまり

\quad\displaystyle{ C_2\:\ (0,0)\quad\overset{(x,y)_{x=y}}{\longrightarrow}\quad(1,1) }

を考えます。

sympy.integrate(u.subs(y, x), (x, 0, 1))-sympy.integrate(v.subs(x, y), (y, 0, 1))

より

\quad\displaystyle{ I(C_2)=-\frac{2}{3} }

さいごに３つめの経路 C_3 としてを原点 (0,0) から y-軸上を [0,1] の区間移動し、そのあと、 x-軸方向に [0,1] の区間に沿って移動する経路をとってみましょう。つまり

\quad\displaystyle{ C_3\:\ (0,0)\quad\overset{(0,y)}{\longrightarrow}\quad(0,1)\quad\overset{(x,1)}{\longrightarrow}\quad(1,1) }

を考えます。すると

sympy.integrate(u.subs(y, 1), (x, 0, 1))-sympy.integrate(v.subs(x, 0), (y, 0, 1))

より、やはり

\quad\displaystyle{ I(C_3)=-\frac{2}{3} }

となります。まとめると、いま３つの異なる経路に対して線積分を考えその値を計算してみましたが、その値は積分経路によらないことがわかりました。もちろん３つしか確認していないので、数学的な証明ではありませんが、別な経路を考えて試してみてください。やはり同じ値になるはずです。実は、これは上で述べたCauchy-Riemannの関係式が成り立っていることに起因する性質です。これは「Cauchyの積分定理」として知られています。

関数を定義する
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{1}{2}\log\,e^{i\pi}=-\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 の周りで７次未満まで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^8) }

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

\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^8) }

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

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(F/G)  # F/G の分母

sympy.numer(F/G)  # F/G の分子

極限・広義積分

極限値を求める
sympy.limit(sympy.sin(x)/x, x, 0)

\quad\displaystyle{ \lim_{x\to0}\frac{\sin(x)}{x}=1 }

広義積分を求める
sympy.integrate(1/(sympy.sqrt(sympy.pi))*sympy.exp(-x**2), (x, -sympy.oo, sympy.oo) )

\quad\displaystyle{ \int_{-\infty}^\infty dx\ \frac{e^{-x^2}}{\sqrt{\pi}}=1 }




行列

行列を入力する
sympy.var('a11, a12, a21, a22, b11, b12, b21, b22')
A = sympy.Matrix([
[a11,a12],
[a21,a22]])
B = sympy.Matrix([
[b11,b12],
[b21,b22]])
行列の要素を参照する
A[0,0]

Python ではインデックスのカウントがゼロから始まるので、数学でいうところの行列の (i,j)-成分が、[i-1,j-1]-成分と読み替える必要があります。慣れないうちは、この点に注意が必要です。

行列の和を求める
sympy.var('c, d')
c*A+d*B
行列の積を求める
A*B

行列の積は通常の数の積と同じ記号です。

行列の冪乗を求める
A**3
行列式を求める
A.det()
逆行列を求める
A**(-1)

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

微分方程式

常微分方程式を解く

１階線形微分方程式
\quad\displaystyle{ \frac{dy}{dx}=-ay }
はPython（SymPy）の組み込み関数 dsolve を使って答えを得ることができます。

sympy.var('a, x')
y = sympy.Function('y')(x)
eq = sympy.Eq( sympy.diff(y, x),  -a*y )
sympy.dsolve(eq)

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

特殊関数

直交多項式

pianofisica.hatenablog.com

超幾何関数

pianofisica.hatenablog.com

入力一覧

pianofisica.hatenablog.com


Fourier級数

Fourier級数展開を求める

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)

数値計算

数式処理ライブラリであるSymPyは数値的な計算をするためにはあまり向いていません。数値計算については、以下の別な記事で取り上げていますので、あわせて読んでみてください。

桁落ち

pianofisica.hatenablog.com

方程式を数値的に解く

pianofisica.hatenablog.com

微分方程式を数値的に解く

pianofisica.hatenablog.com

数値的に積分を求める

pianofisica.hatenablog.com

確率・統計

pianofisica.hatenablog.com
pianofisica.hatenablog.com

モンテカルロ法

pianofisica.hatenablog.com





条件分岐

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

直前の出力結果を参照する

直前の出力結果は" _ "（アンダースコア）で参照できます。いま行列Aの逆行列を出力して、それと行列Aとの積をとってみます。

sympy.var('a11, a12, a21, a22, b11, b12, b21, b22')
A = sympy.Matrix([
[a11,a12],
[a21,a22]])
A**(-1)

と入力したあと、その結果を使う場合に

_*A

などとします。

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.var('x, y')
exp1 = sympy.log((1+x)*(1+y)**2)
exp2 = sympy.expand(exp1)
sympy.print_python( exp2 )

定義した文字を初期化（無効化）する

del(x)

定義した文字をすべて初期化（無効化）する

%reset


さて、以上見てきたように SymPy で関数を扱うには、いろいろな箇所で sympy. という文字列をおまじないのようにつけなければいけません。これが面倒なようであれば、最初のライブラリの読み込みの段階で

from sympy import *

としておくという処方があります。

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

I**2

このように、思いがけず使えない文字があったりするので、SymPy を上記のようにして読み込んだ場合には注意が必要になります。または

import sympy as sp

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


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

pianofisica.hatenablog.com

です。その（数学的にもPyhton的にも）発展編としてBernoulli多項式（ベルヌーイ多項式）をPythonの計算を交えて解説しているのが

pianofisica.hatenablog.com

です。また、より実践的にSymPyを数式処理に活用しているのが次の記事です：

pianofisica.hatenablog.com

さらに、NumPy、Matplotlibライブラリを組み合わせた使い方についての解説は

pianofisica.hatenablog.com

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





キーワード：Python、SymPy、JupyterNotbook

プライバシーポリシー