アンリ・パデ 数学 においてパデ近似 (パデきんじ、英 : Padé approximant )とは、関数を近似する「最良」の有理関数 のこと。たとえば x / ( 1 + 1 2 x ) {\displaystyle x/(1+{\scriptstyle {\frac {1}{2}}}x)} は log (1 + x ) のパデ近似のひとつである:
log ( 1 + x ) = x 1 + 1 2 x + O ( x 3 ) . {\displaystyle \log(1+x)={\frac {x}{1+{\scriptstyle {\frac {1}{2}}}x}}+O(x^{3}).} パデ近似のテイラー級数 は関数のテイラー級数と与えられた次数まで一致する。この近似法は1890年 頃にアンリ・パデ(英語版) が発展させたが、冪級数 の有理関数による近似という考えを始め、その特徴を研究したのはゲオルク・フロベニウス にまで遡る。
多くの場合、パデ近似は、 テイラー級数を有限項で切り捨てるより良い近似を与えるが、テイラー級数が収束 しない場合でも機能する。 これらの理由から、パデ近似はコンピューター計算で広く使用されている。 また、パデ近似は ディオファントス近似 および超越数論において補助関数(英語版) として使用されるが、より正確な評価のためにはパデ近似を応用したアドホックな手法を使うことが一般的である。
また、有理関数であるがために近似として人工的な特異点が発生するおそれがあるが、これはボレル・パデ解析 によって回避することができる。
パデ近似がマクローリン展開 より良い近似になりやすい理由は、多点総和法の観点から見れば明らかである。それは無限遠での漸近展開が0や定数になる例が多いため、「不完全な2点パデ近似」として、通常のパデ近似がマクローリン展開を改良しているものと解釈出来るからである。
定義 滑らかな関数 f (x ) と非負整数 m , n に対して f (x ) の [m /n ] 次パデ近似 とは有理関数
R ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a m x m 1 + b 1 x + b 2 x 2 + ⋯ + b n x n {\displaystyle R(x)={\frac {a_{0}+a_{1}x+a_{2}x^{2}+\dots +a_{m}x^{m}}{1+b_{1}x+b_{2}x^{2}+\dots +b_{n}x^{n}}}} であって
f ( k ) ( 0 ) = R ( k ) ( 0 ) ( 0 ≤ k ≤ m + n ) {\displaystyle f^{(k)}(0)=R^{(k)}(0)\qquad (0\leq k\leq m+n)} を満たすものをいう。(特に n = 0 のときは原点における m 次テイラー多項式 に他ならない。)
有理関数 R (x ) の原点におけるテイラー級数 は先頭の m + n + 1 項が f (x ) のそれと相殺され、
f ( x ) = R ( x ) + O ( x m + n + 1 ) {\displaystyle f(x)=R(x)+O(x^{m+n+1})} となる。
パデ近似は与えられた非負整数 m , n に対して(存在すれば)一意的に決まる——つまり係数 a 0 , a 1 , … , a m , b 1 , … , b n {\displaystyle a_{0},a_{1},\dots ,a_{m},b_{1},\dots ,b_{n}} は一意的に決まる。パデ近似 R (x ) における分母の定数項として 1 を選ぶのはこの一意性のためであり、この正規化をしないとき分母と分子の一意性は定数倍の違いを除いて しか成り立たない。
このように定義されたパデ近似 R (x ) を
[ m / n ] f ( x ) {\displaystyle [m/n]_{f}(x)} と表し(関数 f や変数 x は省略されることもある)、これらを並べた表をパデ表(英語版) という。
指数関数 exp(x ) のパデ表 [m /n ](x )の一部 m \ n 0 1 2 3 0 1 1 {\displaystyle {\frac {1}{1}}} 1 1 − x {\displaystyle {\frac {1}{1-x}}} 1 1 − x + 1 2 x 2 {\displaystyle {\frac {1}{1-x+{\scriptstyle {\frac {1}{2}}}x^{2}}}} 1 1 − x + 1 2 x 2 − 1 6 x 3 {\displaystyle {\frac {1}{1-x+{\scriptstyle {\frac {1}{2}}}x^{2}-{\scriptstyle {\frac {1}{6}}}x^{3}}}} 1 1 + x 1 {\displaystyle {\frac {1+x}{1}}} 1 + 1 2 x 1 − 1 2 x {\displaystyle {\frac {1+{\scriptstyle {\frac {1}{2}}}x}{1-{\scriptstyle {\frac {1}{2}}}x}}} 1 + 1 3 x 1 − 2 3 x + 1 6 x 2 {\displaystyle {\frac {1+{\scriptstyle {\frac {1}{3}}}x}{1-{\scriptstyle {\frac {2}{3}}}x+{\scriptstyle {\frac {1}{6}}}x^{2}}}} 1 + 1 4 x 1 − 3 4 x + 1 4 x 2 − 1 24 x 3 {\displaystyle {\frac {1+{\scriptstyle {\frac {1}{4}}}x}{1-{\scriptstyle {\frac {3}{4}}}x+{\scriptstyle {\frac {1}{4}}}x^{2}-{\scriptstyle {\frac {1}{24}}}x^{3}}}} 2 1 + x + 1 2 x 2 1 {\displaystyle {\frac {1+x+{\scriptstyle {\frac {1}{2}}}x^{2}}{1}}} 1 + 2 3 x + 1 6 x 2 1 − 1 3 x {\displaystyle {\frac {1+{\scriptstyle {\frac {2}{3}}}x+{\scriptstyle {\frac {1}{6}}}x^{2}}{1-{\scriptstyle {\frac {1}{3}}}x}}} 1 + 1 2 x + 1 12 x 2 1 − 1 2 x + 1 12 x 2 {\displaystyle {\frac {1+{\scriptstyle {\frac {1}{2}}}x+{\scriptstyle {\frac {1}{12}}}x^{2}}{1-{\scriptstyle {\frac {1}{2}}}x+{\scriptstyle {\frac {1}{12}}}x^{2}}}} 1 + 2 5 x + 1 20 x 2 1 − 3 5 x + 3 20 x 2 − 1 60 x 3 {\displaystyle {\frac {1+{\scriptstyle {\frac {2}{5}}}x+{\scriptstyle {\frac {1}{20}}}x^{2}}{1-{\scriptstyle {\frac {3}{5}}}x+{\scriptstyle {\frac {3}{20}}}x^{2}-{\scriptstyle {\frac {1}{60}}}x^{3}}}} 3 1 + x + 1 2 x 2 + 1 6 x 3 1 {\displaystyle {\frac {1+x+{\scriptstyle {\frac {1}{2}}}x^{2}+{\scriptstyle {\frac {1}{6}}}x^{3}}{1}}} 1 + 3 4 x + 1 4 x 2 + 1 24 x 3 1 − 1 4 x {\displaystyle {\frac {1+{\scriptstyle {\frac {3}{4}}}x+{\scriptstyle {\frac {1}{4}}}x^{2}+{\scriptstyle {\frac {1}{24}}}x^{3}}{1-{\scriptstyle {\frac {1}{4}}}x}}} 1 + 3 5 x + 3 20 x 2 + 1 60 x 3 1 − 2 5 x + 1 20 x 2 {\displaystyle {\frac {1+{\scriptstyle {\frac {3}{5}}}x+{\scriptstyle {\frac {3}{20}}}x^{2}+{\scriptstyle {\frac {1}{60}}}x^{3}}{1-{\scriptstyle {\frac {2}{5}}}x+{\scriptstyle {\frac {1}{20}}}x^{2}}}} 1 + 1 2 x + 1 10 x 2 + 1 120 x 3 1 − 1 2 x + 1 10 x 2 − 1 120 x 3 {\displaystyle {\frac {1+{\scriptstyle {\frac {1}{2}}}x+{\scriptstyle {\frac {1}{10}}}x^{2}+{\scriptstyle {\frac {1}{120}}}x^{3}}{1-{\scriptstyle {\frac {1}{2}}}x+{\scriptstyle {\frac {1}{10}}}x^{2}-{\scriptstyle {\frac {1}{120}}}x^{3}}}}
たとえば指数関数 exp(x ) の [m /n ] 次パデ近似は一般化された超幾何関数 を用いて
[ m / n ] ( x ) = 1 F 1 ( − m − m − n ; x ) / 1 F 1 ( − n − m − n ; − x ) {\displaystyle [m/n](x)={}_{1}F_{1}\left({\scriptstyle {\begin{matrix}-m\\-m-n\end{matrix}}};x\right){\big /}{}_{1}F_{1}\left({\scriptstyle {\begin{matrix}-n\\-m-n\end{matrix}}};-x\right)} と表される。
計算 指定されたx に対して 、パデ近似はWynnのイプシロン・アルゴリズム[1] によって計算できるが、f のテイラー級数 の部分和
T N ( x ) = c 0 + c 1 x + c 2 x 2 + ⋯ + c N x N {\displaystyle T_{N}(x)=c_{0}+c_{1}x+c_{2}x^{2}+\cdots +c_{N}x^{N}} から数列の変形によって計算することもできる[2] 。ここで
c k = f ( k ) ( 0 ) k ! . {\displaystyle c_{k}={\frac {f^{(k)}(0)}{k!}}.} f は形式的なべき級数 としてもよく、そのため、パデ近似を発散級数 の総和をとるという目的で使用することもできる。
パデ近似を計算する1つの方法は、 多項式最大公約数の拡張ユークリッドアルゴリズムを使用することである。 [3] 関係
R ( x ) = P ( x ) / Q ( x ) = T m + n ( x ) mod x m + n + 1 {\displaystyle R(x)=P(x)/Q(x)=T_{m+n}(x){\text{ mod }}x^{m+n+1}} は、次のような因子K (x ) の存在と同値である:
P ( x ) = Q ( x ) T m + n ( x ) + K ( x ) x m + n + 1 . {\displaystyle P(x)=Q(x)T_{m+n}(x)+K(x)x^{m+n+1}.} これは、 T m + n ( x ) {\displaystyle T_{m+n}(x)} と x m + n + 1 {\displaystyle x^{m+n+1}} の最大公約数を求める計算における1つのステップのベズー恒等式 として解釈できる。
2つの多項式 p と q の最大公約数を計算するには、筆算によって余りの列
r 0 = p , r 1 = q , r k − 1 = q k r k + r k + 1 , deg r k + 1 < deg r k ( k = 1 , 2 , 3... ) {\displaystyle r_{0}=p,\;r_{1}=q,\quad r_{k-1}=q_{k}r_{k}+r_{k+1},\;\deg r_{k+1}<\deg r_{k}\;(k=1,2,3...)} を r k + 1 = 0 {\displaystyle r_{k+1}=0} となるまで計算したことを思い出す。 拡張最大公約数のベズー恒等式では、2つの多項式列
u 0 = 1 , v 0 = 0 , u 1 = 0 , v 1 = 1 , u k + 1 = u k − 1 − q k u k , v k + 1 = v k − 1 − q k v k {\displaystyle u_{0}=1,\;v_{0}=0,\quad u_{1}=0,\;v_{1}=1,\quad u_{k+1}=u_{k-1}-q_{k}u_{k},\;v_{k+1}=v_{k-1}-q_{k}v_{k}} を同時に計算する。これによって、各ステップでベズー恒等式
r k ( x ) = u k ( x ) p ( x ) + v k ( x ) q ( x ) . {\displaystyle r_{k}(x)=u_{k}(x)p(x)+v_{k}(x)q(x).} を得る。
[m /n ]近似の場合、次の拡張ユークリッドアルゴリズムを実行する。
r 0 = x m + n + 1 , r 1 = T m + n ( x ) {\displaystyle r_{0}=x^{m+n+1},\;r_{1}=T_{m+n}(x)} そして v k {\displaystyle v_{k}} の次数がn 以下である最後の段階にそれを停止する。
次に、多項式 P = r k , Q = v k {\displaystyle P=r_{k},\;Q=v_{k}} [ m / n ]パデ近似を与える。拡張最大公約数計算のすべてのステップを計算すると、 パデ表の対角線が得られる。
リーマン・パデゼータ関数 発散級数 の再足し上げ(resummation)、すなわち
∑ z = 1 ∞ f ( z ) {\displaystyle \sum _{z=1}^{\infty }f(z)} を調べるには、パデまたは単に有理ゼータ関数を次のように導入すると便利である。
ζ R ( s ) = ∑ z = 1 ∞ R ( z ) z s , {\displaystyle \zeta _{R}(s)=\sum _{z=1}^{\infty }{\frac {R(z)}{z^{s}}},} ここで、
R ( x ) = [ m / n ] f ( x ) {\displaystyle R(x)=[m/n]_{f}(x)\,} は関数f (x )の次数 (m ,n ) のパデ近似である。 s = 0での ゼータ正則化 値は、発散級数の和と見なされる。
このパデゼータ関数の関数方程式は次のとおりである。
∑ j = 0 n a j ζ R ( s − j ) = ∑ j = 0 m b j ζ 0 ( s − j ) , {\displaystyle \sum _{j=0}^{n}a_{j}\zeta _{R}(s-j)=\sum _{j=0}^{m}b_{j}\zeta _{0}(s-j),} ここで、a j とb j はパデ近似の係数である。下付き文字「0」は、パデが次数[0/0]であることを意味する。したがって、この場合リーマンゼータ関数となる。
DLogPadéメソッド パデ近似を使用すると、関数の臨界点と指数を抽出できる。熱力学では、関数 f (x ) が点x = r の近くで f ( x ) ∼ | x − r | p {\displaystyle f(x)\sim |x-r|^{p}} のように非解析的にふるまうとき、 x =r を臨界点、 p をf の関連する臨界指数と呼ぶ。f の級数展開の十分な項が分かっている場合、パデ近似 [ n / n + 1 ] g ( x ) {\displaystyle [n/n+1]_{g}(x)} の極と残差から臨界点と臨界指数をそれぞれ見積もることができる。 ここで、 g = f ′ f {\displaystyle g={\frac {f'}{f}}} である。
一般化 パデ近似は、1つの変数で関数を近似する。2つの変数による近似は、チザム近似(J.S.R.チザムにちなむ) [4] と呼ばれ、複数の変数による近似はカンタベリー近似(カンタベリーにあるケント大学にいたグレイブス゠モリスにちなむ) [5] と呼ばれる。
2点パデ近似 従来のパデ近似は、マクローリン展開を与えられた次数まで再現するように決定されている。そのため、展開点から離れた箇所での値での近似が悪くなることがある。これを回避するのが多点総和法の1種である2点パデ近似である[6] 。 x = 0 {\displaystyle x=0} で、関数 f ( x ) {\displaystyle f(x)} がある漸近関数 f 0 ( x ) {\displaystyle f_{0}(x)} を用いて、
f ∼ f 0 ( x ) + o ( f 0 ( x ) ) ( x → 0 ) {\displaystyle f\sim f_{0}(x)+o(f_{0}(x))(x\rightarrow 0)}
と表され、更に、 x → ∞ {\displaystyle x\rightarrow \infty } では、ある漸近関数 f ∞ ( x ) {\displaystyle f_{\infty }(x)} を用いて、
f ( x ) ∼ f ∞ ( x ) + o ( f ∞ ( x ) ) ( x → ∞ ) {\displaystyle f(x)\sim f_{\infty }(x)+o(f_{\infty }(x))(x\rightarrow \infty )}
と表される場合を考える。適切に f 0 ( x ) , f ∞ ( x ) {\displaystyle f_{0}(x),f_{\infty }(x)} の主要なふるまいを選び出すことで、パデ近似を拡張して用いることにより、これらの漸近的振舞いを同時に再現する近似関数 F ( x ) {\displaystyle F(x)} を様々な場合に見つけることができる。これにより、通常のパデ近似で近似の精度が最も悪くなる恐れのある x → ∞ {\displaystyle x\rightarrow \infty } で、精度が良くなることが保証される。そのため、2点パデ近似は x = 0 ∼ ∞ {\displaystyle x=0\sim \infty } で大域的に良い近似を与える手法となりうる。
f 0 ( x ) , f ∞ ( x ) {\displaystyle f_{0}(x),f_{\infty }(x)} が多項式や負べきの級数で表される場合や、指数関数、対数関数で表される場合、 x ln x {\displaystyle x\ln x} と表される場合などに適用可能である。これを用いて微分方程式の近似解を精度よく与える方法が存在する[6] 。また、リーマンゼータ関数の非自明な零点についても実軸上の漸近的ふるまいから、最初の非自明な零点をある程度の精度で見積もることができる[6] 。
多点パデ近似 2点パデ近似をさらに拡張したものが多点パデ近似である[6] 。これは、 x = x j ( j = 1 , 2 , 3 ⋯ , N ) {\displaystyle x=x_{j}(j=1,2,3\cdots ,N)} において、近似したい関数 f ( x ) {\displaystyle f(x)} がは指数 n j {\displaystyle n_{j}} で表される特異点
f ( x ) ∼ A j ( x − x j ) n j ( x → x j ) {\displaystyle f(x)\sim {\frac {A_{j}}{(x-x_{j})^{n_{j}}}}(x\rightarrow x_{j})}
を持つ場合に、2点パデ近似の x = 0 , x → ∞ {\displaystyle x=0,x\rightarrow \infty } に加えて、これらの点 x ∼ x j {\displaystyle x\sim x_{j}} で発散するという性質を再現するように近似する方法である。これにより、関数の特異性の情報を取り込むため、さらに良い精度で関数 f ( x ) {\displaystyle f(x)} の近似が可能となる。
また、区間をいくつかの有限または半無限区間に分割することで、その区間を変数変換により、通常の2点パデ近似が適用可能な形式にすることができる。そのようにして、各区間ごとに得られた2点パデ近似を繋ぎ合わせたものを多点パデ近似と呼ぶこともある。
例 sin(x) [7] sin ( x ) ≈ ( 12671 / 4363920 ) x 5 − ( 2363 / 18183 ) x 3 + x 1 + ( 445 / 12122 ) x 2 + ( 601 / 872784 ) x 4 + ( 121 / 16662240 ) x 6 {\displaystyle \sin(x)\approx {\frac {(12671/4363920)x^{5}-(2363/18183)x^{3}+x}{1+(445/12122)x^{2}+(601/872784)x^{4}+(121/16662240)x^{6}}}} exp(x) [8] exp ( x ) ≈ 1 + ( 1 / 2 ) x + ( 1 / 9 ) x 2 + ( 1 / 72 ) x 3 + ( 1 / 1008 ) x 4 + ( 1 / 30240 ) x 5 1 − ( 1 / 2 ) x + ( 1 / 9 ) x 2 − ( 1 / 72 ) x 3 + ( 1 / 1008 ) x 4 − ( 1 / 30240 ) x 5 {\displaystyle \exp(x)\approx {\frac {1+(1/2)x+(1/9)x^{2}+(1/72)x^{3}+(1/1008)x^{4}+(1/30240)x^{5}}{1-(1/2)x+(1/9)x^{2}-(1/72)x^{3}+(1/1008)x^{4}-(1/30240)x^{5}}}} ヤコビの楕円関数 sn(z|3) [9] s n ( z | 3 ) ≈ − ( 9851629 / 283609260 ) z 5 − ( 572744 / 4726821 ) z 3 + z 1 + ( 859490 / 1575607 ) z 2 − ( 5922035 / 56721852 ) z 4 + ( 62531591 / 2977897230 ) z 6 {\displaystyle \mathrm {sn} (z|3)\approx {\frac {-(9851629/283609260)z^{5}-(572744/4726821)z^{3}+z}{1+(859490/1575607)z^{2}-(5922035/56721852)z^{4}+(62531591/2977897230)z^{6}}}} ベッセル関数 J5 (x) J 5 ( x ) ≈ − ( 107 / 28416000 ) x 7 + ( 1 / 3840 ) x 5 1 + ( 151 / 5550 ) x 2 + ( 1453 / 3729600 ) x 4 + ( 1339 / 358041600 ) x 6 + ( 2767 / 120301977600 ) x 8 {\displaystyle J_{5}(x)\approx {\frac {-(107/28416000)x^{7}+(1/3840)x^{5}}{1+(151/5550)x^{2}+(1453/3729600)x^{4}+(1339/358041600)x^{6}+(2767/120301977600)x^{8}}}} erf(x) erf ( x ) ≈ ( 2 / 15 ) ⋅ ( 49140 x + 3570 x 3 + 739 x 5 ) π ⋅ ( 165 x 4 + 1330 x 2 + 3276 ) {\displaystyle \operatorname {erf} (x)\approx {\frac {(2/15)\cdot (49140x+3570x^{3}+739x^{5})}{{\sqrt {\pi }}\cdot (165x^{4}+1330x^{2}+3276)}}} フレネル積分 C(x) C ( x ) ≈ ( 1 / 135 ) ⋅ ( 990791 x 9 π 4 − 147189744 x 5 π 2 + 8714684160 x ) ( 1749 π 4 x 8 + 523536 π 2 x 4 + 64553216 ) {\displaystyle C(x)\approx {\frac {(1/135)\cdot (990791x^{9}\pi ^{4}-147189744x^{5}\pi ^{2}+8714684160x)}{(1749\pi ^{4}x^{8}+523536\pi ^{2}x^{4}+64553216)}}}
脚注 [脚注の使い方 ]
出典 ^ Theorem 1 in Wynn, Peter (Mar 1966), “On the Convergence and Stability of the Epsilon Algorithm”, SIAM Journal on Numerical Analysis 3 (1): 91–122, Bibcode: 1966SJNA....3...91W, doi:10.1137/0703007, JSTOR 2949688, https://jstor.org/stable/2949688 ^ Brezenski, C. (1996), “Extrapolation algorithms and Padé approximations”, Applied Numerical Mathematics 20 (3): 299–318, doi:10.1016/0168-9274(95)00110-7 ^ Problem 5.2b and Algorithm 5.2 (p. 46) in Bini, Dario; Pan, Victor (1994), Polynomial and Matrix computations - Volume 1. Fundamental Algorithms , Progress in Theoretical Computer Science, Birkhäuser, ISBN 978-0-8176-3786-6 ^ Chisholm, J. S. R. (1973). “Rational approximants defined from double power series”. Mathematics of Computation 27 (124): 841–848. doi:10.1090/S0025-5718-1973-0382928-6. ISSN 0025-5718. ^ Graves-Morris, P.R.; Roberts, D.E. (1975). “Calculation of Canterbury approximants”. Computer Physics Communications 10 (4): 234–244. Bibcode: 1975CoPhC..10..234G. doi:10.1016/0010-4655(75)90068-5. ^ a b c d 上岡, 良季. 多点総和法入門 高校生でもわかる!!ココと無限のかなたをつなぐ現代応用数学: テイラー展開から微分方程式の応用まで . https://www.amazon.co.jp/dp/B08DV9TZVD/ ^ “sin(x)のパデ近似”. Wolfram Alpha . 2022年1月16日 閲覧。 ^ “exp(x)のパデ近似”. Wolfram Alpha . 2022年1月16日 閲覧。 ^ “sn(x|3)のパデ近似”. Wolfram Alpha . 2022年1月16日 閲覧。
参考文献 Baker, G. A., Jr.; and Graves-Morris, P. Padé Approximants . Cambridge U.P., 1996 Baker, G. A., Jr. Padé approximant, Scholarpedia , 7(6):9756. doi:10.4249/scholarpedia.9756 Brezinski, C.; and Redivo Zaglia, M. Extrapolation Methods. Theory and Practice . North-Holland, 1991 Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), “Section 5.12 Padé Approximants”, Numerical Recipes: The Art of Scientific Computing (3rd ed.), New York: Cambridge University Press, ISBN 978-0-521-88068-8, http://apps.nrbook.com/empanel/index.html?pg=245 Frobenius, G.; Ueber Relationen zwischen den Näherungsbrüchen von Potenzreihen , [Journal für die reine und angewandte Mathematik (Crelle's Journal)]. Volume 1881, Issue 90, Pages 1–17 Gragg, W.B.; The Pade Table and Its Relation to Certain Algorithms of Numerical Analysis [SIAM Review], Vol. 14, No. 1, 1972, pp. 1–62. Padé, H.; Sur la répresentation approchée d'une fonction par des fractions rationelles , Thesis, [Ann. \'Ecole Nor. (3), 9, 1892, pp. 1–93 supplement. Wynn, P. (1966), “Upon systems of recursions which obtain among the quotients of the Padé table”, Numerische Mathematik 8 (3): 264–269, doi:10.1007/BF02162562
関連項目 パデ表(英語版) バースカラ1世の正弦近似式(英語版)
外部リンク Weisstein, Eric W. "Padé Approximant". mathworld.wolfram.com (英語). Padé Approximants, Oleksandr Pavlyk, The Wolfram Demonstrations Project Data Analysis BriefBook: Pade Approximation, Rudolf K. Bock European Laboratory for Particle Physics, CERN Sinewave, Scott Dattalo, last accessed 2010-11-11. MATLAB function for Pade approximation of models with time delays.