Coverage for chebpy/core/algorithms.py: 99%
194 statements
« prev ^ index » next coverage.py v7.10.2, created at 2025-08-07 10:30 +0000
1"""Numerical algorithms for Chebyshev approximation and manipulation.
3This module provides core numerical algorithms used throughout the ChebPy package,
4including rootfinding, barycentric interpolation, Chebyshev coefficient manipulation,
5and adaptive approximation techniques.
7The algorithms implemented here are based on established numerical methods for
8working with Chebyshev polynomials and approximations, many of which are described
9in Trefethen's "Approximation Theory and Approximation Practice".
10"""
12import warnings
14import numpy as np
15from numpy.fft import fft, ifft
17from .decorators import preandpostprocess
18from .settings import _preferences as prefs
19from .utilities import Interval, infnorm
21# supress numpy division and multiply warnings
22np.seterr(divide="ignore", invalid="ignore")
24# constants
25SPLITPOINT = -0.004849834917525
28# local helpers
29def find(x: np.ndarray) -> np.ndarray:
30 """Find the indices of non-zero elements in an array.
32 A simple wrapper around numpy.where that returns only the indices.
34 Args:
35 x (array-like): Input array.
37 Returns:
38 numpy.ndarray: Indices of non-zero elements in the input array.
39 """
40 return np.where(x)[0]
43def rootsunit(ak: np.ndarray, htol: float = None) -> np.ndarray:
44 """Compute the roots of a function on [-1,1] using Chebyshev coefficients.
46 This function finds the real roots of a function on the interval [-1,1]
47 using the coefficients in its Chebyshev series representation. For large
48 degree polynomials, it uses a recursive subdivision approach.
50 Args:
51 ak (numpy.ndarray): Coefficients of the Chebyshev series.
52 htol (float, optional): Tolerance for determining which roots to keep.
53 Defaults to 100 * machine epsilon.
55 Returns:
56 numpy.ndarray: Array of roots in the interval [-1,1], sorted in
57 ascending order.
59 References:
60 I. J. Good, "The colleague matrix, a Chebyshev analogue of the
61 companion matrix", Quarterly Journal of Mathematics 12 (1961).
63 J. A. Boyd, "Computing zeros on a real interval through
64 Chebyshev expansion and polynomial rootfinding", SIAM Journal on
65 Numerical Analysis 40 (2002).
67 L. N. Trefethen, Approximation Theory and Approximation
68 Practice, SIAM, 2013, chapter 18.
69 """
70 htol = htol if htol is not None else 1e2 * prefs.eps
71 n = standard_chop(ak, tol=htol)
72 ak = ak[:n]
74 # if n > 50, we split and recurse
75 if n > 50:
76 chebpts = chebpts2(ak.size)
77 lmap = Interval(-1, SPLITPOINT)
78 rmap = Interval(SPLITPOINT, 1)
79 lpts = lmap(chebpts)
80 rpts = rmap(chebpts)
81 lval = clenshaw(lpts, ak)
82 rval = clenshaw(rpts, ak)
83 lcfs = vals2coeffs2(lval)
84 rcfs = vals2coeffs2(rval)
85 lrts = rootsunit(lcfs, 2 * htol)
86 rrts = rootsunit(rcfs, 2 * htol)
87 return np.append(lmap(lrts), rmap(rrts))
89 # trivial base case
90 if n <= 1:
91 return np.array([])
93 # nontrivial base case: either compute directly or solve
94 # a Colleague Matrix eigenvalue problem
95 if n == 2:
96 rts = np.array([-ak[0] / ak[1]])
97 elif n <= 50:
98 v = 0.5 * np.ones(n - 2)
99 colleague_matrix = np.diag(v, -1) + np.diag(v, 1)
100 colleague_matrix[0, 1] = 1
101 coeffs_matrix = np.zeros(colleague_matrix.shape, dtype=ak.dtype)
102 coeffs_matrix[-1, :] = ak[:-1]
103 eigenvalue_matrix = colleague_matrix - 0.5 * 1.0 / ak[-1] * coeffs_matrix
104 rts = np.linalg.eigvals(eigenvalue_matrix)
106 # discard values with large imaginary part and treat the remaining
107 # ones as real; then sort and retain only the roots inside [-1,1]
108 mask = abs(np.imag(rts)) < htol
109 rts = np.real(rts[mask])
110 rts = rts[abs(rts) <= 1.0 + htol]
111 rts = np.sort(rts)
112 if rts.size >= 2:
113 rts[0] = max([rts[0], -1])
114 rts[-1] = min([rts[-1], 1])
115 return rts
118@preandpostprocess
119def bary(xx: np.ndarray, fk: np.ndarray, xk: np.ndarray, vk: np.ndarray) -> np.ndarray:
120 """Evaluate a function using the barycentric interpolation formula.
122 This function implements the barycentric interpolation formula for evaluating
123 a function at arbitrary points given its values at a set of nodes. It uses
124 an efficient algorithm that switches between two implementations based on
125 the number of evaluation points.
127 Args:
128 xx (numpy.ndarray): Array of evaluation points.
129 fk (numpy.ndarray): Array of function values at the interpolation nodes xk.
130 xk (numpy.ndarray): Array of interpolation nodes.
131 vk (numpy.ndarray): Barycentric weights corresponding to the interpolation nodes xk.
133 Returns:
134 numpy.ndarray: Function values at the evaluation points xx.
136 References:
137 J.P. Berrut, L.N. Trefethen, Barycentric Lagrange Interpolation, SIAM
138 Review (2004)
139 """
140 # either iterate over the evaluation points, or ...
141 if xx.size < 4 * xk.size:
142 out = np.zeros(xx.size)
143 for i in range(xx.size):
144 tt = vk / (xx[i] - xk)
145 out[i] = np.dot(tt, fk) / tt.sum()
147 # ... iterate over the barycenters
148 else:
149 numer = np.zeros(xx.size)
150 denom = np.zeros(xx.size)
151 for j in range(xk.size):
152 temp = vk[j] / (xx - xk[j])
153 numer = numer + temp * fk[j]
154 denom = denom + temp
155 out = numer / denom
157 # replace NaNs
158 for k in find(np.isnan(out)):
159 idx = find(xx[k] == xk)
160 if idx.size > 0:
161 out[k] = fk[idx[0]]
163 return out
166@preandpostprocess
167def clenshaw(xx: np.ndarray, ak: np.ndarray) -> np.ndarray:
168 """Evaluate a Chebyshev series using Clenshaw's algorithm.
170 This function implements Clenshaw's algorithm for the evaluation of a
171 first-kind Chebyshev series expansion at an array of points.
173 Args:
174 xx (numpy.ndarray): Array of points at which to evaluate the series.
175 ak (numpy.ndarray): Coefficients of the Chebyshev series.
177 Returns:
178 numpy.ndarray: Values of the Chebyshev series at the points xx.
180 References:
181 C. W. Clenshaw, "A note on the summation of Chebyshev series",
182 Mathematics of Computation, Vol. 9, No. 51, 1955, pp. 118-120.
183 """
184 bk1 = 0 * xx
185 bk2 = 0 * xx
186 xx = 2 * xx
187 idx = range(ak.size)
188 for k in idx[ak.size : 1 : -2]:
189 bk2 = ak[k] + xx * bk1 - bk2
190 bk1 = ak[k - 1] + xx * bk2 - bk1
191 if np.mod(ak.size - 1, 2) == 1:
192 bk1, bk2 = ak[1] + xx * bk1 - bk2, bk1
193 out = ak[0] + 0.5 * xx * bk1 - bk2
194 return out
197def standard_chop(coeffs: np.ndarray, tol: float = None) -> int:
198 """Determine where to truncate a Chebyshev series based on coefficient decay.
200 This function determines an appropriate cutoff point for a Chebyshev series
201 by analyzing the decay of its coefficients. It implements the algorithm
202 described by Aurentz and Trefethen.
204 Args:
205 coeffs (numpy.ndarray): Coefficients of the Chebyshev series.
206 tol (float, optional): Tolerance for determining the cutoff point.
207 Defaults to machine epsilon from preferences.
209 Returns:
210 int: Index at which to truncate the series.
212 References:
213 J. Aurentz and L.N. Trefethen, "Chopping a Chebyshev series" (2015)
214 (http://arxiv.org/pdf/1512.01803v1.pdf)
215 """
216 # check magnitude of tol:
217 tol = tol if tol is not None else prefs.eps
218 if tol >= 1: # pragma: no cover
219 cutoff = 1
220 return cutoff
222 # ensure length at least 17:
223 n = coeffs.size
224 cutoff = n
225 if n < 17:
226 return cutoff
228 # Step 1: Convert coeffs input to a new monotonically nonincreasing
229 # vector (envelope) normalized to begin with the value 1.
230 b = np.flipud(np.abs(coeffs))
231 m = np.flipud(np.maximum.accumulate(b))
232 if m[0] == 0.0:
233 # TODO: check this
234 cutoff = 1 # cutoff = 0
235 return cutoff
236 envelope = m / m[0]
238 # Step 2: Scan envelope for a value plateauPoint, the first point, if any,
239 # that is followed by a plateau
240 for j in np.arange(1, n):
241 j2 = round(1.25 * j + 5)
242 if j2 > n - 1:
243 # there is no plateau: exit
244 return cutoff
245 e1 = envelope[j]
246 e2 = envelope[int(j2)]
247 r = 3 * (1 - np.log(e1) / np.log(tol))
248 plateau = (e1 == 0.0) | (e2 / e1 > r)
249 if plateau:
250 # a plateau has been found: go to Step 3
251 plateau_point = j
252 break
254 # Step 3: Fix cutoff at a point where envelope, plus a linear function
255 # included to bias the result towards the left end, is minimal.
256 if envelope[int(plateau_point)] == 0.0:
257 cutoff = plateau_point
258 else:
259 j3 = sum(envelope >= tol ** (7.0 / 6.0))
260 if j3 < j2:
261 j2 = j3 + 1
262 envelope[j2] = tol ** (7.0 / 6.0)
263 cc = np.log10(envelope[: int(j2)])
264 cc = cc + np.linspace(0, (-1.0 / 3.0) * np.log10(tol), int(j2))
265 d = np.argmin(cc)
266 # TODO: check this
267 cutoff = d # + 2
268 return min((cutoff, n - 1))
271def adaptive(cls: type, fun: callable, hscale: float = 1, maxpow2: int = None) -> np.ndarray:
272 """Adaptively determine the number of points needed to represent a function.
274 This function implements an adaptive algorithm to determine the appropriate
275 number of points needed to represent a function to a specified tolerance.
276 It cycles over powers of two, evaluating the function at Chebyshev points
277 and checking if the resulting coefficients can be truncated.
279 Args:
280 cls: The class that provides the _chebpts and _vals2coeffs methods.
281 fun (callable): The function to be approximated.
282 hscale (float, optional): Scale factor for the tolerance. Defaults to 1.
283 maxpow2 (int, optional): Maximum power of 2 to try. If None, uses the
284 value from preferences.
286 Returns:
287 numpy.ndarray: Coefficients of the Chebyshev series representing the function.
289 Warns:
290 UserWarning: If the constructor does not converge within the maximum
291 number of iterations.
292 """
293 minpow2 = 4 # 17 points
294 maxpow2 = maxpow2 if maxpow2 is not None else prefs.maxpow2
295 for k in range(minpow2, max(minpow2, maxpow2) + 1):
296 n = 2**k + 1
297 points = cls._chebpts(n)
298 values = fun(points)
299 coeffs = cls._vals2coeffs(values)
300 eps = prefs.eps
301 tol = eps * max(hscale, 1) # scale (decrease) tolerance by hscale
302 chplen = standard_chop(coeffs, tol=tol)
303 if chplen < coeffs.size:
304 coeffs = coeffs[:chplen]
305 break
306 if k == maxpow2:
307 warnings.warn(f"The {cls.__name__} constructor did not converge: using {n} points")
308 break
309 return coeffs
312def coeffmult(fc: np.ndarray, gc: np.ndarray) -> np.ndarray:
313 """Multiply two Chebyshev series in coefficient space.
315 This function performs multiplication of two Chebyshev series represented by
316 their coefficients. It uses FFT-based convolution for efficiency.
318 Args:
319 fc (numpy.ndarray): Coefficients of the first Chebyshev series.
320 gc (numpy.ndarray): Coefficients of the second Chebyshev series.
322 Returns:
323 numpy.ndarray: Coefficients of the product series.
325 Note:
326 The input series must have the same length.
327 """
328 fc_extended = np.append(2.0 * fc[:1], (fc[1:], fc[:0:-1]))
329 gc_extended = np.append(2.0 * gc[:1], (gc[1:], gc[:0:-1]))
330 ak = ifft(fft(fc_extended) * fft(gc_extended))
331 ak = np.append(ak[:1], ak[1:] + ak[:0:-1]) * 0.25
332 ak = ak[: fc.size]
333 inputcfs = np.append(fc, gc)
334 out = np.real(ak) if np.isreal(inputcfs).all() else ak
335 return out
338def barywts2(n: int) -> np.ndarray:
339 """Compute barycentric weights for Chebyshev points of the second kind.
341 This function calculates the barycentric weights used in the barycentric
342 interpolation formula for Chebyshev points of the second kind.
344 Args:
345 n (int): Number of points (n+1 weights will be computed).
347 Returns:
348 numpy.ndarray: Array of barycentric weights.
350 Note:
351 For Chebyshev points of the second kind, the weights have a simple
352 explicit formula with alternating signs.
353 """
354 if n == 0:
355 wts = np.array([])
356 elif n == 1:
357 wts = np.array([1])
358 else:
359 wts = np.append(np.ones(n - 1), 0.5)
360 wts[n - 2 :: -2] = -1
361 wts[0] = 0.5 * wts[0]
362 return wts
365def chebpts2(n: int) -> np.ndarray:
366 """Compute Chebyshev points of the second kind.
368 This function calculates the n Chebyshev points of the second kind in the
369 interval [-1, 1], which are the extrema of the Chebyshev polynomial T_{n-1}
370 together with the endpoints ±1.
372 Args:
373 n (int): Number of points to compute.
375 Returns:
376 numpy.ndarray: Array of n Chebyshev points of the second kind.
378 Note:
379 The points are ordered from left to right on the interval [-1, 1].
380 """
381 if n == 1:
382 pts = np.array([0.0])
383 else:
384 nn = np.arange(n)
385 pts = np.cos(nn[::-1] * np.pi / (n - 1))
386 return pts
389def vals2coeffs2(vals: np.ndarray) -> np.ndarray:
390 """Convert function values to Chebyshev coefficients.
392 This function maps function values at Chebyshev points of the second kind
393 to coefficients of the corresponding first-kind Chebyshev polynomial expansion.
394 It uses an FFT-based algorithm for efficiency.
396 Args:
397 vals (numpy.ndarray): Function values at Chebyshev points of the second kind.
399 Returns:
400 numpy.ndarray: Coefficients of the first-kind Chebyshev polynomial expansion.
402 Note:
403 This transformation is the discrete cosine transform of type I (DCT-I),
404 which is implemented here using FFT for efficiency.
405 """
406 n = vals.size
407 if n <= 1:
408 coeffs = vals
409 return coeffs
410 tmp = np.append(vals[::-1], vals[1:-1])
411 if np.isreal(vals).all():
412 coeffs = ifft(tmp)
413 coeffs = np.real(coeffs)
414 elif np.isreal(1j * vals).all(): # pragma: no cover
415 coeffs = ifft(np.imag(tmp))
416 coeffs = 1j * np.real(coeffs)
417 else:
418 coeffs = ifft(tmp)
419 coeffs = coeffs[:n]
420 coeffs[1 : n - 1] = 2 * coeffs[1 : n - 1]
421 return coeffs
424def coeffs2vals2(coeffs: np.ndarray) -> np.ndarray:
425 """Convert Chebyshev coefficients to function values.
427 This function maps coefficients of a first-kind Chebyshev polynomial expansion
428 to function values at Chebyshev points of the second kind. It uses an FFT-based
429 algorithm for efficiency.
431 Args:
432 coeffs (numpy.ndarray): Coefficients of the first-kind Chebyshev polynomial expansion.
434 Returns:
435 numpy.ndarray: Function values at Chebyshev points of the second kind.
437 Note:
438 This transformation is the inverse discrete cosine transform of type I (IDCT-I),
439 which is implemented here using FFT for efficiency. It is the inverse of vals2coeffs2.
440 """
441 n = coeffs.size
442 if n <= 1:
443 vals = coeffs
444 return vals
445 coeffs = coeffs.copy()
446 coeffs[1 : n - 1] = 0.5 * coeffs[1 : n - 1]
447 tmp = np.append(coeffs, coeffs[n - 2 : 0 : -1])
448 if np.isreal(coeffs).all():
449 vals = fft(tmp)
450 vals = np.real(vals)
451 elif np.isreal(1j * coeffs).all(): # pragma: no cover
452 vals = fft(np.imag(tmp))
453 vals = 1j * np.real(vals)
454 else:
455 vals = fft(tmp)
456 vals = vals[n - 1 :: -1]
457 return vals
460def newtonroots(fun: callable, rts: np.ndarray, tol: float = None, maxiter: int = None) -> np.ndarray:
461 """Refine root approximations using Newton's method.
463 This function applies Newton's method to refine the approximations of roots
464 for a callable and differentiable function. It is typically used to polish
465 already computed roots to higher accuracy.
467 Args:
468 fun (callable): A callable and differentiable function.
469 rts (numpy.ndarray): Initial approximations of the roots.
470 tol (float, optional): Tolerance for convergence. Defaults to 2 * machine epsilon.
471 maxiter (int, optional): Maximum number of iterations. Defaults to value from preferences.
473 Returns:
474 numpy.ndarray: Refined approximations of the roots.
476 Note:
477 The function must support differentiation via a .diff() method that returns
478 the derivative function.
479 """
480 tol = tol if tol is not None else 2 * prefs.eps
481 maxiter = maxiter if maxiter is not None else prefs.maxiter
482 if rts.size > 0:
483 dfun = fun.diff()
484 prv = np.inf * rts
485 count = 0
486 while (infnorm(rts - prv) > tol) & (count <= maxiter):
487 count += 1
488 prv = rts
489 rts = rts - fun(rts) / dfun(rts)
490 return rts