Coverage for src/chebpy/algorithms.py: 100%

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

13from collections.abc import Callable

14from typing import Any

16import numpy as np

17from numpy.fft import fft, ifft

19from .decorators import preandpostprocess

20from .settings import _preferences as prefs

21from .utilities import Interval, infnorm

23# supress numpy division and multiply warnings

24np.seterr(divide="ignore", invalid="ignore")

26# constants

27SPLITPOINT = -0.004849834917525

30# local helpers

31def find(x: np.ndarray) -> np.ndarray:

32 """Find the indices of non-zero elements in an array.

34 A simple wrapper around numpy.where that returns only the indices.

36 Args:

37 x (array-like): Input array.

39 Returns:

40 numpy.ndarray: Indices of non-zero elements in the input array.

41 """

42 return np.where(x)[0]

45def rootsunit(ak: np.ndarray, htol: float | None = None) -> np.ndarray:

46 """Compute the roots of a function on [-1,1] using Chebyshev coefficients.

48 This function finds the real roots of a function on the interval [-1,1]

49 using the coefficients in its Chebyshev series representation. For large

50 degree polynomials, it uses a recursive subdivision approach.

52 Args:

53 ak (numpy.ndarray): Coefficients of the Chebyshev series.

54 htol (float, optional): Tolerance for determining which roots to keep.

55 Defaults to 100 * machine epsilon.

57 Returns:

58 numpy.ndarray: Array of roots in the interval [-1,1], sorted in

59 ascending order.

61 References:

62 I. J. Good, "The colleague matrix, a Chebyshev analogue of the

63 companion matrix", Quarterly Journal of Mathematics 12 (1961).

65 J. A. Boyd, "Computing zeros on a real interval through

66 Chebyshev expansion and polynomial rootfinding", SIAM Journal on

67 Numerical Analysis 40 (2002).

69 L. N. Trefethen, Approximation Theory and Approximation

70 Practice, SIAM, 2013, chapter 18.

71 """

72 htol = htol if htol is not None else 1e2 * prefs.eps

73 n = standard_chop(ak, tol=htol)

74 ak = ak[:n]

76 # if n > 50, we split and recurse

77 if n > 50:

78 chebpts = chebpts2(ak.size)

79 lmap = Interval(-1, SPLITPOINT)

80 rmap = Interval(SPLITPOINT, 1)

81 lpts = lmap(chebpts)

82 rpts = rmap(chebpts)

83 lval = clenshaw(lpts, ak)

84 rval = clenshaw(rpts, ak)

85 lcfs = vals2coeffs2(lval)

86 rcfs = vals2coeffs2(rval)

87 lrts = rootsunit(lcfs, 2 * htol)

88 rrts = rootsunit(rcfs, 2 * htol)

89 return np.append(lmap(lrts), rmap(rrts))

91 # trivial base case

92 if n <= 1:

93 return np.array([])

95 # nontrivial base case: either compute directly or solve

96 # a Colleague Matrix eigenvalue problem

97 if n == 2:

98 rts = np.array([-ak[0] / ak[1]])

99 elif n <= 50:

100 v = 0.5 * np.ones(n - 2)

101 colleague_matrix = np.diag(v, -1) + np.diag(v, 1)

102 colleague_matrix[0, 1] = 1

103 coeffs_matrix = np.zeros(colleague_matrix.shape, dtype=ak.dtype)

104 coeffs_matrix[-1, :] = ak[:-1]

105 eigenvalue_matrix = colleague_matrix - 0.5 * 1.0 / ak[-1] * coeffs_matrix

106 rts = np.linalg.eigvals(eigenvalue_matrix)

107

108 # discard values with large imaginary part and treat the remaining

109 # ones as real; then sort and retain only the roots inside [-1,1]

110 mask = abs(np.imag(rts)) < htol

111 rts = np.real(rts[mask])

112 rts = rts[abs(rts) <= 1.0 + htol]

113 rts = np.sort(rts)

114 if rts.size >= 2:

115 rts[0] = max([rts[0], -1])

116 rts[-1] = min([rts[-1], 1])

117 return rts

118

119

120@preandpostprocess

121def bary(xx: np.ndarray, fk: np.ndarray, xk: np.ndarray, vk: np.ndarray) -> np.ndarray:

122 """Evaluate a function using the barycentric interpolation formula.

123

124 This function implements the barycentric interpolation formula for evaluating

125 a function at arbitrary points given its values at a set of nodes. It uses

126 an efficient algorithm that switches between two implementations based on

127 the number of evaluation points.

128

129 Args:

130 xx (numpy.ndarray): Array of evaluation points.

131 fk (numpy.ndarray): Array of function values at the interpolation nodes xk.

132 xk (numpy.ndarray): Array of interpolation nodes.

133 vk (numpy.ndarray): Barycentric weights corresponding to the interpolation nodes xk.

134

135 Returns:

136 numpy.ndarray: Function values at the evaluation points xx.

137

138 References:

139 J.P. Berrut, L.N. Trefethen, Barycentric Lagrange Interpolation, SIAM

140 Review (2004)

141 """

142 # either iterate over the evaluation points, or ...

143 if xx.size < 4 * xk.size:

144 out = np.zeros(xx.size)

145 for i in range(xx.size):

146 tt = vk / (xx[i] - xk)

147 out[i] = np.dot(tt, fk) / tt.sum()

148

149 # ... iterate over the barycenters

150 else:

151 numer = np.zeros(xx.size)

152 denom = np.zeros(xx.size)

153 for j in range(xk.size):

154 temp = vk[j] / (xx - xk[j])

155 numer = numer + temp * fk[j]

156 denom = denom + temp

157 out = numer / denom

158

159 # replace NaNs

160 for k in find(np.isnan(out)):

161 idx = find(xx[k] == xk)

162 if idx.size > 0:

163 out[k] = fk[idx[0]]

164

165 return out

166

167

168@preandpostprocess

169def clenshaw(xx: np.ndarray, ak: np.ndarray) -> np.ndarray:

170 """Evaluate a Chebyshev series using Clenshaw's algorithm.

171

172 This function implements Clenshaw's algorithm for the evaluation of a

173 first-kind Chebyshev series expansion at an array of points.

174

175 Args:

176 xx (numpy.ndarray): Array of points at which to evaluate the series.

177 ak (numpy.ndarray): Coefficients of the Chebyshev series.

178

179 Returns:

180 numpy.ndarray: Values of the Chebyshev series at the points xx.

181

182 References:

183 C. W. Clenshaw, "A note on the summation of Chebyshev series",

184 Mathematics of Computation, Vol. 9, No. 51, 1955, pp. 118-120.

185

186 Examples:

187 >>> import numpy as np

188 >>> coeffs = np.array([1.0])

189 >>> x = np.array([0.0])

190 >>> result = clenshaw(x, coeffs)

191 >>> bool(abs(float(result[0]) - 1.0) < 1e-10)

192 True

193 """

194 bk1 = 0 * xx

195 bk2 = 0 * xx

196 xx = 2 * xx

197 idx = range(ak.size)

198 for k in idx[ak.size : 1 : -2]:

199 bk2 = ak[k] + xx * bk1 - bk2

200 bk1 = ak[k - 1] + xx * bk2 - bk1

201 if np.mod(ak.size - 1, 2) == 1:

202 bk1, bk2 = ak[1] + xx * bk1 - bk2, bk1

203 out: np.ndarray = ak[0] + 0.5 * xx * bk1 - bk2

204 return out

205

206

207def standard_chop(coeffs: np.ndarray, tol: float | None = None) -> int:

208 """Determine where to truncate a Chebyshev series based on coefficient decay.

209

210 This function determines an appropriate cutoff point for a Chebyshev series

211 by analyzing the decay of its coefficients. It implements the algorithm

212 described by Aurentz and Trefethen.

213

214 Args:

215 coeffs (numpy.ndarray): Coefficients of the Chebyshev series.

216 tol (float, optional): Tolerance for determining the cutoff point.

217 Defaults to machine epsilon from preferences.

218

219 Returns:

220 int: Index at which to truncate the series.

221

222 References:

223 J. Aurentz and L.N. Trefethen, "Chopping a Chebyshev series" (2015)

224 (http://arxiv.org/pdf/1512.01803v1.pdf)

225 """

226 # check magnitude of tol:

227 tol = tol if tol is not None else prefs.eps

228 if tol >= 1: # pragma: no cover

229 cutoff = 1

230 return cutoff

231

232 # ensure length at least 17:

233 n = coeffs.size

234 cutoff = n

235 if n < 17:

236 return cutoff

237

238 # Step 1: Convert coeffs input to a new monotonically nonincreasing

239 # vector (envelope) normalized to begin with the value 1.

240 b = np.flipud(np.abs(coeffs))

241 m = np.flipud(np.maximum.accumulate(b))

242 if m[0] == 0.0:

243 cutoff = 1

244 return cutoff

245 envelope = m / m[0]

246

247 # Step 2: Scan envelope for a value plateauPoint, the first point J-1,

248 # if any, that is followed by a plateau. Uses 1-based j to match the

249 # MATLAB reference implementation; envelope is indexed with [j-1].

250 for j in range(2, n + 1):

251 j2 = round(1.25 * j + 5)

252 if j2 > n:

253 # there is no plateau: exit

254 return cutoff

255 e1 = envelope[j - 1]

256 e2 = envelope[int(j2) - 1]

257 r = 3 * (1 - np.log(e1) / np.log(tol))

258 plateau = (e1 == 0.0) | (e2 / e1 > r)

259 if plateau:

260 # a plateau has been found: go to Step 3

261 plateau_point = j - 1

262 break

263

264 # Step 3: Fix cutoff at a point where envelope, plus a linear function

265 # included to bias the result towards the left end, is minimal.

266 if envelope[plateau_point - 1] == 0.0: # pragma: no cover

267 cutoff = plateau_point

268 else:

269 j3 = int(np.sum(envelope >= tol ** (7.0 / 6.0)))

270 if j3 < j2:

271 j2 = j3 + 1

272 envelope[int(j2) - 1] = tol ** (7.0 / 6.0)

273 cc = np.log10(envelope[: int(j2)])

274 cc = cc + np.linspace(0, (-1.0 / 3.0) * np.log10(tol), int(j2))

275 d = np.argmin(cc)

276 cutoff = max(int(d), 1)

277 return cutoff

278

279

280def adaptive(cls: Any, fun: Callable[..., Any], hscale: float = 1, maxpow2: int | None = None) -> np.ndarray:

281 """Adaptively determine the number of points needed to represent a function.

282

283 This function implements an adaptive algorithm to determine the appropriate

284 number of points needed to represent a function to a specified tolerance.

285 It cycles over powers of two, evaluating the function at Chebyshev points

286 and checking if the resulting coefficients can be truncated.

287

288 Args:

289 cls: The class that provides the _chebpts and _vals2coeffs methods.

290 fun (callable): The function to be approximated.

291 hscale (float, optional): Scale factor for the tolerance. Defaults to 1.

292 maxpow2 (int, optional): Maximum power of 2 to try. If None, uses the

293 value from preferences.

294

295 Returns:

296 numpy.ndarray: Coefficients of the Chebyshev series representing the function.

297

298 Warns:

299 UserWarning: If the constructor does not converge within the maximum

300 number of iterations.

301 """

302 minpow2 = 4 # 17 points

303 maxpow2 = maxpow2 if maxpow2 is not None else prefs.maxpow2

304 tol = prefs.eps * max(hscale, 1)

305 coeffs: np.ndarray = np.array([])

306 for k in range(minpow2, max(minpow2, maxpow2) + 1):

307 n = 2**k + 1

308 points = cls._chebpts(n)

309 values = fun(points)

310 coeffs = cls._vals2coeffs(values)

311 # If function values are at or below tolerance the function is

312 # indistinguishable from zero (cf. classicCheck.m vscale==0 guard).

313 vscale = np.max(np.abs(values))

314 if vscale <= tol:

315 coeffs = np.array([0.0])

316 break

317 chplen = standard_chop(coeffs, tol=tol)

318 if chplen < coeffs.size:

319 coeffs = coeffs[:chplen]

320 break

321 if k == maxpow2:

322 warnings.warn(f"The {cls.__name__} constructor did not converge: using {n} points", stacklevel=2)

323 break

324 return coeffs

325

326

327def coeffmult(fc: np.ndarray, gc: np.ndarray) -> np.ndarray:

328 """Multiply two Chebyshev series in coefficient space.

329

330 This function performs multiplication of two Chebyshev series represented by

331 their coefficients. It uses FFT-based convolution for efficiency.

332

333 Args:

334 fc (numpy.ndarray): Coefficients of the first Chebyshev series.

335 gc (numpy.ndarray): Coefficients of the second Chebyshev series.

336

337 Returns:

338 numpy.ndarray: Coefficients of the product series.

339

340 Note:

341 The input series must have the same length.

342 """

343 fc_extended = np.append(2.0 * fc[:1], (fc[1:], fc[:0:-1]))

344 gc_extended = np.append(2.0 * gc[:1], (gc[1:], gc[:0:-1]))

345 ak = ifft(fft(fc_extended) * fft(gc_extended))

346 ak = np.append(ak[:1], ak[1:] + ak[:0:-1]) * 0.25

347 ak = ak[: fc.size]

348 inputcfs = np.append(fc, gc)

349 out = np.real(ak) if np.isreal(inputcfs).all() else ak

350 return out

351

352

353def barywts2(n: int) -> np.ndarray:

354 """Compute barycentric weights for Chebyshev points of the second kind.

355

356 This function calculates the barycentric weights used in the barycentric

357 interpolation formula for Chebyshev points of the second kind.

358

359 Args:

360 n (int): Number of points (n+1 weights will be computed).

361

362 Returns:

363 numpy.ndarray: Array of barycentric weights.

364

365 Note:

366 For Chebyshev points of the second kind, the weights have a simple

367 explicit formula with alternating signs.

368 """

369 if n == 0:

370 wts = np.array([])

371 elif n == 1:

372 wts = np.array([1])

373 else:

374 wts = np.append(np.ones(n - 1), 0.5)

375 wts[n - 2 :: -2] = -1

376 wts[0] = 0.5 * wts[0]

377 return wts

378

379

380def chebpts2(n: int) -> np.ndarray:

381 """Compute Chebyshev points of the second kind.

382

383 This function calculates the n Chebyshev points of the second kind in the

384 interval [-1, 1], which are the extrema of the Chebyshev polynomial T_{n-1}

385 together with the endpoints ±1.

386

387 Args:

388 n (int): Number of points to compute.

389

390 Returns:

391 numpy.ndarray: Array of n Chebyshev points of the second kind.

392

393 Note:

394 The points are ordered from left to right on the interval [-1, 1].

395 """

396 if n == 1:

397 pts = np.array([0.0])

398 else:

399 nn = np.arange(n)

400 pts = np.cos(nn[::-1] * np.pi / (n - 1))

401 return pts

402

403

404def vals2coeffs2(vals: np.ndarray) -> np.ndarray:

405 """Convert function values to Chebyshev coefficients.

406

407 This function maps function values at Chebyshev points of the second kind

408 to coefficients of the corresponding first-kind Chebyshev polynomial expansion.

409 It uses an FFT-based algorithm for efficiency.

410

411 Args:

412 vals (numpy.ndarray): Function values at Chebyshev points of the second kind.

413

414 Returns:

415 numpy.ndarray: Coefficients of the first-kind Chebyshev polynomial expansion.

416

417 Note:

418 This transformation is the discrete cosine transform of type I (DCT-I),

419 which is implemented here using FFT for efficiency.

420 """

421 n = vals.size

422 if n <= 1:

423 coeffs = vals

424 return coeffs

425 tmp = np.append(vals[::-1], vals[1:-1])

426 if np.isreal(vals).all():

427 coeffs = ifft(tmp)

428 coeffs = np.real(coeffs)

429 elif np.isreal(1j * vals).all(): # pragma: no cover

430 coeffs = ifft(np.imag(tmp))

431 coeffs = 1j * np.real(coeffs)

432 else:

433 coeffs = ifft(tmp)

434 coeffs = coeffs[:n]

435 coeffs[1 : n - 1] = 2 * coeffs[1 : n - 1]

436 return coeffs

437

438

439def coeffs2vals2(coeffs: np.ndarray) -> np.ndarray:

440 """Convert Chebyshev coefficients to function values.

441

442 This function maps coefficients of a first-kind Chebyshev polynomial expansion

443 to function values at Chebyshev points of the second kind. It uses an FFT-based

444 algorithm for efficiency.

445

446 Args:

447 coeffs (numpy.ndarray): Coefficients of the first-kind Chebyshev polynomial expansion.

448

449 Returns:

450 numpy.ndarray: Function values at Chebyshev points of the second kind.

451

452 Note:

453 This transformation is the inverse discrete cosine transform of type I (IDCT-I),

454 which is implemented here using FFT for efficiency. It is the inverse of vals2coeffs2.

455 """

456 n = coeffs.size

457 if n <= 1:

458 vals = coeffs

459 return vals

460 coeffs = coeffs.copy()

461 coeffs[1 : n - 1] = 0.5 * coeffs[1 : n - 1]

462 tmp = np.append(coeffs, coeffs[n - 2 : 0 : -1])

463 if np.isreal(coeffs).all():

464 vals = fft(tmp)

465 vals = np.real(vals)

466 elif np.isreal(1j * coeffs).all(): # pragma: no cover

467 vals = fft(np.imag(tmp))

468 vals = 1j * np.real(vals)

469 else:

470 vals = fft(tmp)

471 vals = vals[n - 1 :: -1]

472 return vals

473

474

475def cheb2leg(c: np.ndarray) -> np.ndarray:

476 """Convert Chebyshev coefficients to Legendre coefficients.

477

478 Converts the vector ``c`` of Chebyshev coefficients to a vector of Legendre

479 coefficients such that::

480

481 c[0]*T_0 + c[1]*T_1 + ... = l[0]*P_0 + l[1]*P_1 + ...

482

483 Uses a stable O(n²) three-term recurrence derived from the Chebyshev

484 recurrence ``T_n = 2x T_{n-1} - T_{n-2}``.

485

486 Args:

487 c (array-like): Chebyshev coefficients.

488

489 Returns:

490 numpy.ndarray: Legendre coefficients of the same polynomial.

491 """

492 c = np.asarray(c, dtype=float)

493 n = c.size

494 if n <= 1:

495 return c.copy()

496

497 # Build Legendre coefficients via the recurrence:

498 # M[j, col] = coeff of P_j in T_col

499 # Recurrence: M[j,col] = 2j/(2j-1)*M[j-1,col-1]

500 # + 2(j+1)/(2j+3)*M[j+1,col-1]

501 # - M[j,col-2]

502 # Initial columns: M[:,0] = [1,0,...], M[:,1] = [0,1,0,...]

503 leg_coeffs = np.zeros(n)

504

505 prev_prev = np.zeros(n)

506 prev_prev[0] = 1.0 # T_0 = P_0

507 leg_coeffs += c[0] * prev_prev

508

509 prev = np.zeros(n)

510 prev[1] = 1.0 # T_1 = P_1

511 leg_coeffs += c[1] * prev

512

513 j = np.arange(n)

514 for col in range(2, n):

515 curr = np.zeros(n)

516 # 2j/(2j-1) * prev[j-1] (for j >= 1)

517 curr[1:] += 2.0 * j[1:] / (2.0 * j[1:] - 1.0) * prev[:-1]

518 # 2(j+1)/(2j+3) * prev[j+1] (for j+1 <= n-1)

519 curr[:-1] += 2.0 * (j[:-1] + 1.0) / (2.0 * j[:-1] + 3.0) * prev[1:]

520 curr -= prev_prev

521 leg_coeffs += c[col] * curr

522 prev_prev = prev

523 prev = curr

524

525 return leg_coeffs

526

527

528def leg2cheb(c: np.ndarray) -> np.ndarray:

529 """Convert Legendre coefficients to Chebyshev coefficients.

530

531 Converts the vector ``c`` of Legendre coefficients to a vector of Chebyshev

532 coefficients such that::

533

534 c[0]*P_0 + c[1]*P_1 + ... = l[0]*T_0 + l[1]*T_1 + ...

535

536 Uses a stable O(n²) three-term recurrence derived from the Legendre

537 recurrence ``(n+1) P_{n+1} = (2n+1) x P_n - n P_{n-1}``.

538

539 Args:

540 c (array-like): Legendre coefficients.

541

542 Returns:

543 numpy.ndarray: Chebyshev coefficients of the same polynomial.

544 """

545 c = np.asarray(c, dtype=float)

546 n = c.size

547 if n == 0:

548 return np.zeros(0)

549 if n == 1:

550 return np.array([c[0]])

551

552 # Build Chebyshev coefficients via the Legendre recurrence.

553 # The Chebyshev representation of P_j is computed column by column.

554 # Multiplication by x in Chebyshev basis:

555 # (x*f)[0] = f[1]/2

556 # (x*f)[1] = f[0] + f[2]/2

557 # (x*f)[k] = (f[k-1] + f[k+1])/2 for k >= 2

558 result = np.zeros(n)

559

560 prev_prev = np.zeros(n)

561 prev_prev[0] = 1.0 # P_0 = T_0

562 result += c[0] * prev_prev

563

564 prev = np.zeros(n)

565 prev[1] = 1.0 # P_1 = T_1

566 result += c[1] * prev

567

568 for j in range(2, n):

569 # x * prev in Chebyshev basis

570 xprev = np.zeros(n)

571 xprev[1] += prev[0] # from x*T_0 = T_1

572 xprev[: n - 1] += prev[1:] / 2.0 # T_{k-1} from x*T_k for k>=1

573 xprev[2:] += prev[1 : n - 1] / 2.0 # T_{k+1} from x*T_k for k>=1

574

575 curr = ((2 * j - 1) * xprev - (j - 1) * prev_prev) / j

576 result += c[j] * curr

577 prev_prev = prev

578 prev = curr

579

580 return result

581

582

583def _conv_legendre(a: np.ndarray, b: np.ndarray) -> tuple[np.ndarray, np.ndarray]:

584 """Convolve two Legendre series using the Hale-Townsend algorithm.

585

586 Computes the convolution of two functions expressed as Legendre series on

587 [-1, 1]. The result is a piecewise polynomial on [-2, 2], split into a

588 left piece on [-2, 0] and a right piece on [0, 2]. The Legendre

589 coefficients of each piece (with respect to the linear map of the piece

590 to [-1, 1]) are returned.

591

592 The algorithm is based on:

593 N. Hale and A. Townsend, "An algorithm for the convolution of Legendre

594 series", SIAM J. Sci. Comput., 36(3), A1207-A1220, 2014.

595

596 Args:

597 a (array-like): Legendre coefficients of the first function on [-1, 1].

598 b (array-like): Legendre coefficients of the second function on [-1, 1].

599

600 Returns:

601 tuple: (gamma_left, gamma_right) where each element is a 1-D array of

602 Legendre coefficients for the left [-2, 0] and right [0, 2] pieces

603 respectively.

604 """

605 a = np.asarray(a, dtype=float).ravel()

606 b = np.asarray(b, dtype=float).ravel()

607

608 # Ensure a has the higher (or equal) degree

609 if len(b) > len(a):

610 a, b = b, a

611

612 na, nb = len(a), len(b)

613 mn = na + nb

614

615 # Pad a to length mn

616 alpha = np.zeros(mn)

617 alpha[:na] = a

618

619 # Build the tridiagonal S matrix (mn x mn), the Legendre cumulative-integral

620 # operator (S f)(x) = ∫_{-1}^{x} f(t) dt. Entries (using column index n):

621 # S[0, 0] = 1

622 # S[n+1, n] = 1/(2n+1) for n >= 0 (sub-diagonal)

623 # S[n-1, n] = -1/(2n+1) for n >= 1 (super-diagonal)

624 k = np.arange(mn)

625 main = np.zeros(mn)

626 main[0] = 1.0

627 sub = 1.0 / (2.0 * k[:-1] + 1.0) # [1, 1/3, 1/5, ...], length mn-1

628 supra = -1.0 / (2.0 * k[1:] + 1.0) # [-1/3, -1/5, -1/7, ...], length mn-1

629

630 def _s_apply(v: np.ndarray) -> np.ndarray:

631 """Apply the S matrix to vector v."""

632 res = main * v

633 res[1:] += sub * v[:-1]

634 res[:-1] += supra * v[1:]

635 return res

636

637 def _rec(alpha_arg: np.ndarray, beta: np.ndarray, sgn: float, s00: float) -> np.ndarray:

638 """Compute Legendre coefficients of the convolution on one piece.

639

640 Uses the recurrence from Theorem 4.1 of Hale & Townsend (2014).

641 """

642 n_beta = len(beta)

643 # Save / restore main[0] for S

644 save_main0 = main[0]

645 main[0] = s00

646

647 # scl[k] = (-1)^k / (2k-1) for k=1,...,n_beta (1-indexed)

648 scl = np.ones(n_beta) / (2.0 * np.arange(1, n_beta + 1) - 1.0)

649 scl[1::2] = -scl[1::2]

650

651 # First column

652 v_new = _s_apply(alpha_arg)

653 v = v_new.copy()

654 gamma = beta[0] * v_new.copy()

655 beta_scl = scl * beta

656 beta_scl[0] = 0.0

657 gamma[0] += float(v_new[:n_beta].dot(beta_scl))

658

659 if n_beta > 1:

660 # Second column

661 v_new = _s_apply(v) + sgn * v

662 v_old = v.copy()

663 v = v_new.copy()

664 v_new[0] = 0.0

665 gamma += beta[1] * v_new

666 beta_scl = -beta_scl * (2.0 - 0.5) / (2.0 - 1.5)

667 beta_scl[1] = 0.0

668 gamma[1] += float(v_new[:n_beta].dot(beta_scl))

669

670 # Remaining columns

671 for nn in range(3, n_beta + 1):

672 v_new = (2 * nn - 3) * _s_apply(v) + v_old

673 v_new[: nn - 1] = 0.0

674 gamma += v_new * beta[nn - 1]

675 beta_scl = -beta_scl * (nn - 0.5) / (nn - 1.5)

676 beta_scl[nn - 1] = 0.0

677 gamma[nn - 1] += float(v_new[:n_beta].dot(beta_scl))

678 v_old = v.copy()

679 v = v_new.copy()

680

681 # Restore

682 main[0] = save_main0

683

684 # Trim trailing near-zeros

685 ag = np.abs(gamma)

686 mg = np.max(ag) if ag.size > 0 else 0.0

687 if mg > 0:

688 loc = np.where(ag > np.finfo(float).eps * mg)[0]

689 gamma = gamma[: loc[-1] + 1] if loc.size > 0 else gamma[:1]

690 else:

691 gamma = gamma[:1]

692 return gamma

693

694 gamma_left = _rec(alpha.copy(), b, -1.0, 1.0)

695 gamma_right = _rec(-alpha.copy(), b, 1.0, -1.0)

696

697 return gamma_left, gamma_right

698

699

700def newtonroots(fun: Any, rts: np.ndarray, tol: float | None = None, maxiter: int | None = None) -> np.ndarray:

701 """Refine root approximations using Newton's method.

702

703 This function applies Newton's method to refine the approximations of roots

704 for a callable and differentiable function. It is typically used to polish

705 already computed roots to higher accuracy.

706

707 Args:

708 fun (callable): A callable and differentiable function.

709 rts (numpy.ndarray): Initial approximations of the roots.

710 tol (float, optional): Tolerance for convergence. Defaults to 2 * machine epsilon.

711 maxiter (int, optional): Maximum number of iterations. Defaults to value from preferences.

712

713 Returns:

714 numpy.ndarray: Refined approximations of the roots.

715

716 Note:

717 The function must support differentiation via a .diff() method that returns

718 the derivative function.

719 """

720 tol = tol if tol is not None else 2 * prefs.eps

721 maxiter = maxiter if maxiter is not None else prefs.maxiter

722 if rts.size > 0:

723 dfun = fun.diff()

724 prv = np.inf * rts

725 count = 0

726 while (infnorm(rts - prv) > tol) & (count <= maxiter):

727 count += 1

728 prv = rts

729 rts = rts - fun(rts) / dfun(rts)

730 return rts

Coverage for src / chebpy / algorithms.py: 100%

304 statements