MyEnigma

とある自律移動システムエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

Pythonの数式処理ライブラリSymPyをWolfram Alpha(Mathematica, Maxima)の代わりに使う方法

目次

 

Wolfram Alphaの問題点

以前、Wolfram Alphaがすごいという記事を書きましたが、

myenigma.hatenablog.com

myenigma.hatenablog.com

不満もあります。

それは複数行のコードを読みこませることができなかったり、

以前の計算の結果を使って、再計算するようなことができないことです。

また、計算されたデータなどを自由にグラフ化することもできません。

 

上記のような問題を解決するためには、

Wolfram Alphaの元になっているMathematicaや、

Mathematicaのフリー版であるMaximaを使うのがいいですが、

それぞれ独自の言語実装なので、

少しハードルが高くなってしまいます。

 

そこで今回おすすめなのが

Pythonの数式処理ライブラリSymPyです。

f:id:meison_amsl:20151121222658p:plain

 

今回はこのSymPyについて説明したいと思います。

 

Sympyとは?

SympyはPython製の数式処理モジュールです。

 

式の展開から、微積分、

方程式の解など、

いわゆるMathematicaやWolfram Alphaができる

代数処理をPythonで実現することができます。

 

また、SymPyはpythonのデフォルト機能のみで作成されており、

他のライブラリに依存していないため、

利用しやすいという特徴があります。

 

ブラウザでSympyを使う

SymPyのHPには、

Sympyを試すことができるブラウザツールがあります。

 

これが非常に便利で、

少し数式を確認したいだけであれば、

このブラウザSympyのみで十分だと思います。

f:id:meison_amsl:20151121205625p:plain

 

ブラウザ内でフルスクリーンモードにしたり、

コマンドラインを初期化したりできます。

一番うれしいのは上記の図のように、

計算された数式をLatexを使って描画してくれるところです。

スクリーンショットを撮れば、

プレゼンなどにも使えるのでいいですね。

 

ちなみにこのブラウザSymPyは

Google App Engine上で動いているみたいです。

 

JupyterでSymPyを使う

下記の記事のように、

Pythonのインタラクティブ実行環境Jupyter Notebookで

SymPyを使うこともできます。

myenigma.hatenablog.com

 

インストール

基本的にはpipでインストールできます。

pip install sympy

 

また、pythonのライブラリをひとまとめにした

Anacondaにもsympyが入っているので、

こちらをインストールしても、SymPyが使えるようになります。

Download Anaconda now! | Continuum

 

Sympyの数式処理の例題

Sympyを使った数式処理の例を説明します。

ちなみに下記のサンプルコードでは、

from sympy import *
x = Symbol('x')
y = Symbol('y')

のようにSymPyがインポートされ、

x,yなどの変数が指定されていることを前提にしています。

前述のブラウザSympyでは、上記のインポートは不要です。  

多項式展開

expandを使うと式を展開してくれます。

expr = (x + y)**5
expand(expr)

f:id:meison_amsl:20151121212823p:plain

 

方程式をある変数で解く

solve関数で、解きたい変数の指定することで、

指定の変数で式を変形してくれます。

expr = cos(x) * ln(y) + 2/y
solve(expr, x)

f:id:meison_amsl:20151122103359p:plain

連立方程式を解く

前述のsolve関数を使うのですが、

複数の連立方程式を解きたい場合は、

それぞれの方程式と、解きたい変数をリストにして

solveに渡すことで解くことができます。

x, y, z = symbols('x y z')
eq1=x * y * z + 234
eq2=x + y + z - 20
eq3=5 * x - y + 2 * z - 85
solve([eq1,eq2,eq3], [x,y,z])

f:id:meison_amsl:20160504215513p:plain

 

数値を代入する

subs関数を使うことで、

変数に数値を代入し、その結果を表示できます。

f = x**2 + 3*x + 2 
f1 = f.subs([(x, 1)])

f:id:meison_amsl:20151122110120p:plain

 

微分

diffという関数を使います。

x,y = symbols('x y')
f = x**2 / y + 2 * x - ln(y)
diff(f,x)

f:id:meison_amsl:20151121221226p:plain

 

積分

integrate関数を使うと積分してくれます。

expr = cos(x) * ln(y) + 2/y
integrate(expr, x)

f:id:meison_amsl:20151122104548p:plain

 

テイラー展開

series関数を使うと、

指定した式のテイラー展開をすることができます。

sin(x).series(x)

f:id:meison_amsl:20151121215124p:plain

 

極限

limit関数を使います。

一つ目の引数が極限を求めたい関数

二つ目の引数が先ほどの関数の中で、極限を求めたい変数

三つ目の引数が、先ほどの変数をどの値に限りなく近づけるかです。

ちなみに無限大∞を表すには、アルファベットのo(オー)を2つ並べて表します。

limit(sin(x)/x, x, 0)
limit(sin(x)/x, x, oo)

f:id:meison_amsl:20151122212422p:plain

 

その他

その他のSymPyの使い方に関しては、

下記のSymPyのドキュメントや

下記のサイトではSymPyの使い方をランダムに表示してくれるので、

使い方の勉強ができます。

 

Sympyを使ったエンジニアリング例題

実際にSympyを使った数式計算をしてみたいと思います。

三次元の回転行列を計算する

下記の記事における三次元の3つ回転行列を掛けあわせた行列を計算してみます。

myenigma.hatenablog.com

 

下記のコードを実行すると、

theta=Symbol("theta")
chi=Symbol("chi")
phi=Symbol("phi")
Rx=Matrix([[1,0,0],[0,cos(theta),-sin(theta)],[0,sin(theta),cos(theta)]])
Ry=Matrix([[cos(chi),0,sin(chi)],[0,1,0],[-sin(chi),0,cos(chi)]])
Rz=Matrix([[cos(phi),-sin(phi),0],[sin(phi),cos(phi),0],[0,0,1]])
print(Rx*Ry*Rz)

下記が表示されます。

Matrix([[cos(chi)cos(phi), -sin(phi)cos(chi), sin(chi)], [sin(chi)sin(theta)cos(phi) + sin(phi)cos(theta), -sin(chi)sin(phi)sin(theta) + cos(phi)cos(theta), -sin(theta)cos(chi)], [-sin(chi)cos(phi)cos(theta) + sin(phi)sin(theta), sin(chi)sin(phi)cos(theta) + sin(theta)cos(phi), cos(chi)cos(theta)]])

先ほどのSymPy Liveを使うとこのように数式が表示されます。

(行列の一つの項が長いので、かぶってしまっていますが。。。)

f:id:meison_amsl:20151121210818p:plain

 

先ほどの記事内での資料と同じ数式になっているのがわかります。

これで三次元空間での回転行列の式が計算できました。

 

SymPyによるCコードの出力

実は、SymPyには、

SymPyで表現された数式で計算を実施する

Cコードを自動生成する機能があります。

Codgenというモジュールを使います。

  

下記がサンプルコードです。

前述のブラウザSymPyで計算するときにも、

モジュールのインポートが必要になります。

 

from sympy.utilities.codegen import *
[(c_name, c_code), (h_name, c_header)] = codegen(("f", sqrt(x*x+y*y+z*z)), "C", "test")
print c_code
print c_header

 

c_codeとc_headerに生成されたCコードが入っています。

表示させると下記のようになるはずです。

f:id:meison_amsl:20151122114204p:plain

f:id:meison_amsl:20151122114210p:plain

  

codegenの二つ目の引数がSymPyによる数式

4つ目の引数がクラス名になります。

 

SymPyはJupyterと一緒に使うとかなり便利

SymPyはターミナル上で実行しても良いですが、

Jupyter Notebook上で実行すると、

計算の過程などがすべて残って非常に便利です。

f:id:meison_amsl:20170513152411p:plain

  

Jupyter Notebookに関しては下記を参照下さい。

myenigma.hatenablog.com

 

一つ注意点として、

Jupyter Noterbook上でSymPyを使用して

計算された数式などをJupyter上に表示させたい場合は、

Jupyter Notebook上で下記のコマンドを実行しておく必要があります。

import sympy
from sympy import init_printing
init_printing()

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

 

MyEnigma Supporters

もしこの記事が参考になり、

ブログをサポートしたいと思われた方は、

こちらからよろしくお願いします。

myenigma.hatenablog.com