Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/interpolate/_pade.py: 15%

26 statements  

« prev     ^ index     » next       coverage.py v7.3.2, created at 2023-12-12 06:31 +0000

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)-1-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 >>> import numpy as np 

29 >>> from scipy.interpolate import pade 

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

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

32 

33 >>> e_exp.reverse() 

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

35 

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

37 

38 >>> e_poly(1) 

39 2.7166666666666668 

40 

41 >>> p(1)/q(1) 

42 2.7179487179487181 

43 

44 """ 

45 an = asarray(an) 

46 if n is None: 

47 n = len(an) - 1 - m 

48 if n < 0: 

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

50 if n < 0: 

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

52 N = m + n 

53 if N > len(an)-1: 

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

55 an = an[:N+1] 

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

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

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

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

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

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

62 C = hstack((Akj, Bkj)) 

63 pq = linalg.solve(C, an) 

64 p = pq[:n+1] 

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

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

67