You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
69 lines
1.8 KiB
Python
69 lines
1.8 KiB
Python
6 years ago
|
from __future__ import division, print_function, absolute_import
|
||
|
|
|
||
|
|
from numpy import zeros, asarray, eye, poly1d, hstack, r_
|
||
|
|
from scipy import linalg
|
||
|
|
|
||
|
|
__all__ = ["pade"]
|
||
|
|
|
||
|
|
def pade(an, m, n=None):
|
||
|
|
"""
|
||
|
|
Return Pade approximation to a polynomial as the ratio of two polynomials.
|
||
|
|
|
||
|
|
Parameters
|
||
|
|
----------
|
||
|
|
an : (N,) array_like
|
||
|
|
Taylor series coefficients.
|
||
|
|
m : int
|
||
|
|
The order of the returned approximating polynomial `q`.
|
||
|
|
n : int, optional
|
||
|
|
The order of the returned approximating polynomial `p`. By default,
|
||
|
|
the order is ``len(an)-m``.
|
||
|
|
|
||
|
|
Returns
|
||
|
|
-------
|
||
|
|
p, q : Polynomial class
|
||
|
|
The Pade approximation of the polynomial defined by `an` is
|
||
|
|
``p(x)/q(x)``.
|
||
|
|
|
||
|
|
Examples
|
||
|
|
--------
|
||
|
|
>>> from scipy.interpolate import pade
|
||
|
|
>>> e_exp = [1.0, 1.0, 1.0/2.0, 1.0/6.0, 1.0/24.0, 1.0/120.0]
|
||
|
|
>>> p, q = pade(e_exp, 2)
|
||
|
|
|
||
|
|
>>> e_exp.reverse()
|
||
|
|
>>> e_poly = np.poly1d(e_exp)
|
||
|
|
|
||
|
|
Compare ``e_poly(x)`` and the Pade approximation ``p(x)/q(x)``
|
||
|
|
|
||
|
|
>>> e_poly(1)
|
||
|
|
2.7166666666666668
|
||
|
|
|
||
|
|
>>> p(1)/q(1)
|
||
|
|
2.7179487179487181
|
||
|
|
|
||
|
|
"""
|
||
|
|
an = asarray(an)
|
||
|
|
if n is None:
|
||
|
|
n = len(an) - 1 - m
|
||
|
|
if n < 0:
|
||
|
|
raise ValueError("Order of q <m> must be smaller than len(an)-1.")
|
||
|
|
if n < 0:
|
||
|
|
raise ValueError("Order of p <n> must be greater than 0.")
|
||
|
|
N = m + n
|
||
|
|
if N > len(an)-1:
|
||
|
|
raise ValueError("Order of q+p <m+n> must be smaller than len(an).")
|
||
|
|
an = an[:N+1]
|
||
|
|
Akj = eye(N+1, n+1)
|
||
|
|
Bkj = zeros((N+1, m), 'd')
|
||
|
|
for row in range(1, m+1):
|
||
|
|
Bkj[row,:row] = -(an[:row])[::-1]
|
||
|
|
for row in range(m+1, N+1):
|
||
|
|
Bkj[row,:] = -(an[row-m:row])[::-1]
|
||
|
|
C = hstack((Akj, Bkj))
|
||
|
|
pq = linalg.solve(C, an)
|
||
|
|
p = pq[:n+1]
|
||
|
|
q = r_[1.0, pq[n+1:]]
|
||
|
|
return poly1d(p[::-1]), poly1d(q[::-1])
|
||
|
|
|