Pade 近似の係数を求めるプログラムを C で書いてみた

Pade 近似とは

関数 $f(x)$ を有理関数 $[m, n](x):=P_m(x)/Q_n(x)$ *1で近似する方法。
$m+n$ 次までのマクロリーン展開の係数が両者で一致することが条件。
詳しくは参考 URL まで。

アルゴリズム

$f(x)$ の $i$ 次のマクロリーン展開係数を $a_i$ と書く。
また、 $P_m(x)$ と $Q_n(x)$ の多項式の係数をそれぞれ $p_i$ 、 $q_i$ とする。*2
Pade 近似の条件は
$p_k-\sum_{j=1}^{k}a_{k-j}q_j=a_k~~~(k=1,2,\cdots,m+n)$
となります。
なお、有理式の形なので $q_0=1$ とおいてよく、それにより $p_0=a_0$ も導かれます。

条件は $m+n$ 元連立方程式 $B\mathbf{x}=\mathbf{c}$ を解くことと同値です。
ここで、 $\mathbf{x}={}^{\mathrm T}(p_1,p_2,\cdots,p_m,q_1,q_2,\cdots,q_n)$ 、 $\mathbf{c}={}^{\mathrm T}(a_1,a_2,\cdots,a_{m+n})$ 、および
$B=\left(\begin{array}{cccccccc} 1 & & & & & -a_0 & & \\ & \ddots & & & & \vdots & \ddots & \\ & & \ddots & & & -a_{n-1} & \cdots & -a_0 \\ & & & \ddots & & -a_{n} & \cdots & -a_1 \\ & & & & 1 & & & \\ & & & & & \vdots & \ddots & \vdots \\ & & & & & & & \\ & & & & & -a_{m+n-1} & \cdots & -a_m \end{array}\right)$
です。あとは拡大係数行列を定義して掃き出し法使ってもいいし、 $B$ の逆行列求めてしまってもいい。*3

実装

C で書きました。行列計算の部分は部分 Pivot 法を用いた掃き出し法を使いました。
実用上、係数読み込みの部分が若干長くなってます。Pade 近似の部分は後半。

コンパイルするときは、

$ gcc pade.c -o pade

とすれば OK。
微係数 or 展開係数をファイルから読み込むときはそれをコマンドライン引数に加える。
ファイル形式は、係数が , で区切られて一行で書かれているもの。
空行と行先頭が # なものは無視するようにしています。

外部ファイルに出力したければ標準出力で。

例

# input.dat
# f(x) = e^x
# Differential Coefficients

1.0, 1.0, 1.0, 1.0, 1.0

というファイルを作って、

$ ./pade input.dat > output.dat

とすればいい。

[補足] Global Pade 近似

参考 URL に Global Pade 近似なるものを発見。 $x$ が小さいところは普通の Pade 近似で合わせて、大きいところは漸近形に従うように定義するもの。

格子比熱に関するデバイの式や電気抵抗率に関するグリューナイゼンの式と言った積分形式で書かれたものを Pade 近似して Fitting を行うときにも用いられている。詳しくは以下の論文を見て下さい。

Structural, Thermal, Magnetic and Electronic Transport Properties of the LaNi₂(Ge_1-xP_x)₂ System
R. J. Goetsch et al. : Phys. Rev. B 85, 054517 (2012)., arXiv:1112.1864

なお、ねがてぃぶろぐ Scilabで関数フィッティング: 金属の電気抵抗によると積分の形式のまま再帰的に Fitting 出来るらしい。Scilab を使っているらしいがどういった計算が中でされているかは分かりません。*4

参考になるもの

http://next1.msi.sk.shibaura-it.ac.jp/MULTIMEDIA/numeanal2/node19.html
Pade 近似の詳しい説明があります。ただ、「Pade の有理関数近似アルゴリズム」には以下の間違いがあります。

(4) の変数 $j$ は $i$
(6) の代入式は $b_{ij}\leftarrow 1$

まあ正直、一旦全部の要素を 0 にしてしまって、上のような行列を作るという指針（私の pade.c はこういう指針。）のほうがいいと思います。for ループは確かに増えますが正直大きいループではないし。
あと、このサイトの行列 $(b_{i,j})_{1\leq i\leq N, 1\leq j\leq N+1}$ は拡大係数行列です。