Coverage for pySDC/core/Nodes.py: 100%
78 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1import numpy as np
2from scipy.linalg import eigh_tridiagonal
4NODE_TYPES = ['EQUID', 'LEGENDRE', 'CHEBY-1', 'CHEBY-2', 'CHEBY-3', 'CHEBY-4']
6QUAD_TYPES = ['GAUSS', 'RADAU-LEFT', 'RADAU-RIGHT', 'LOBATTO']
9class NodesError(Exception):
10 """Exception class to handle error in NodesGenerator class"""
12 pass
15class NodesGenerator(object):
16 """
17 Class that can be used to generate generic distribution of nodes derived
18 from Gauss quadrature rule.
19 Its implementation is fully inspired from a `book of W. Gautschi <https://doi.org/10.1093/oso/9780198506720.001.0001>`_.
21 Attributes
22 ----------
23 node_type : str
24 The type of node distribution
25 quad_type : str
26 The quadrature type
28 """
30 def __init__(self, node_type='LEGENDRE', quad_type='LOBATTO'):
31 """
33 Parameters
34 ----------
35 node_type : str, optional
36 The type of node distribution, can be
38 - EQUID : equidistant nodes
39 - LEGENDRE : node distribution from Legendre polynomials
40 - CHEBY-1 : node distribution from Chebychev polynomials (1st kind)
41 - CHEBY-2 : node distribution from Chebychev polynomials (2nd kind)
42 - CHEBY-3 : node distribution from Chebychev polynomials (3rd kind)
43 - CHEBY-4 : node distribution from Chebychev polynomials (4th kind)
45 The default is 'LEGENDRE'.
47 quad_type : str, optional
48 The quadrature type, can be
50 - GAUSS : inner point only, no node at boundary
51 - RADAU-LEFT : only left boundary as node
52 - RADAU-RIGHT : only right boundary as node
53 - LOBATTO : left and right boundary as node
55 The default is 'LOBATTO'.
56 """
58 # Check argument validity
59 for arg, vals in zip(['node_type', 'quad_type'], [NODE_TYPES, QUAD_TYPES]):
60 val = eval(arg)
61 if val not in vals:
62 raise NodesError(f"{arg}='{val}' not implemented, must be in {vals}")
64 # Store attributes
65 self.node_type = node_type
66 self.quad_type = quad_type
68 def getNodes(self, num_nodes):
69 """
70 Computes a given number of quadrature nodes.
72 Parameters
73 ----------
74 num_nodes : int
75 Number of nodes to compute.
77 Returns
78 -------
79 nodes : np.1darray
80 Nodes located in [-1, 1], in increasing order.
81 """
82 # Check number of nodes
83 if self.quad_type in ['LOBATTO', 'RADAU-LEFT'] and num_nodes < 2:
84 raise NodesError(f"num_nodes must be larger than 2 for {self.quad_type}, but got {num_nodes}")
85 elif num_nodes < 1:
86 raise NodesError("you surely want at least one node ;)")
88 # Equidistant nodes
89 if self.node_type == 'EQUID':
90 if self.quad_type == 'GAUSS':
91 return np.linspace(-1, 1, num=num_nodes + 2)[1:-1]
92 elif self.quad_type == 'LOBATTO':
93 return np.linspace(-1, 1, num=num_nodes)
94 elif self.quad_type == 'RADAU-RIGHT':
95 return np.linspace(-1, 1, num=num_nodes + 1)[1:]
96 elif self.quad_type == 'RADAU-LEFT':
97 return np.linspace(-1, 1, num=num_nodes + 1)[:-1]
99 # Quadrature nodes linked to orthogonal polynomials
100 alpha, beta = self.getTridiagCoefficients(num_nodes)
101 nodes = eigh_tridiagonal(alpha, np.sqrt(beta[1:]))[0]
102 nodes.sort()
104 return nodes
106 def getOrthogPolyCoefficients(self, num_coeff):
107 """
108 Produces a given number of analytic three-term recurrence coefficients.
110 Parameters
111 ----------
112 num_coeff : int
113 Number of coefficients to compute.
115 Returns
116 -------
117 alpha : np.1darray
118 The alpha coefficients of the three-term recurrence.
119 beta : np.1darray
120 The beta coefficients of the three-term recurrence.
121 """
122 if self.node_type == 'LEGENDRE':
123 k = np.arange(num_coeff, dtype=float)
124 alpha = 0 * k
125 beta = k**2 / (4 * k**2 - 1)
126 beta[0] = 2
127 elif self.node_type == 'CHEBY-1':
128 alpha = np.zeros(num_coeff)
129 beta = np.full(num_coeff, 0.25)
130 beta[0] = np.pi
131 if num_coeff > 1:
132 beta[1] = 0.5
133 elif self.node_type == 'CHEBY-2':
134 alpha = np.zeros(num_coeff)
135 beta = np.full(num_coeff, 0.25)
136 beta[0] = np.pi / 2
137 elif self.node_type == 'CHEBY-3':
138 alpha = np.zeros(num_coeff)
139 alpha[0] = 0.5
140 beta = np.full(num_coeff, 0.25)
141 beta[0] = np.pi
142 elif self.node_type == 'CHEBY-4':
143 alpha = np.zeros(num_coeff)
144 alpha[0] = -0.5
145 beta = np.full(num_coeff, 0.25)
146 beta[0] = np.pi
147 return alpha, beta
149 def evalOrthogPoly(self, t, alpha, beta):
150 """
151 Evaluates the two higher order orthogonal polynomials corresponding
152 to the given (alpha,beta) coefficients.
154 Parameters
155 ----------
156 t : float or np.1darray
157 The point where to evaluate the orthogonal polynomials.
158 alpha : np.1darray
159 The alpha coefficients of the three-term recurrence.
160 beta : np.1darray
161 The beta coefficients of the three-term recurrence.
163 Returns
164 -------
165 pi[0] : float or np.1darray
166 The second higher order orthogonal polynomial evaluation.
167 pi[1] : float or np.1darray
168 The higher oder orthogonal polynomial evaluation.
169 """
170 t = np.asarray(t, dtype=float)
171 pi = np.array([np.zeros_like(t) for i in range(3)])
172 pi[1:] += 1
173 for alpha_j, beta_j in zip(alpha, beta):
174 pi[2] *= t - alpha_j
175 pi[0] *= beta_j
176 pi[2] -= pi[0]
177 pi[0] = pi[1]
178 pi[1] = pi[2]
179 return pi[0], pi[1]
181 def getTridiagCoefficients(self, num_nodes):
182 """
183 Computes recurrence coefficients for the tridiagonal Jacobian matrix,
184 taking into account the quadrature type.
186 Parameters
187 ----------
188 num_nodes : int
189 Number of nodes that should be computed from those coefficients.
191 Returns
192 -------
193 alpha : np.1darray
194 The modified alpha coefficients of the three-term recurrence.
195 beta : np.1darray
196 The modified beta coefficients of the three-term recurrence.
197 """
198 # Coefficients for Gauss quadrature type
199 alpha, beta = self.getOrthogPolyCoefficients(num_nodes)
201 # If not Gauss quadrature type, modify the alpha/beta coefficients
202 if self.quad_type.startswith('RADAU'):
203 b = -1.0 if self.quad_type.endswith('LEFT') else 1.0
204 b1, b2 = self.evalOrthogPoly(b, alpha[:-1], beta[:-1])[:2]
205 alpha[-1] = b - beta[-1] * b1 / b2
206 elif self.quad_type == 'LOBATTO':
207 a, b = -1.0, 1.0
208 a2, a1 = self.evalOrthogPoly(a, alpha[:-1], beta[:-1])[:2]
209 b2, b1 = self.evalOrthogPoly(b, alpha[:-1], beta[:-1])[:2]
210 alpha[-1], beta[-1] = np.linalg.solve([[a1, a2], [b1, b2]], [a * a1, b * b1])
211 return alpha, beta