Coverage for src / chebpy / algorithms.py: 99%
194 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-20 05:55 +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.
184 Examples:
185 >>> import numpy as np
186 >>> coeffs = np.array([1.0])
187 >>> x = np.array([0.0])
188 >>> result = clenshaw(x, coeffs)
189 >>> bool(abs(float(result[0]) - 1.0) < 1e-10)
190 True
191 """
192 bk1 = 0 * xx
193 bk2 = 0 * xx
194 xx = 2 * xx
195 idx = range(ak.size)
196 for k in idx[ak.size : 1 : -2]:
197 bk2 = ak[k] + xx * bk1 - bk2
198 bk1 = ak[k - 1] + xx * bk2 - bk1
199 if np.mod(ak.size - 1, 2) == 1:
200 bk1, bk2 = ak[1] + xx * bk1 - bk2, bk1
201 out = ak[0] + 0.5 * xx * bk1 - bk2
202 return out
205def standard_chop(coeffs: np.ndarray, tol: float = None) -> int:
206 """Determine where to truncate a Chebyshev series based on coefficient decay.
208 This function determines an appropriate cutoff point for a Chebyshev series
209 by analyzing the decay of its coefficients. It implements the algorithm
210 described by Aurentz and Trefethen.
212 Args:
213 coeffs (numpy.ndarray): Coefficients of the Chebyshev series.
214 tol (float, optional): Tolerance for determining the cutoff point.
215 Defaults to machine epsilon from preferences.
217 Returns:
218 int: Index at which to truncate the series.
220 References:
221 J. Aurentz and L.N. Trefethen, "Chopping a Chebyshev series" (2015)
222 (http://arxiv.org/pdf/1512.01803v1.pdf)
223 """
224 # check magnitude of tol:
225 tol = tol if tol is not None else prefs.eps
226 if tol >= 1: # pragma: no cover
227 cutoff = 1
228 return cutoff
230 # ensure length at least 17:
231 n = coeffs.size
232 cutoff = n
233 if n < 17:
234 return cutoff
236 # Step 1: Convert coeffs input to a new monotonically nonincreasing
237 # vector (envelope) normalized to begin with the value 1.
238 b = np.flipud(np.abs(coeffs))
239 m = np.flipud(np.maximum.accumulate(b))
240 if m[0] == 0.0:
241 # TODO: check this
242 cutoff = 1 # cutoff = 0
243 return cutoff
244 envelope = m / m[0]
246 # Step 2: Scan envelope for a value plateauPoint, the first point, if any,
247 # that is followed by a plateau
248 for j in np.arange(1, n):
249 j2 = round(1.25 * j + 5)
250 if j2 > n - 1:
251 # there is no plateau: exit
252 return cutoff
253 e1 = envelope[j]
254 e2 = envelope[int(j2)]
255 r = 3 * (1 - np.log(e1) / np.log(tol))
256 plateau = (e1 == 0.0) | (e2 / e1 > r)
257 if plateau:
258 # a plateau has been found: go to Step 3
259 plateau_point = j
260 break
262 # Step 3: Fix cutoff at a point where envelope, plus a linear function
263 # included to bias the result towards the left end, is minimal.
264 if envelope[int(plateau_point)] == 0.0:
265 cutoff = plateau_point
266 else:
267 j3 = sum(envelope >= tol ** (7.0 / 6.0))
268 if j3 < j2:
269 j2 = j3 + 1
270 envelope[j2] = tol ** (7.0 / 6.0)
271 cc = np.log10(envelope[: int(j2)])
272 cc = cc + np.linspace(0, (-1.0 / 3.0) * np.log10(tol), int(j2))
273 d = np.argmin(cc)
274 # TODO: check this
275 cutoff = d # + 2
276 return min((cutoff, n - 1))
279def adaptive(cls: type, fun: callable, hscale: float = 1, maxpow2: int = None) -> np.ndarray:
280 """Adaptively determine the number of points needed to represent a function.
282 This function implements an adaptive algorithm to determine the appropriate
283 number of points needed to represent a function to a specified tolerance.
284 It cycles over powers of two, evaluating the function at Chebyshev points
285 and checking if the resulting coefficients can be truncated.
287 Args:
288 cls: The class that provides the _chebpts and _vals2coeffs methods.
289 fun (callable): The function to be approximated.
290 hscale (float, optional): Scale factor for the tolerance. Defaults to 1.
291 maxpow2 (int, optional): Maximum power of 2 to try. If None, uses the
292 value from preferences.
294 Returns:
295 numpy.ndarray: Coefficients of the Chebyshev series representing the function.
297 Warns:
298 UserWarning: If the constructor does not converge within the maximum
299 number of iterations.
300 """
301 minpow2 = 4 # 17 points
302 maxpow2 = maxpow2 if maxpow2 is not None else prefs.maxpow2
303 for k in range(minpow2, max(minpow2, maxpow2) + 1):
304 n = 2**k + 1
305 points = cls._chebpts(n)
306 values = fun(points)
307 coeffs = cls._vals2coeffs(values)
308 eps = prefs.eps
309 tol = eps * max(hscale, 1) # scale (decrease) tolerance by hscale
310 chplen = standard_chop(coeffs, tol=tol)
311 if chplen < coeffs.size:
312 coeffs = coeffs[:chplen]
313 break
314 if k == maxpow2:
315 warnings.warn(f"The {cls.__name__} constructor did not converge: using {n} points")
316 break
317 return coeffs
320def coeffmult(fc: np.ndarray, gc: np.ndarray) -> np.ndarray:
321 """Multiply two Chebyshev series in coefficient space.
323 This function performs multiplication of two Chebyshev series represented by
324 their coefficients. It uses FFT-based convolution for efficiency.
326 Args:
327 fc (numpy.ndarray): Coefficients of the first Chebyshev series.
328 gc (numpy.ndarray): Coefficients of the second Chebyshev series.
330 Returns:
331 numpy.ndarray: Coefficients of the product series.
333 Note:
334 The input series must have the same length.
335 """
336 fc_extended = np.append(2.0 * fc[:1], (fc[1:], fc[:0:-1]))
337 gc_extended = np.append(2.0 * gc[:1], (gc[1:], gc[:0:-1]))
338 ak = ifft(fft(fc_extended) * fft(gc_extended))
339 ak = np.append(ak[:1], ak[1:] + ak[:0:-1]) * 0.25
340 ak = ak[: fc.size]
341 inputcfs = np.append(fc, gc)
342 out = np.real(ak) if np.isreal(inputcfs).all() else ak
343 return out
346def barywts2(n: int) -> np.ndarray:
347 """Compute barycentric weights for Chebyshev points of the second kind.
349 This function calculates the barycentric weights used in the barycentric
350 interpolation formula for Chebyshev points of the second kind.
352 Args:
353 n (int): Number of points (n+1 weights will be computed).
355 Returns:
356 numpy.ndarray: Array of barycentric weights.
358 Note:
359 For Chebyshev points of the second kind, the weights have a simple
360 explicit formula with alternating signs.
361 """
362 if n == 0:
363 wts = np.array([])
364 elif n == 1:
365 wts = np.array([1])
366 else:
367 wts = np.append(np.ones(n - 1), 0.5)
368 wts[n - 2 :: -2] = -1
369 wts[0] = 0.5 * wts[0]
370 return wts
373def chebpts2(n: int) -> np.ndarray:
374 """Compute Chebyshev points of the second kind.
376 This function calculates the n Chebyshev points of the second kind in the
377 interval [-1, 1], which are the extrema of the Chebyshev polynomial T_{n-1}
378 together with the endpoints ±1.
380 Args:
381 n (int): Number of points to compute.
383 Returns:
384 numpy.ndarray: Array of n Chebyshev points of the second kind.
386 Note:
387 The points are ordered from left to right on the interval [-1, 1].
388 """
389 if n == 1:
390 pts = np.array([0.0])
391 else:
392 nn = np.arange(n)
393 pts = np.cos(nn[::-1] * np.pi / (n - 1))
394 return pts
397def vals2coeffs2(vals: np.ndarray) -> np.ndarray:
398 """Convert function values to Chebyshev coefficients.
400 This function maps function values at Chebyshev points of the second kind
401 to coefficients of the corresponding first-kind Chebyshev polynomial expansion.
402 It uses an FFT-based algorithm for efficiency.
404 Args:
405 vals (numpy.ndarray): Function values at Chebyshev points of the second kind.
407 Returns:
408 numpy.ndarray: Coefficients of the first-kind Chebyshev polynomial expansion.
410 Note:
411 This transformation is the discrete cosine transform of type I (DCT-I),
412 which is implemented here using FFT for efficiency.
413 """
414 n = vals.size
415 if n <= 1:
416 coeffs = vals
417 return coeffs
418 tmp = np.append(vals[::-1], vals[1:-1])
419 if np.isreal(vals).all():
420 coeffs = ifft(tmp)
421 coeffs = np.real(coeffs)
422 elif np.isreal(1j * vals).all(): # pragma: no cover
423 coeffs = ifft(np.imag(tmp))
424 coeffs = 1j * np.real(coeffs)
425 else:
426 coeffs = ifft(tmp)
427 coeffs = coeffs[:n]
428 coeffs[1 : n - 1] = 2 * coeffs[1 : n - 1]
429 return coeffs
432def coeffs2vals2(coeffs: np.ndarray) -> np.ndarray:
433 """Convert Chebyshev coefficients to function values.
435 This function maps coefficients of a first-kind Chebyshev polynomial expansion
436 to function values at Chebyshev points of the second kind. It uses an FFT-based
437 algorithm for efficiency.
439 Args:
440 coeffs (numpy.ndarray): Coefficients of the first-kind Chebyshev polynomial expansion.
442 Returns:
443 numpy.ndarray: Function values at Chebyshev points of the second kind.
445 Note:
446 This transformation is the inverse discrete cosine transform of type I (IDCT-I),
447 which is implemented here using FFT for efficiency. It is the inverse of vals2coeffs2.
448 """
449 n = coeffs.size
450 if n <= 1:
451 vals = coeffs
452 return vals
453 coeffs = coeffs.copy()
454 coeffs[1 : n - 1] = 0.5 * coeffs[1 : n - 1]
455 tmp = np.append(coeffs, coeffs[n - 2 : 0 : -1])
456 if np.isreal(coeffs).all():
457 vals = fft(tmp)
458 vals = np.real(vals)
459 elif np.isreal(1j * coeffs).all(): # pragma: no cover
460 vals = fft(np.imag(tmp))
461 vals = 1j * np.real(vals)
462 else:
463 vals = fft(tmp)
464 vals = vals[n - 1 :: -1]
465 return vals
468def newtonroots(fun: callable, rts: np.ndarray, tol: float = None, maxiter: int = None) -> np.ndarray:
469 """Refine root approximations using Newton's method.
471 This function applies Newton's method to refine the approximations of roots
472 for a callable and differentiable function. It is typically used to polish
473 already computed roots to higher accuracy.
475 Args:
476 fun (callable): A callable and differentiable function.
477 rts (numpy.ndarray): Initial approximations of the roots.
478 tol (float, optional): Tolerance for convergence. Defaults to 2 * machine epsilon.
479 maxiter (int, optional): Maximum number of iterations. Defaults to value from preferences.
481 Returns:
482 numpy.ndarray: Refined approximations of the roots.
484 Note:
485 The function must support differentiation via a .diff() method that returns
486 the derivative function.
487 """
488 tol = tol if tol is not None else 2 * prefs.eps
489 maxiter = maxiter if maxiter is not None else prefs.maxiter
490 if rts.size > 0:
491 dfun = fun.diff()
492 prv = np.inf * rts
493 count = 0
494 while (infnorm(rts - prv) > tol) & (count <= maxiter):
495 count += 1
496 prv = rts
497 rts = rts - fun(rts) / dfun(rts)
498 return rts