テイラー展開、パデ近似をpythonで簡単に表せる日が来た

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

f:id:koideforest:20171019204501p:plain

ここで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()