Coverage for pySDC/projects/RDC/equidistant_RDC.py: 93%

1import numpy as np 

2from scipy.special import roots_legendre 

3from scipy.interpolate import BarycentricInterpolator 

4 

5from pySDC.core.Errors import CollocationError, ParameterError 

6from pySDC.core.Collocation import CollBase 

7 

8 

9class MyBarycentricInterpolator(BarycentricInterpolator): 

10 """ 

11 Overwrite BarycentricInterolator to inject custom weights 

12 """ 

13 

14 def __init__(self, xi, yi=None, weights=None, axis=0): 

15 super(MyBarycentricInterpolator, self).__init__(xi, yi, axis) 

16 self.wi = weights 

17 

18 

19class Equidistant_RDC(CollBase): 

20 """ 

21 Implements equidistant nodes with blended barycentric interpolation 

22 

23 Attributes: 

24 fh_weights: blended FH weights for barycentric interpolation 

25 """ 

26 

27 def __init__(self, num_nodes, tleft=0, tright=1, **kwargs): 

28 """ 

29 Initialization 

30 

31 Args: 

32 num_nodes: number of nodes 

33 tleft (float): left interval boundary (usually 0) 

34 tright (float): right interval boundary (usually 1) 

35 """ 

36 

37 if type(num_nodes) is int: 

38 max_d = 15 

39 nnodes = num_nodes 

40 else: 

41 if type(num_nodes) is not tuple: 

42 raise ParameterError('Expecting int or tuple for num_nodes parameter, got %s' % type(num_nodes)) 

43 if len(num_nodes) != 2: 

44 raise ParameterError('Expecting 1 or 2 arguments for num_nodes, got %s' % num_nodes) 

45 if type(num_nodes[0]) is not int: 

46 raise ParameterError('Expecting int type for first num_nodes argument, got %s' % type(num_nodes[0])) 

47 if type(num_nodes[1]) is not int: 

48 raise ParameterError('Expecting int type for second num_nodes argument, got %s' % type(num_nodes[1])) 

49 max_d = num_nodes[1] 

50 nnodes = num_nodes[0] 

51 

52 if nnodes < 2: 

53 raise CollocationError("Number of nodes should be at least 2 for equidistant, but is %d" % num_nodes) 

54 

55 try: 

56 super(Equidistant_RDC, self).__init__(num_nodes=nnodes, node_type='EQUID', quad_type='LOBATTO', **kwargs) 

57 except AttributeError: 

58 pass 

59 

60 self.order = self.num_nodes 

61 self.nodes = self._getNodes 

62 

63 d = min(self.num_nodes - 1, max_d) 

64 self.fh_weights = self._getFHWeights(d) 

65 self.weights = self._getWeights(tleft, tright) 

66 

67 self.Qmat = self._gen_Qmatrix 

68 self.Smat = self._gen_Smatrix 

69 self.delta_m = self._gen_deltas 

70 self.left_is_node = True 

71 self.right_is_node = True 

72 

73 def _getFHWeights(self, d): 

74 """ 

75 Computes blended FH weights for barycentric interpolation 

76 

77 This method is ported from Georges Klein's matlab function 

78 

79 Args: 

80 d (int): blending parameter 

81 

82 Returns: 

83 numpy.ndarray: weights 

84 """ 

85 

86 n = self.num_nodes - 1 

87 w = np.zeros(n + 1) 

88 

89 for k in range(0, n + 1): 

90 ji = max(k - d, 0) 

91 jf = min(k, n - d) 

92 sumcoeff = [] 

93 for i in range(ji, jf + 1): 

94 prodterm = [] 

95 for j in range(i, i + d + 1): 

96 if j == k: 

97 prodterm.append(1) 

98 else: 

99 prodterm.append(self.nodes[k] - self.nodes[j]) 

100 product = 1.0 / np.prod(prodterm) 

101 sumcoeff.append((-1) ** (i - 1) * product) 

102 y = sorted(sumcoeff, key=abs) 

103 w[k] = np.sum(y) 

104 

105 return w 

106 

107 def _getWeights(self, a, b): 

108 """ 

109 Computes weights using custom barycentric interpolation 

110 

111 Args: 

112 a (float): left interval boundary 

113 b (float): right interval boundary 

114 

115 Returns: 

116 numpy.ndarray: weights of the collocation formula given by the nodes 

117 """ 

118 if self.nodes is None: 

119 raise CollocationError("Need nodes before computing weights, got %s" % self.nodes) 

120 

121 circ_one = np.zeros(self.num_nodes) 

122 circ_one[0] = 1.0 

123 tcks = [] 

124 for i in range(self.num_nodes): 

125 tcks.append(MyBarycentricInterpolator(self.nodes, np.roll(circ_one, i), self.fh_weights)) 

126 

127 # Generate evaluation points for quadrature 

128 tau, omega = roots_legendre(self.num_nodes) 

129 phi = (b - a) / 2 * tau + (b + a) / 2 

130 

131 weights = [np.sum((b - a) / 2 * omega * p(phi)) for p in tcks] 

132 weights = np.array(weights) 

133 

134 return weights 

135 

136 @property 

137 def _gen_Qmatrix(self): 

138 """ 

139 Compute tleft-to-node integration matrix for later use in collocation formulation 

140 

141 Returns: 

142 numpy.ndarray: matrix containing the weights for tleft to node 

143 """ 

144 if self.nodes is None: 

145 raise CollocationError(f"Need nodes before computing weights, got {self.nodes}") 

146 M = self.num_nodes 

147 Q = np.zeros([M + 1, M + 1]) 

148 

149 # Generate Lagrange polynomials associated to the nodes 

150 circ_one = np.zeros(self.num_nodes) 

151 circ_one[0] = 1.0 

152 tcks = [] 

153 for i in range(M): 

154 tcks.append(MyBarycentricInterpolator(self.nodes, np.roll(circ_one, i), self.fh_weights)) 

155 

156 # Generate evaluation points for quadrature 

157 a, b = self.tleft, self.nodes[:, None] 

158 tau, omega = roots_legendre(self.num_nodes) 

159 tau, omega = tau[None, :], omega[None, :] 

160 phi = (b - a) / 2 * tau + (b + a) / 2 

161 

162 # Compute quadrature 

163 intQ = np.array([np.sum((b - a) / 2 * omega * p(phi), axis=-1) for p in tcks]) 

164 

165 # Store into Q matrix 

166 Q[1:, 1:] = intQ.T 

167 

168 return Q