コンパクトスキームはなぜPadé（パデ）スキームと呼ばれるか

概要

微分を数値的に計算するコンパクトスキームは，Padé（パデ）スキームと呼ばれることもある．なぜPadéスキームと呼ばれるのかをまとめた．

差分法による微分の計算

数値シミュレーションにおいて微分値が解析的に求められることは滅多にないため，計算機を用いて近似的に計算することになる．このような手続きを，数値的に微分を計算するという．
$x$の関数$f(x)$がなめらかであれば，$f$の微分は次式で定義される．

f'=\frac{\partial f}{\partial x}=\lim_{\Delta x \rightarrow 0} \frac{f(x+\Delta x)-f(x)}{\Delta x}

関数$f$に対して，$x$方向に$\Delta x$だけ離れた2点を考え，その2点の間隔を無限小に近づけた時の極限が$f$の勾配である．しかし，この無限小というのが曲者で，計算機では表現できない．そのため何らかの方法で近似的に計算する必要がある．

このとき，$\Delta x$に対する制約を緩め，無限小ではなく関数$f$の変化に対して十分小さいと見なすと，微分の定義式は減算と除算を組み合わせた差分で表現できるようになる．

f'=\frac{\partial f}{\partial x}\approx \frac{f(x+\Delta x)-f(x)}{\Delta x}

記述の簡略化のために，座標値$x$を離散点の番号$i$で表し，その地点での関数値$f(x)$を$f_i$と書く．$x\pm n\Delta x$は$i\pm n$と表される．

点の取り方にはいくつもの選択肢がある．代表的な式として，2次，4次，6次精度の中心差分式を示す．

\begin{eqnarray}
f'_i&\approx&\frac{f_{i+1}-f_{i-1}}{2\Delta x}\\
&\approx&\frac{-f_{i+2}+8f_{i+1}-8f_{i-1}+f_{i-2}}{12\Delta x}\\
&\approx&\frac{f_{i+3}-9f_{i+2}+45f_{i+1}-45f_{i-1}+9f_{i-2}+f_{i-3}}{60\Delta x}
\end{eqnarray}

1点の微分値を計算するのに要する離散点やその幾何学的配置をステンシルという．下図からわかるように，この差分の計算法では，$f$の値がわかれば，その微分値が直ちに計算できる．このような方法を陽的差分（あるいは陽的スキーム）とよぶ．2次，4次，6次精度の陽的中心差分は，それぞれ3点，5点，7点ステンシルである．精度を上げる度にステンシルの幅が広くなっていく．また，$n$次精度を達成するためには，どうやら$n+1$点を使う必要がある．

コンパクトスキーム

陽的な差分では，$n$次精度を達成するために$n+1$点の値を使う必要があった．陽的な差分法とは異なる方法で，差分を計算する方法（差分スキーム）として，本稿の主題であるコンパクトスキーム（コンパクト差分）がある．コンパクトスキームは，近傍点の差分値も未知数として差分近似式を構成する．コンパクトスキームは古くから興味の対象となり，色々と調べられてきた．Lele¹が統一的なまとめ方をしたので，Leleのコンパクト差分とも呼ばれることもあるが，それは正確ではない．実際，Lele¹の論文の中では

ここで示すスキームは，一般化したPadéスキームである．（著者訳）

と述べている．

Leleが一般化した，変数$f$の1階微分$f'$に対する7点ステンシル，5パラメータのコンパクト差分は次式で表される¹．

\beta f'_{i-2}+\alpha f'_{i-1}+f'_i+\alpha f'_{i+1}+\beta f'_{i+2}=a\frac{f_{i+1}-f_{i-1}}{2\Delta x}+b\frac{f_{i+2}-f_{i-2}}{4\Delta x}+c\frac{f_{i+3}-f_{i-3}}{6\Delta x}

ここで，

\alpha=\frac{1}{4},\beta=0,a=\frac{6}{4},b=0,c=0

を採用すると，ステンシルが3点のコンパクト差分を作ることができる．

話がずれてきているので，本稿の目的を再確認すると，コンパクトスキームはなぜPadéスキームと呼ばれるかを説明する事であった．Leleの論文には，Padéスキームについて何カ所かで言及されている．そのうちの一つを見ると，

$\alpha$が0になると，よく知られた4次精度の中心差分になる. 同様に，$\alpha = 1/4$で古典的Padéスキームが得られる．さらに，$\alpha = 1/3$では誤差の主要項が消え，このスキームは形式的に6次精度になる．
中略
$\alpha = 1/4$および$\alpha = 1/3$として得られる3重対角スキームは，Collatz²によって与えられた．（著者訳）

とある．本稿の対象は，$\alpha = 1/4$としたときに得られる古典的Padéスキームである．古典的Padéスキームは，Kopal³によってCollatzよりも早い1961年に与えられている．Kopalによる導出過程を追いながら，どこでPadé近似を用いるか確認しよう．

古典的Padéスキームの導出

下準備として，物理量の表し方を整理しておく．基準となる位置$x_0$から，$\Delta x$間隔で等間隔に離散点$x_i= x_0 + i\Delta x$を設け，その地点で関数値$f(x_i)$を定義する．その地点以外では，関数値は定義されず，点と点の繋がりを示す情報もない．$f(x_i)$は簡単のために$f_i$と書く．

\begin{eqnarray}
x_i &=& x_0 + i\Delta x\\
f_i &=& f(x_i)
\end{eqnarray}

演算子の定義

演算子$\mathcal{F}$を次のように定義する．

\mathcal{F} = e^{\Delta x \frac{d}{dx}}

この演算子をTaylor展開すると

e^{\Delta x \frac{d}{dx}} = \mathcal{I}+\Delta x\frac{d}{dx}+\frac{\Delta x^2}{2!}\frac{d^2}{dx^2}+\frac{\Delta x^3}{3!}\frac{d^3}{dx^3}+\cdots

である．$\mathcal{I}$は自分自身を得る演算子である．

\mathcal{I}f_i = f_i

実際に関数に適用すると

\begin{eqnarray}
\mathcal{F}f_i &=& f_i+\Delta x\frac{d f_i}{dx}+\frac{\Delta x^2}{2!}\frac{d^2 f_i}{dx^2}+\frac{\Delta x^3}{3!}\frac{d^3 f_i}{dx^3}+\cdots\\
&=& f_{i+1}.
\end{eqnarray}

したがって，演算子$\mathcal{F}$は$\Delta x$ずれた点の値を返す前進演算子である．

$f_{i-1}$は

\begin{eqnarray}
\mathcal{F}^{-1}f_i&=&e^{-\Delta x \frac{d}{dx}}f_i\\
&=& f_{i-1},
\end{eqnarray}

$f_{i+\frac{1}{2}}$は

\begin{eqnarray}
\mathcal{F}^\frac{1}{2}f_i&=&e^{\frac{\Delta x}{2} \frac{d}{dx}}f_i\\
&=&f_{i+\frac{1}{2}}
\end{eqnarray}

として求める．

次に，$i$の近傍で$f$の平均を取る演算子$\mathcal{A}$と差を取る演算子$\mathcal{S}$を定義する．

\begin{eqnarray}
\mathcal{A}f_i &=& \frac{1}{2}(f_{i+\frac{1}{2}}+f_{i-\frac{1}{2}})\\
\mathcal{S}f_i &=& f_{i+\frac{1}{2}}-f_{i-\frac{1}{2}}
\end{eqnarray}

ここで，

\begin{eqnarray}
\mathcal{A}^2f_i &=& \frac{1}{4}(f_{i+1}+2f_{i}+f_{i-1})\\
\mathcal{S}^2f_i &=& f_{i+1}-2f_{i}+f_{i-1}
\end{eqnarray}

であるから，$\mathcal{A}$と$\mathcal{S}$には

\mathcal{A} = \sqrt{1+\frac{1}{4}\mathcal{S}^2}

なる関係があることが示される．

今，離散化された点上での微分を考えているので，$i\pm \frac{1}{2}$の点では値が定義されていない．近傍点から補間をすれば$i\pm \frac{1}{2}$での値$f_{i\pm\frac{1}{2}}$も求められるだろうが，それでは微分の精度が補間の精度に制限されてしまう．

そこで，演算子を組み合わせて$i$上での$f$の値を使いつつ，精度を制御できる式を構築する．

微分演算子の導出

Kopal³の手順に沿って，微分演算子を導出する．先に定義した前進演算子$\mathcal{F}$を組み合わせ，$d/dx$を表す式を得る．

\frac{1}{2}\mathcal{S}f_i = \frac{1}{2}(f_{i+\frac{1}{2}}-f_{i-\frac{1}{2}})\\
= \frac{1}{2}\left[ e^{\frac{\Delta x}{2} \frac{d}{dx}}-e^{-\frac{\Delta x}{2} \frac{d}{dx}}\right]f_i

変形には演算子$\mathcal{F}$が使われている．ここで，

\sinh x= \frac{e^x-e^{-x}}{2}

を用いて変形すると

\frac{1}{2}\mathcal{S}f_i = \left[\sinh\left(\frac{\Delta x}{2}\frac{d}{dx}\right)\right]f_i

が得られるので，演算子$\mathcal{S}$は次のように書かれる．

\frac{1}{2}\mathcal{S} = \left[\sinh\left(\frac{\Delta x}{2}\frac{d}{dx}\right)\right]

逆三角関数を用いて微分演算子を求めると

\frac{d}{dx}=\frac{2}{\Delta x}\sinh^{-1}\left(\frac{\mathcal{S}}{2}\right)

となる．

微分演算子の近似

$\sinh^{-1}x$のTaylor展開は

\sinh^{-1}x = \sum_{n=0}^{\infty}{\frac{(-1)^n(2n)!}{2^{2n}(n!)^2}\frac{x^{2n+1}}{2n+1}}

と書かれるので，適当な項まで展開すると⁴

\sinh^{-1}x = x-\frac{1}{2}\frac{x^3}{3}+\frac{3}{8}\frac{x^5}{5}-\frac{15}{48}\frac{x^7}{7}+\cdots.

$x$に$\mathcal{S}/2$を代入すると，近似された微分演算子が得られる．

 \frac{d}{dx}=\frac{1}{\Delta x}\left( \mathcal{S}-\frac{1}{3}\frac{\mathcal{S}^3}{2^3}+\frac{3}{20}\frac{\mathcal{S}^5}{2^5}-\frac{15}{168}\frac{\mathcal{S}^7}{2^7}+\cdots\right).

$\mathcal{S}$の奇数乗は，$i$ではなく$i\pm \frac{1}{2}$上の点を用いた式になるので，Kopalは$i$における値を使うように上式を修正した³．

 \frac{1}{\mathcal{A}\mathcal{S}}\frac{d}{dx}=\frac{1}{\mathcal{A}\mathcal{S}}\frac{1}{\Delta x}\left( \mathcal{S}-\frac{1}{3}\frac{\mathcal{S}^3}{2^3}+\frac{3}{20}\frac{\mathcal{S}^5}{2^5}-\frac{15}{168}\frac{\mathcal{S}^7}{2^7}+\cdots\right).

ここで，

\mathcal{A} = \sqrt{1+\frac{1}{4}\mathcal{S}^2}

なる関係があったことを思い出すと，

\begin{eqnarray}
 \frac{1}{\mathcal{A}\mathcal{S}}\frac{d}{dx}&=&\frac{1}{\mathcal{A}\mathcal{S}}\frac{1}{\Delta x}\sinh^{-1}\left(\frac{\mathcal{S}}{2}\right)\\
 &=&\frac{1}{\mathcal{S}}\frac{1}{\Delta x}\sinh^{-1}\left(\frac{\mathcal{S}}{2}\right)\left(1+\frac{1}{4}\mathcal{S}^2\right)^{-1/2}
\end{eqnarray}

$\mathcal{A}^{-1}$のTaylor展開は

 \mathcal{A}^{-1}=\left(1+\frac{1}{4}\mathcal{S}^2\right)^{-1/2}=1-\frac{1}{8}\mathcal{S}^2+\frac{3}{128}\mathcal{S}^4-\frac{5}{1024}\mathcal{S}^6+\frac{35}{32768}\mathcal{S}^8+\cdots

なので，修正された微分演算子を得るには

\frac{1}{\mathcal{S}}\frac{1}{\Delta x}\left( \mathcal{S}-\frac{1}{3}\frac{\mathcal{S}^3}{2^3}+\frac{3}{20}\frac{\mathcal{S}^5}{2^5}-\frac{15}{168}\frac{\mathcal{S}^7}{2^7}+\cdots\right)\left(1-\frac{1}{8}\mathcal{S}^2+\frac{3}{128}\mathcal{S}^4-\frac{5}{1024}\mathcal{S}^6+\frac{35}{32768}\mathcal{S}^8+\cdots\right)

を計算することになる．A4用紙1/2ページほどの計算をするか，あるいはWolfram Alphaにお願いすると直ちに答えが得られる．

 \frac{1}{\mathcal{A}\mathcal{S}}\frac{d}{dx}=\frac{1}{\Delta x}\left(1-\frac{1}{6}\mathcal{S}^2+\frac{1}{30}\mathcal{S}^4-\frac{1}{140}\mathcal{S}^6+\frac{1}{630}\mathcal{S}^8+\cdots\right)

$d/dx=$の形に変形すると，陽的差分の公式が求められる．演算子を$f_i$に作用させ，

 \frac{d}{dx}f_i=\frac{\mathcal{A}\mathcal{S}}{\Delta x}\left(1-\frac{1}{6}\mathcal{S}^2+\frac{1}{30}\mathcal{S}^4-\frac{1}{140}\mathcal{S}^6+\frac{1}{630}\mathcal{S}^8+\cdots\right)f_i

括弧内を第1項目まで採用すれば，

\begin{eqnarray}
\mathcal{A}\mathcal{S}f_i&=&\frac{1}{2}\left[\left(\mathcal{S}f\right)_{i+\frac{1}{2}}+\left(\mathcal{S}f\right)_{i-\frac{1}{2}}\right]\\
&=&\frac{1}{2}\left[(f_{i+1}-f_{i})+(f_{i}-f_{i-1})\right]\\
&=&\frac{1}{2}(f_{i+1}-f_{i-1})
\end{eqnarray}

なので

\frac{d}{dx}f_i=\frac{f_{i+1}-f_{i-1}}{2\Delta x}

となって，これは2次精度の中心差分と一致する．また，括弧内を第2項目まで採用すれば，

\mathcal{A}\mathcal{S}f_i-\frac{1}{6}\mathcal{A}\mathcal{S}^3f_i=\frac{1}{2}(f_{i+1}-f_{i-1})-\frac{1}{12}(f_{i+2}-2f_{i+1}+2f_{i-1}-f_{i-2})

となり，項をまとめて4次精度の中心差分

\frac{d}{dx}f_i=\frac{-f_{i+2}+8f_{i+1}-8f_{i-1}+f_{i-2}}{12\Delta x}

が得られる．

高次の項まで採用すれば高次精度の中心差分式が得られることは明かであるが，それに伴ってステンシル幅が広くなる．この問題に対して，狭いステンシル幅で高次精度を達成するために，多項式

 1-\frac{1}{6}\mathcal{S}^2+\frac{1}{30}\mathcal{S}^4-\frac{1}{140}\mathcal{S}^6+\frac{1}{630}\mathcal{S}^8+\cdots

をPadé近似するのである．

Padé近似

物理の授業よろしく，ここで一度横道に逸れる．多項式のPadé近似がどのような手法かを理解するために，Padé近似を簡単に説明する．

ざっくり言うと，Padé近似は関数を次のような多項式で近似する手法である．

r(x) \approx R(x)= \frac{\sum_{i=0}^{m}a_ix^i}{1+\sum_{i=0}^{n}b_ix^i}

$m\ge 0,n\ge 1$である．$r(x)$のTaylor展開が$m+n$次の項まで$R(x)$と一致するとき，$R(x)$は$r(x)$のPadé近似と呼ばれる．つまり，Padé近似を用いれば，$m+n$次のTaylor展開と同じ精度を得るのに，分子が$m$次多項式，分母が$n$次多項式となる．本稿でPadé近似したい多項式では，$\mathcal{S}$の次数がステンシル幅に影響するので，低い次数の多項式の組合せで式を置き換えられるのであれば，狭いステンシル幅で高次精度を達成できるようになるのでとても好都合である．

例として，$r(x)=e^x$のPadé近似を求めてみよう．$m=2, n=3$を選ぶとき，$r$のTaylor展開と5次の項まで一致すればよいので，

1+x+\frac{x^2}{2}+\frac{x^3}{6}+\frac{x^4}{24}+\frac{x^5}{120}= \frac{a_0+a_1x+a_2x^2}{1+b_1x+b_2x^2+b_3x^3}

となる$a_i,b_i$を求めれば，Padé近似できたということになる．

\left(1+x+\frac{x^2}{2}+\frac{x^3}{6}+\frac{x^4}{24}+\frac{x^5}{120}\right)\left(1+b_1x+b_2x^2+b_3x^3\right) =\left(a_0+a_1x+a_2x^2\right)

の両辺に対して$x$の次数で係数を揃えると，以下のような関係が得られる．

$x$の次数	式
0	$a_0=1$
1	$a_1=1+b_1$
2	$a_2=\frac{1}{2}+b_1+b_2$
3	$\frac{1}{6}+\frac{b_1}{2}+b_2+b_3=0$
4	$\frac{1}{24}+\frac{b_1}{6}+\frac{b_2}{2}+b_3=0$
5	$\frac{1}{120}+\frac{b_1}{24}+\frac{b_2}{6}+\frac{b_3}{2}=0$

これらを基に連立方程式

\begin{pmatrix}
    1&0&0& 0& 0& 0\\
    0&1&0&-1& 0& 0\\
    0&0&2&-2&-2& 0\\
    0&0&0& 3& 6& 6\\
    0&0&0& 4&12&24\\
    0&0&0& 5&20&60\\
\end{pmatrix}
\begin{pmatrix}
    a_0\\
    a_1\\
    a_2\\
    b_1\\
    b_2\\
    b_3\\
\end{pmatrix}
=
\begin{pmatrix}
    1\\
    1\\
    1\\
    -1\\
    -1\\
    -1\\
\end{pmatrix}

を作り，解くことで

\begin{pmatrix}
    a_0\\
    a_1\\
    a_2\\
    b_1\\
    b_2\\
    b_3\\
\end{pmatrix}
=
\begin{pmatrix}
    1\\
    \frac{2}{5}\\
    \frac{1}{20}\\
    -\frac{3}{5}\\
    \frac{3}{20}\\
    -\frac{1}{60}\\
\end{pmatrix}

を得る．これらより，$e^x$の$2+3$次のPadé近似

e^x\approx \frac{1+a_1x+\frac{2}{5}x^2}{1-\frac{3}{5}x+\frac{3}{20}x^2-\frac{1}{60}x^3}+O(x^6)

が得られる．

微分演算子のPadé近似

Padéスキームの話に戻ろう．幸いなことに，本稿でPadé近似したい関数$r(\mathcal{S})$（正確には演算子だが）はTaylor展開の形で求められているので，これを任意の$m,n$を用いてPadé近似すればよい．つまり，

1-\frac{1}{6}\mathcal{S}^2+\frac{1}{30}\mathcal{S}^4-\frac{1}{140}\mathcal{S}^6+\frac{1}{630}\mathcal{S}^8+\cdots=\frac{\sum_{i=0}^{m}a_i\mathcal{S}^i}{1+\sum_{i=0}^{n}b_i\mathcal{S}^i}.

$\mathcal{S}$は偶数次の項しかないので，最も項が少なくなるように$m=0,n=2$を選べば次式が得られる．

1-\frac{1}{6}\mathcal{S}^2=\frac{1}{1+b_2\mathcal{S}^2}.

この式は

\left(1-\frac{1}{6}\mathcal{S}^2\right)\left(1+b_2\mathcal{S}^2\right)=1

であるから，$\mathcal{S}^2$の係数は

-\frac{1}{6}+b_2=0

より$b_2=1/6$が定まる．前置きは長かったが，これで終わりである．

Padé近似された微分演算子

 \frac{d}{dx}=\frac{\mathcal{A}\mathcal{S}}{\Delta x}\left(\frac{1}{1+\frac{1}{6}\mathcal{S}^2}\right)

を$f_i$に作用させると

\begin{eqnarray}
 \frac{d}{dx}f_i&=&\frac{\mathcal{A}\mathcal{S}}{\Delta x}\left(\frac{1}{1+\frac{1}{6}\mathcal{S}^2}\right)f_i\\
&=&\left(\frac{1}{1+\frac{1}{6}\mathcal{S}^2}\right)\frac{f_{i+1}-f_{i-1}}{2\Delta x}
\end{eqnarray}

が得られる．$df_i/dx$を$f'_i$と書き，右辺の演算子を左辺に持って行く．

\left(1+\frac{1}{6}\mathcal{S}^2\right)f'_i=\frac{f_{i+1}-f_{i-1}}{2\Delta x}.

左辺の演算子を展開し，

\begin{eqnarray}
\left(1+\frac{1}{6}\mathcal{S}^2\right)f'_i&=&f'_i+\frac{1}{6}\left(f'_{i-1}-2f'_i+f'_i\right)\\
&=&\frac{1}{6}f'_{i-1}+\frac{2}{3}f'_i+\frac{1}{6}f'_{i+1}
\end{eqnarray}

$f'_i$の係数が$1$となるように変形すると，古典的Padéスキーム

\frac{1}{4}f'_{i-1}+f'_i+\frac{1}{4}f'_{i+1}=\frac{3}{2}\frac{f_{i+1}-f_{i-1}}{2\Delta x}

が得られた．

高次近似

当然，$m=2,n=2$を選べばより高次精度のPadéスキームが得られる．

1-\frac{1}{6}\mathcal{S}^2+\frac{1}{30}\mathcal{S}^4=\frac{a_0+a_2\mathcal{S}^2}{1+b_2\mathcal{S}^2}.

この式は

\left(1-\frac{1}{6}\mathcal{S}^2+\frac{1}{30}\mathcal{S}^4\right)\left(1+b_2\mathcal{S}^2\right)=a_0+a_2\mathcal{S}^2

であるから，$\mathcal{S}$の次数で係数を整理すると

\begin{eqnarray}
1-a_0&=&0\\
b_2-\frac{1}{6}-a_2&=&0\\
-\frac{1}{6}b_2+\frac{1}{30}&=&0
\end{eqnarray}

より$a_0=1, a_2=1/30, b_2=1/5$が定まる．

微分演算子

 \frac{d}{dx}=\frac{\mathcal{A}\mathcal{S}}{\Delta x}\left(\frac{1+\frac{1}{30}\mathcal{S}^2}{1+\frac{1}{5}\mathcal{S}^2}\right)

を$f_i$に作用させ，右辺の分母を左辺にもっていくと，

\left(1+\frac{1}{5}\mathcal{S}^2\right)f'_i=\frac{1}{\Delta x}\left(1+\frac{1}{30}\mathcal{S}^2\right)\mathcal{A}\mathcal{S}f_i

が得られる．演算子を全て展開し，$f'_i$の係数が$1$になるように変形すれば，6次精度のPadéスキームが得られる．

\frac{1}{3}f'_{i-1}+f'_i+\frac{1}{3}f'_{i+1}=\frac{14}{9}\frac{f_{i+1}-f_{i-1}}{2\Delta x}+\frac{1}{9}\frac{f_{i+2}-f_{i-2}}{4\Delta x}

文献³にはPadé近似の係数表が用意されているので，係数をいちいち求めなくても表を見るだけでよい⁵．

まとめ

Q. コンパクトスキームはなぜPadéスキームと呼ばれるか？
A. 微分の離散式の高次精度化において，ステンシル幅を狭めるためにPadé近似を使っているから．

Q. 演算子でこんな計算していいの？数学的に保証されてる？
A. Heaviside「私は消化のプロセスを知らないからといって食事をしないわけではない」

Padé近似の係数を計算するPythonコード

import mpmath as mp
import sympy as sp

def f(x):
    return 2/x*sp.asinh(x/2)*1/sp.sqrt(1+x**2/4)

def PadeApproximant(r,x,m,n):
    
    #Taylor展開で0からm+n次の項まで求める
    R = sp.series(r, x, 0, m+n+1).removeO() #ランダウの記号が入っていると係数を正しく取り出せないので削除する
    
    #Taylor展開した各項の係数の取り出し
    Rcoef = []
    for i in R.as_poly().all_coeffs()[::]:  
        Rcoef.append(mp.mpf(i))

    #高次項から順に格納されているので低次項を先に持ってくる  
    Rcoef = Rcoef[::-1]
    
    #all_coeffs()は最高次の係数が0のときに省略するっぽいので，リストの長さが1だけ短いときは0を追加する
    if len(Rcoef) == m+n:
        Rcoef.append(mp.mpf(0))
    print(Rcoef)
    
    #Pade近似
    p, q = mp.pade(Rcoef, m, n)
    return p,q

def main():
    #A(x)Q(x)=P(x)+O(x**(M+N+1)), Qᵢ(i=0 to N),Pᵢ(i=0 to M)
    M=0
    N=2
    
    x = sp.Symbol("x")
    P, Q = PadeApproximant(f(x),x,M,N)
    print(P[::],Q[::])
    # M=0,N=2のとき[mpf('1.0')] [mpf('1.0'), mpf('0.0'), mpf('0.16666666666666666')]
    # M=2,N=2のとき[mpf('1.0'), mpf('0.0'), mpf('0.033333333333333354')] [mpf('1.0'), mpf('0.0'), mpf('0.20000000000000001')]

if __name__ == "__main__":
    main()

1. Lele, S.K. (1992): Compact finite difference schemes with spectral-like resolution. Journal of Computational Physics 103, pp. 16–42. ↩ ↩² ↩³

2. Collatz, L. (1966):The numerical treatment of differential equations. Springer-Verlag, New York. ↩

3. Kopal, Z. (1961): Numerical analysis. Wiley, New York. ↩ ↩² ↩³ ↩⁴

4. Wolfram AlphaにSeries sinh^-1(x)を計算してもらえばよい．Sympyでもできる． ↩

5. 文献として参照するために海外の図書館の払い下げ品を購入したのだが，とてもカビ臭い． ↩