級数 ¶

モジュール sympy.series には級数およびその周辺の構成要素に関係する機能が実装されている。本稿ではそれらのうち、なんとか私の理解の範囲内にあるものについて見ていく。

Note

本文中のすべての IPython セッション中のサンプルコードで、以下のインポートおよび出力書式設定が済んでいるものとする。

init_printing(pretty_print=False)

極限 ¶

モジュール sympy.series.limits が提供する機能について記す。チュートリアルで見たこと以上のものはないようだ。

関数 limit(e, z, z0, dir='+')

変数 z が z0 に近づくときの関数 e(z) の極限値を求める。

引数 e は SymPy オブジェクト
引数 z は e の極限を考える変数。
引数 z0 は z を何に近づけるかを値で指定する。
キーワード引数 dir で右極限か左極限を指定する。プラスが右。

クラス Limit

関数版と同じことをするが、いわば遅延評価版。実際に極限を求める処理を実行するには、メソッド doit(**hint) を呼び出す。

サンプル ¶

よくある例を試す。

In [1]: limit(x * sin(1 / x), x, 0)
Out[1]: 0

In [2]: limit(sin(x), x, oo)
Out[2]: sin(oo)

In [3]: limit(gamma(x + 1 / 2)/gamma(x)/sqrt(x), x, oo)
Out[3]: 1

In [4]: limit(gamma(x) - 1 / (E**x - 1), x, 0)
Out[4]: -EulerGamma + 1/2

シンボルだけの数式の極限を計算する。

In [1]: limit(f(x), x, 7)
Out[1]: f(7)

In [2]: h = symbols('h')

In [3]: limit((f(x + h) - f(x)) / h, h, 0)
Out[3]: Subs(Derivative(f(_xi_1), _xi_1), (_xi_1,), (x,))

数列の極限も計算できる。

In [1]: limit((1 + x/n)**n, n, +oo)
Out[1]: exp(x)

二変数関数の極限を試す。近づけ方の指定がよくわからない。

In [1]: limit((x * y) / (x ** 2 + y ** 2), x, 0)
Out[1]: 0

In [2]: limit((x * y) / (x ** 2 + y ** 2), y, 0)
Out[2]: 0

In [3]: limit((x * y) / (x ** 2 + y ** 2), x, y)
Out[3]: 1/2

うまくいかない例を挙げる。これは SymPy の改良を期待できるだろうか。

In [1]: limit(fibonacci(k + 1)/fibonacci(k), k, oo)
Out[1]: Limit(fibonacci(k + 1)/fibonacci(k), k, oo, dir='-')

べき級数展開 ¶

モジュール sympy.series.series が提供する機能について記す。ここには関数 series がある。

関数 series(expr, x=None, x0=0, n=6, dir='+'): 関数 expr の点 x = x0 の周りでの n 次のべき級数展開を求める。各引数の意味は、これまで見てきた SymPy 関数での対応するものと同様。

サンプル ¶

だいたい想像通りに動作するようだ。

In [1]: series(exp(x))
Out[1]: 1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)

In [2]: series(exp(x), n=10)
Out[2]: 1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + x**6/720 + x**7/5040 + x**8/40320 + x**9/362880 + O(x**10)

In [3]: a = symbols('a')

In [4]: series(f(x), x, a, 3)
Out[4]: f(a) + (-a + x)*Subs(Derivative(f(_xi_1), _xi_1), (_xi_1,), (a,)) + (-a + x)**2*Subs(Derivative(f(_xi_1), _xi_1, _xi_1), (_xi_1,), (a,))/2 + O((-a + x)**3, (x, a))

In [5]: series(x**x, n=4)
Out[5]: 1 + x*log(x) + x**2*log(x)**2/2 + x**3*log(x)**3/6 + O(x**4*log(x)**4)

In [6]: series(gamma(x), n=2)
Out[6]: 1/x - EulerGamma + x*(EulerGamma**2/2 + pi**2/12) + O(x**2)

各出力の n 次以上の項にある O(x) については次節に記す。

ビッグオー記法 ¶

モジュール sympy.series.order が提供する機能について記す。

クラス Order

ビッグオー記法オブジェクトを表現するためのクラス。

変数がどの値に近づくときの極限の振る舞いの評価なのかを意識するとよい。アルゴリズムプログラミングの議論での記法は、常に「ある自然数が無限大に発散する」ときの計算量なり何なりのオーダーである。
リファレンスの、つまり docstring の説明がよくできている。一読するとよい。
デフォルトでは変数が 0 に近づくときの極限の評価を与える。無限大の評価をする場合は明示的に変数と oo とのペアを与える必要がある。
実はエイリアス O = Order が定義されている。当然これを積極的に利用する。
O(f(x), x) は内部で自動的に O(f(x).as_leading_term(x), x) と変形されるとのこと。これは一変数の場合で、それ以上の多変数の場合はこの規則がパラメーターのシンボル順に適用される。
ビッグオー記法オブジェクトは概念としては集合なので、同じ変数と関係する関数と Order オブジェクトに対しては演算子 in が適用できる。

サンプル ¶

何度も言うが、どこへの極限を考えているのかを常に意識するのがよい。

In [1]: 1 + x**10 + O(x**5)
Out[1]: 1 + O(x**5)

In [2]: O((x - 1)**4, (x, oo))
Out[2]: O(x**4, (x, oo))

In [3]: lg = lambda x: log(x, 2)

In [4]: O(9*lg(n) + 5 * (lg(n))**3 + 3*n**2 + 2*n**3, (n, oo))
Out[4]: O(n**3, (n, oo))

In [5]: from itertools import islice

In [6]: L = (O(i, (n, oo)) for i in (1, lg(n), n, n*lg(n), n**2, factorial(n)))

In [7]: all(i in j for i, j in zip(L, islice(L, 1, None)))
Out[7]: True

最初の例は、高次多項式を小さい次数に省略して表示するひとつの方法を示している。
上の例の lg の定義はこのデモにとっては必要ない。ビッグオー記法の定義により、対数に関するオーダーの比較は 2 でも 10 でも自然対数の底でも同じだ。

留数 ¶

モジュール sympy.series.residues が提供する機能について記す。

関数 residue(f, z, z0)

関数 f の点 z = z0 の周りで展開した Laurent 級数の -1 次の項の係数を返す。

本関数は留数定理で計算する前提の積分に対して大いに有用だ。留数を求めるには、いきなり級数が与えられているのならばともかく、ふつうは線積分や極限計算をすることで求める。そのような紙上での計算では手間がかかるので、これが使えるのはありがたい。

例 ¶

複素関数の入門書にある積分の演習問題から留数計算のみを抜き出して residue() を呼び出す。実際には極のうち、指定された領域に含まれるものを拾い出す行程が別途生じる。

>>> from sympy.abc import z
>>> residue(1/(z**2 + 1)**3, z, -I)
3*I/16
>>> residue(1/(z**2 + 1)**3, z, I)
-3*I/16

>>> residue(cos(z)/(1 + z**2), z, I)
-I*cosh(1)/2
>>> residue(cos(z)/(1 + z**2), z, -I)
I*cosh(1)/2

>>> [residue(1/(1 + z**6)), z, exp((2*k + 1) * pi*I/6) for k in range(6)]
[-(-1)**(1/6)/6, -I/6, -(-1)**(5/6)/6, (-1)**(1/6)/6, I/6, (-1)**(5/6)/6]

>>> residue(1/((z**2 - 1)*(z**3 - 1)), z, 1)
-1/4
>>> residue(1/((z**2 - 1)*(z**3 - 1)), z, -1)
1/4
>>> solve(z**2 + z + 1, z)
[-1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2]
>>> alpha, beta = _[0], _[1]
>>> residue(1/((z**2 - 1)*(z**3 - 1)), z, alpha)
sqrt(3)*I/9
>>> residue(1/((z**2 - 1)*(z**3 - 1)), z, beta)
-sqrt(3)*I/9

\((-1)^{1/6}\) などが目立つが、複素数では指数関数は多価なのでこう表現されるのだろう。

数列 ¶

モジュール sympy.series.sequences が提供する機能について記す。ここには一連の数列表現・操作クラス群と、その生成関数の定義がある。

関数 sequence(seq, limits=None)

数列を生成する。後述のクラス SeqExpr のサブクラスのオブジェクトを返す。

戻り値の型は引数 seq の型で決まるらしい。リファレンスによると次のような指定方法がある。
- 数列の各項を tuple オブジェクトで指定する。このときは周期的な数列が得られる。例 (1, 2, 3)
- 一般項をシンボルの式で表現する。例 n ** 2
キーワード引数 limits は数列の定義域を指定するために使う。
- デフォルトは区間 [0, oo) 上の添字が許される数列を生成する。
- 有限数列でも無限数列でも構わないが、整数全体を区間として指定することはできない。例外 ValueError が送出する。

クラス図 ¶

このモジュールに定義されているクラスを継承関係だけ図示するとこういう感じになる：

クラス EmptySequence は空の数列を表現するためのもの。オブジェクトを S.EmptySequence として参照できる。

クラス `SeqBase`¶

クラス SeqBase はすべての数列クラスのスーパークラス。次のようなプロパティーが用意されている：

gen: 数列の生成要素。状況によって肩が異なる。
interval, start, stop, length: 数列の定義域に関する情報。
variables: 数列の添字シンボル。一般に多重数列であることもあるので tuple 型である。

また、数列オブジェクトに対して次のような演算子が使える：

単項演算子 -
二項演算子 +=, -=, *=, +, -, *
- なぜかスカラー倍がサポートされていない。
角括弧 a[i]

数列が有限であるときに限り、このオブジェクトは iterable である。

メソッドは次のようなものがある：

coeff(pt)

数列の第 pt 項を求める。

定義域内の任意の実数で評価してくれるようだ。

coeff_mul(other):

これが数列のスカラー倍を実現する。引数 other にスカラーを指定するとよい。

クラス `SeqExpr` とそのサブクラス ¶

クラス SeqExpr は数列オブジェクトを表現するための基礎部分を提供する。スーパークラスで宣言されているプロパティーを実装するためだけにある。

クラス SeqFormula は数列を一般項で表現する。スーパークラスのものに加えて、次のプロパティーが利用できる。

formula: 普通は一般項を表現する式を参照する。実装は単に self.gen を参照する。

クラス SeqPer は周期的な数列を表現する。スーパークラスのものに加えて、次のプロパティーが利用できる。

period: 数列の周期を返す。実装は単に len(self.gen) を参照する。
periodical: 数列の周期部分を返す。実装は単に self.gen を参照する。

クラス `SeqExprOp` とそのサブクラス ¶

このクラスは数列オブジェクトを（オブジェクト指向プログラミングの意味で）合成するための素となる。サブクラス SeqAdd と SeqMul は数列オブジェクト同士の和と積をそれぞれ表現する。

ユーザーは明示的にこれらのクラスからオブジェクトを生成することはなさそうだ。

Fourier 級数 ¶

モジュール sympy.series.fourier が提供する機能について記す。それはもちろん Fourier 級数を計算するためのものだ。

関数 fourier_series(f, limits=None)

関数 f の Fourier 級数を求める。戻り値の型は次に記すクラス FourierSeries のオブジェクトだ。

クラス FourierSeries

関数版と同じことをするが、例によって遅延評価版。次のメソッドがある：

scale(s)

f(x) から s * f(x) の Fourier 級数を返す。

scalex(s)

同様に f(s * x) を返す。

shift(s)

同様に f(x) + s を返す。

shiftx(s)

同様に f(x + s) を返す。

trancate(n=3)

Fourier 級数の最初の n 項までを返す。

n=None の場合は特別に、初項から順に項を yield する。

ただし、

与えられた関数を f とした。
最初の引数 self は省略した。
s は変数 x と独立な項とする。

サンプル ¶

一般に Fourier 級数を計算するのは時間がかかるので、利用する側が計算量を減らす工夫をするのが腕の見せどころとなる。上述の s を引数に取るメソッドが応用できる場合は積極的にそうするのだ。

In [1]: fourier_series(t/2).truncate(3)
Out[1]: sin(t) - sin(2*t)/2 + sin(3*t)/3

In [2]: fourier_series(t/2, (t, -3 * pi, 3 * pi)).truncate(3)
Out[2]: 3*sin(t/3) - 3*sin(2*t/3)/2 + sin(t)

In [3]: fourier_series(abs(t)).truncate(2)
Out[3]: cos(t)*Integral(cos(t)*Abs(t), (t, -pi, pi))/pi + Integral(Abs(t), (t, -pi, pi))/(2*pi)

In [4]: fourier_series(exp(I * t))
Out[4]: FourierSeries(exp(I*t), (t, -pi, pi), (0, SeqFormula(Piecewise((pi, Or(Eq(_n, -1), Eq(_n, 1))), (-2*_n*sin(_n*pi)/(_n**2 - 1), True))*cos(_n*t)/pi, (_n, 1, oo)), SeqFormula(Piecewise((-I*pi, Eq(_n, -1)), (I*pi, Eq(_n, 1)), (-2*I*sin(_n*pi)/(_n**2 - 1), True))*sin(_n*t)/pi, (_n, 1, oo))))

In [5]: fourier_series(exp(-t)).truncate(3)
Out[5]: (-exp(pi)/2 + exp(-pi)/2)*sin(t)/pi + (-2*exp(-pi)/5 + 2*exp(pi)/5)*sin(2*t)/pi + (-exp(pi)/2 + exp(-pi)/2)*cos(t)/pi + (-exp(-pi)/5 + exp(pi)/5)*cos(2*t)/pi + (-exp(-pi) + exp(pi))/(2*pi)

In [6]: fourier_series(exp(-t**2)).truncate(2)

[1][2] 単純な一次関数を異なる区間指定で Fourier 級数を計算させた。結果も異なる。

下の図を参照。カーブが大きいほうが区間指定なしのほうだ。
[3] 絶対値の二次の Fourier 級数を計算した。定積分が含まれているままだが、ここは簡単にしてくれないのか。
[4] これは何かおかしい。なぜこのように長い式が出力されるのだろうか。
[5] 指数関数の三次の Fourier 級数を計算した。これは正しいだろう。
[6] 簡単な Gauss 関数の Fourier 級数は、しかし私の環境では関数呼び出しから戻ってこなかった。