mpmathというモジュールの中にテイラー展開をしてくれmpmath.taylorがいらっしゃる。
Function approximation — mpmath 1.0.0 documentation
mpmath.taylorは展開係数の入ったリストを返してくれる。リストの中身は昇べきの順で並んでいる。
一方、mpmath.polyvalは係数のリストを引数にして多項式の値を返してくれる。
Polynomials — SymPy v0.7.0 documentation
ただし、リストの順番は次数に対して降べきの順である。
したがって、list[::-1]でリストの順番を逆にすればtaylor展開の値がpolyvalですぐに求まる。
#!/usr/bin/env python import numpy as np import mpmath as mp import matplotlib.pyplot as plt def taylor( f, x, x0, n ): coefficients = mp.taylor( f, x0, n) return mp.polyval( coefficients[ ::-1], x - x0 )
主要部分はたったこれだけである。
例えばsin関数のマクローリン展開を20乗の項まで展開すると(奇数乗の展開係数はゼロであるため、偶数の次数だけ入力すれば違いを見るのには十分)、
x = np.linspace( 0., 10., 100 ) f = mp.sin x0 = 0. for i in range( 10 ): plt.plot( x, [ taylor( f, xi, x0, 2 * i ) for xi in x ] ) plt.xlim( 0., 10. ) plt.ylim( -5., 5. ) plt.plot( x, [ f( xi ) for xi in x ] ) plt.show()
ここでsin関数はmpmathのものでないと動かない。numpy.sinは少なくとも動かなかった。
更に、パデ近似(ペデという人も居るかも)も似たノリで簡単に出来る。
Function approximation — mpmath 1.0.0 documentation
Taylor展開よりも良いと聞かされて来たが、展開係数を求めるのが億劫で、なかなかどんなもんなのか感覚を掴めずにいた。
def pade( f, x, x0, nt, np, nq ) c_taylor = mp.taylor( f, x0, nt ) c_p, c_q = mp.pade( c_taylor, np, nq ) p = mp.polyval( c_p[::-1], x - x0 ) q = mp.polyval( c_q[::-1], x - x0 ) return p / q
パデ近似はTaylor展開の次数に合わせて近似を行い、np + nq = ntの制約がかかっている。
とりあえず適当にsin関数を近似してみると、
plt.plot( x, [ taylor( f, xi, x0, 10 ) ], label="taylor" ) plt.plot( x, [ pade( f, xi, x0, 10, 5, 5 ) for xi in x ], label="pade(5,5)" ) plt.plot( x, [ pade( f, xi, x0, 10, 3, 7 ) for xi in x ], label="pade(3,7)" ) plt.plot( x, [ pade( f, xi, x0, 10, 7, 3 ) for xi in x ], label="pade(7,3)" ) plt.plor( x, [ f( xi ) for xi in x ], label="exact" ) plt.legend() plt.show()
sin関数だとあまり大きな改善は見られない。
少なくとも、これでいちいち手打ちしなくても近似形が得られるのでだいぶ心がスッキリ。