hans

hans

【Python】四、Sympy——《Pythonを使った科学計算》


SymPy は Python の数学記号計算ライブラリで、これを使うことで数学式の記号的推導ができます。呼び出しを便利にするために、以下のすべての例プログラムでは事前に sympy ライブラリからすべての内容をインポートしたと仮定します:

>>> from sympy import *

* はすべてを意味します


1. 古典的な公式#

1668632086147.jpg
これはオイラーの等式と呼ばれ、ここで e は自然指数の底、i は虚数単位、pi は円周率です。この公式は数学で最も驚くべき公式と称賛されており、5 つの基本的な数学定数を加法、乗法、累乗演算で結びつけています。以下に SymPy を使ってこの公式を検証してみましょう。

読み込まれた記号の中で、E は自然指数の底、I は虚数単位、pi は円周率を表すため、上記の公式は次のように直接計算できます:

>>> E**(I*pi)+1
0

オイラーの等式は以下の公式で計算できます、
1668632127887.jpg
上記の公式を SymPy で証明するために、変数 x を導入する必要があります。SymPy では、数学記号は Symbol クラスのオブジェクトであるため、使用する前に作成する必要があります:

>>> x = Symbol('x')

expand 関数は公式を展開することができ、これを使って E**(I*x) を展開してみましょう:

>>> expand( E**(I*x) )
exp(I*x)

成功しませんでした、ただ別の書き方に変わっただけです。ここでの exp は math.exp や numpy.exp ではなく、sympy.exp であり、自然指数関数を表すためのクラスです。


math.exp(x)の指数を返します。x は数でなければなりません。
numpy.exp(x)の指数を返します。x は配列でも数でもかまいません。

>>> import math
>>> import numpy as np

>>> x = 2
>>> y = np.array([1,2,3,4])
>>> z = np.array([2])

>>> math.exp(x)
7.38905609893065
>>> math.exp(z)
7.38905609893065
>>> math.exp(y)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: only length-1 arrays can be converted to Python scalars

>>> np.exp(x)
7.3890560989306504
>>> np.exp(z)
7.3890560989306504
>>> np.exp(y)
array([  2.71828183,   7.3890561 ,  20.08553692,  54.59815003])

numpy は Python の配列演算パッケージです!


expand 関数にはキーワード引数 complex(虚部を持つパラメータタイプ、実質的には複素数、NumPy で説明済み)があり、これが True の場合、expand は公式を実数と虚数の 2 つの部分に分けます:

>>> expand(exp(I*x), complex=True)
I*exp(-im(x))*sin(re(x)) + cos(re(x))*exp(-im(x))

今回得られた結果は非常に複雑で、sin、cos、re、im はすべて sympy で定義されたクラスであり、re は実数部分を取得し、im は虚数部分を取得します。明らかにここでの演算は記号 x を複素数として扱っています。記号 x を実数として指定するために、次のように記号 x を再定義する必要があります:

>>> x = Symbol("x", real=True)
>>> expand(exp(I*x), complex=True)
I*sin(x) + cos(x)

Python におけるシングルクォート '' とダブルクォート "" の使い方:
これらには明確な違いや特別な使い方はなく、ほとんどの場合は互換性があります。しかし、注意が必要なケースがあります:

>>> str = "Hello,world!"
>>> str1 = 'Hello,world!'
>>> str2 = "Let's go!"
>>> str3 = 'Let\'s go!' # \はエスケープ文字で、システムにその後の'の役割が引用符ではないことを伝えます。
>>> str4 = "Say:\"Hello!\""
>>> str5 = 'Say:"Hello"'
>>> print str
Hello,world!
>>> print str1
Hello,world!
>>> print str2
Let's go!
>>> print str3
Let's go!
>>> print str4
Say:"Hello!"
>>> print str5
Say:"Hello"

以下のようにエスケープ文字なしで書くとエラーが発生します:

>>> str6 ='Let's go!'
  File "<stdin>", line 1
    str6 ='Let's go!'
               ^
SyntaxError: invalid syntax

2. 積分演算#

SymPy の記号積分関数 integrate は、記号的な積分を行うのに役立ちます。integrate は不定積分を行うことができます:

>>> integrate(x*sin(x), x)
-x*cos(x) + sin(x)

x の値の範囲を指定すると、integrate は定積分演算を行います:

>>> integrate(x*sin(x), (x, 0, 2*pi))
-2*pi

球体の体積を計算するために、まず円の面積を計算する方法を見てみましょう。円の半径を r と仮定すると、円上の任意の点の Y 座標関数は次のようになります:

>>> x, y, r = symbols('x,y,r')
>>> 2 * integrate(sqrt(r*r-x**2), (x, -r, r))
2*Integral((r**2 - x**2)**(1/2), (x, -r, r))

残念ながら、integrate 関数は結果を計算せず、入力した式をそのまま返しました。これは SymPy が r が 0 より大きいことを知らないためで、以下のように r を再定義すれば正しい答えが得られます:

>>> r = symbols('r', positive=True)
>>> circle_area = 2 * integrate(sqrt(r**2-x**2), (x, -r, r))
>>> circle_area
pi*r**2

次にこの面積公式に対して定積分を行うことで球体の体積を得ることができますが、X 軸の座標が変化するにつれて、対応する切断面の半径も変化します。今、X 軸の座標を x、球体の半径を r と仮定すると、x での切断面の半径は前の公式 y (x) を使って計算できます。
1668632183810.jpg
したがって、circle_area の変数 r を置き換える必要があります:

>>> circle_area = circle_area.subs(r, sqrt(r**2-x**2))
>>> circle_area
pi*(r**2 - x**2)

次に circle_area の変数 x を区間 - r から r まで定積分し、球体の体積公式を得ます:

>>> integrate(circle_area, (x, -r, r))
4*pi*r**3/3

subs を使った式の置換
subs 関数は式の中の記号を置き換えることができ、3 つの呼び出し方法があります:
expression.subs (x, y) : 式の中の x を y に置き換えます
expression.subs({x,u}) : 辞書を使って複数回置き換えます
expression.subs ([(x,y),(u,v)]) : リストを使って複数回置き換えます
注意してください、複数回の置き換えは順次実行されるため:
expression.subs([(x,y),(y,x)])
は 2 つの記号 x,y を交換することはできません。


ここでいくつかの例を挙げてみましょう。

>>> from sympy import *
>>> a,b,c = symbols('a,b,c') #小文字の最初の文字
>>> y = a + 10*b + 100*c
>>> y = y.subs(a,b)
>>> y
11*b + 100*c
>>> y = a + 10*b + 100*c
>>> y = y.subs([(a,b),(b,c)])
>>> y
111*c
>>> y = a + 10*b + 100*c
>>> y = y.subs([(a,b),(b,a)])
>>> y
11*a + 100*c
>>> x=Symbol('x') #この最初の文字は大文字に注意
>>> y = a + 10*b + 100*c
>>> y = y.subs([(a,x),(b,a),(x,b)])
>>> y
10*a + b + 100*c

最後にこの順序で手書きします。


完全な内容を見たい場合は、原文のアドレスにアクセスしてください: http://old.sebug.net/paper/books/scipydoc/sympy_intro.html

読み込み中...
文章は、創作者によって署名され、ブロックチェーンに安全に保存されています。