Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1from numpy import zeros, asarray, eye, poly1d, hstack, r_ 

2from scipy import linalg 

3 

4__all__ = ["pade"] 

5 

6def pade(an, m, n=None): 

7 """ 

8 Return Pade approximation to a polynomial as the ratio of two polynomials. 

9 

10 Parameters 

11 ---------- 

12 an : (N,) array_like 

13 Taylor series coefficients. 

14 m : int 

15 The order of the returned approximating polynomial `q`. 

16 n : int, optional 

17 The order of the returned approximating polynomial `p`. By default, 

18 the order is ``len(an)-m``. 

19 

20 Returns 

21 ------- 

22 p, q : Polynomial class 

23 The Pade approximation of the polynomial defined by `an` is 

24 ``p(x)/q(x)``. 

25 

26 Examples 

27 -------- 

28 >>> from scipy.interpolate import pade 

29 >>> e_exp = [1.0, 1.0, 1.0/2.0, 1.0/6.0, 1.0/24.0, 1.0/120.0] 

30 >>> p, q = pade(e_exp, 2) 

31 

32 >>> e_exp.reverse() 

33 >>> e_poly = np.poly1d(e_exp) 

34 

35 Compare ``e_poly(x)`` and the Pade approximation ``p(x)/q(x)`` 

36 

37 >>> e_poly(1) 

38 2.7166666666666668 

39 

40 >>> p(1)/q(1) 

41 2.7179487179487181 

42 

43 """ 

44 an = asarray(an) 

45 if n is None: 

46 n = len(an) - 1 - m 

47 if n < 0: 

48 raise ValueError("Order of q <m> must be smaller than len(an)-1.") 

49 if n < 0: 

50 raise ValueError("Order of p <n> must be greater than 0.") 

51 N = m + n 

52 if N > len(an)-1: 

53 raise ValueError("Order of q+p <m+n> must be smaller than len(an).") 

54 an = an[:N+1] 

55 Akj = eye(N+1, n+1, dtype=an.dtype) 

56 Bkj = zeros((N+1, m), dtype=an.dtype) 

57 for row in range(1, m+1): 

58 Bkj[row,:row] = -(an[:row])[::-1] 

59 for row in range(m+1, N+1): 

60 Bkj[row,:] = -(an[row-m:row])[::-1] 

61 C = hstack((Akj, Bkj)) 

62 pq = linalg.solve(C, an) 

63 p = pq[:n+1] 

64 q = r_[1.0, pq[n+1:]] 

65 return poly1d(p[::-1]), poly1d(q[::-1]) 

66