Coverage for src/methodsnm/recpol.py: 55%

1from numpy import array 

2from methodsnm.fe import * 

3 

4import numpy as np 

5from abc import ABC, abstractmethod 

6 

7import matplotlib.pyplot as plt 

8 

9 

10class RecursivePolynomial(ABC): 

11 """ 

12 Abstract base class for recursive polynomial evaluation. 

13 

14 Polynomials are evaluated using the following recursion: 

15 * d "starter" (given functions) 

16 * P_i(x) = s_i(x) for i < d 

17 * recursion formulate for i >= d: 

18 * P_i(x) = sum_{j=0}^{d-1} c_j(i,x) P_{i-d+j}(x) 

19 

20 Attributes: 

21 ----------- 

22 starter : list 

23 List of starter functions. 

24 rec_coeff : list 

25 List of recursive coefficients. 

26 """ 

27 starter = None 

28 rec_coeff = None 

29 

30 @abstractmethod 

31 def __init__(self): 

32 pass 

33 

34 def evaluate_all(self, x, n): 

35 """ 

36 Evaluates the recursive polynomial for all values of x up to n. 

37 

38 Parameters: 

39 ----------- 

40 x : array_like 

41 Array of values to evaluate the polynomial at. 

42 n : int 

43 Maximum degree of the polynomial to evaluate. 

44 

45 Returns: 

46 -------- 

47 vals : ndarray 

48 Array of shape (len(x), n+1 ) containing the values of the polynomial 

49 evaluated at each value of x up to degree n. 

50 """ 

51 vals = np.empty((len(x), n+1), dtype=x.dtype) 

52 d = len(self.rec_coeff) 

53 for si, s in enumerate(self.starter): 

54 vals[:,si] = s(x) 

55 for i in range(d, n+1): 

56 vals[:,i] = sum([self.rec_coeff[j](i,x)*vals[:,i-d+j] for j in range(d)]) 

57 return vals 

58 

59 def evaluate(self, x, n): 

60 """ 

61 Evaluate the recursive polynomial of degree n for all values of x. 

62 

63 Parameters: 

64 ----------- 

65 x : array_like 

66 Array of values to evaluate the polynomial at. 

67 n : int 

68 Degree of the polynomial to evaluate. 

69 

70 Returns: 

71 -------- 

72 vals : ndarray 

73 Array of length len(x) containing the values of the polynomial 

74 evaluated at each value of x. 

75 """ 

76 d = len(self.rec_coeff) 

77 vals = np.empty((len(x), d+1)) 

78 for si, s in enumerate(self.starter): 

79 vals[:,si] = s(x) 

80 for i in range(d, n+1): 

81 if i != d: 

82 for j in range(d): 

83 vals[:,j] = vals[:,j+1] 

84 vals[:,d] = sum([self.rec_coeff[j](i,x)*vals[:,j] for j in range(d)]) 

85 return vals[:,d] 

86 

87 def plot_all(self, x, n): 

88 """ 

89 Plots the recursive polynomial for all values of x up to n. 

90 

91 Parameters: 

92 ----------- 

93 x : array_like 

94 Array of values to evaluate the polynomial at. 

95 n : int 

96 Maximum degree of the polynomial to evaluate. 

97 """ 

98 vals = self.evaluate_all(x, n) 

99 for i in range(n+1): 

100 plt.plot(x, vals[:,i], label="P_{}".format(i)) 

101 plt.legend() 

102 plt.show() 

103 

104 

105class Monomials(RecursivePolynomial): 

106 def __init__(self): 

107 self.starter = [lambda x: np.ones_like(x)] 

108 self.rec_coeff = [lambda n,x: x] 

109 

110class LegendrePolynomials(RecursivePolynomial): 

111 def __init__(self): 

112 self.starter = [lambda x: x**0, lambda x: x] 

113 self.rec_coeff = [lambda n,x: -(n-1)/n, lambda n,x: (2*n-1)/n*x] 

114 

115class JacobiPolynomials(RecursivePolynomial): 

116 alpha = None 

117 beta = None 

118 def __init__(self, alpha, beta): 

119 self.alpha = alpha 

120 self.beta = beta 

121 self.starter = [lambda x: np.ones_like(x), lambda x: alpha+1+(alpha+beta+2)/2*(x-1)] 

122 def c0(n,x): 

123 alpha, beta = self.alpha, self.beta 

124 return -2*(n+alpha-1)*(n+beta-1)*(2*n+alpha+beta)/((2*n+alpha+beta-2)*(n+alpha+beta)*2*n) 

125 def c1(n,x): 

126 alpha, beta = self.alpha, self.beta 

127 return (2*n+alpha+beta-1)*((2*n+alpha+beta)*(2*n+alpha+beta-2)*x+alpha**2-beta**2)/((2*n+alpha+beta-2)*(n+alpha+beta)*2*n) 

128 self.rec_coeff = [c0, c1] 

129 

130class IntegratedLegendrePolynomials(RecursivePolynomial): 

131 def __init__(self): 

132 self.starter = [lambda x: x, lambda x: 0.5*(x**2-1)] 

133 self.rec_coeff = [lambda n,x: -(n-2)/n, lambda n,x: (2*n-3)/n*x]