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

325 statements  

« prev     ^ index     » next       coverage.py v7.16.2, created at 2026-10-01 13:43 +0000

1"""Implementation of the Chebfun class for piecewise function approximation. 

2 

3This module provides the Chebfun class, which is the main user-facing class in the 

4ChebPy package. It represents functions using piecewise polynomial approximations 

5on arbitrary intervals, allowing for operations such as integration, differentiation, 

6root-finding, and more. 

7 

8The Chebfun class is inspired by the MATLAB package of the same name and provides 

9similar functionality for working with functions rather than numbers. 

10""" 

11 

12from __future__ import annotations 

13 

14import operator 

15from collections.abc import Callable, Iterator 

16from typing import Any, cast 

17 

18import numpy as np 

19from matplotlib.axes import Axes 

20 

21from . import _pointwise 

22from ._construction import generate_funs 

23from ._convolution import convolve 

24from ._singular_construction import generate_singular_funs 

25from ._ufuncs import register_ufuncs 

26from .bndfun import Bndfun 

27from .decorators import cache, cast_arg_to_chebfun, float_argument, self_empty 

28from .exceptions import BadFunLengthArgument 

29from .plotting import plot_chebfun, plotcoeffs_chebfun 

30from .settings import _preferences as prefs 

31from .utilities import Domain, check_funs, compute_breakdata 

32 

33 

34class Chebfun: 

35 """Main class for representing and manipulating functions in ChebPy. 

36 

37 The Chebfun class represents functions using piecewise polynomial approximations 

38 on arbitrary intervals. It provides a comprehensive set of operations for working 

39 with these function representations, including: 

40 

41 - Function evaluation at arbitrary points 

42 - Algebraic operations (addition, multiplication, etc.) 

43 - Calculus operations (differentiation, integration, etc.) 

44 - Rootfinding 

45 - Plotting 

46 

47 Chebfun objects can be created from callable functions, constant values, or 

48 directly from function pieces. The class supports both adaptive and fixed-length 

49 approximations, allowing for efficient representation of functions with varying 

50 complexity across different intervals. 

51 

52 Attributes: 

53 funs (numpy.ndarray): Array of function pieces that make up the Chebfun. 

54 breakdata (OrderedDict): Mapping of breakpoints to function values. 

55 transposed (bool): Flag indicating if the Chebfun is transposed. 

56 """ 

57 

58 def __init__(self, funs: Any) -> None: 

59 """Initialize a Chebfun object. 

60 

61 Args: 

62 funs (list): List of function objects to be included in the Chebfun. 

63 These will be checked and sorted using check_funs. 

64 """ 

65 self.funs = check_funs(funs) 

66 self.breakdata = compute_breakdata(self.funs) 

67 self.transposed = False 

68 

69 @classmethod 

70 def initempty(cls) -> Chebfun: 

71 """Initialize an empty Chebfun. 

72 

73 Returns: 

74 Chebfun: An empty Chebfun object with no functions. 

75 

76 Examples: 

77 >>> f = Chebfun.initempty() 

78 >>> f.isempty 

79 True 

80 """ 

81 return cls([]) 

82 

83 @classmethod 

84 def initidentity(cls, domain: Any = None) -> Chebfun: 

85 """Initialize a Chebfun representing the identity function f(x) = x. 

86 

87 Args: 

88 domain (array-like, optional): Domain on which to define the identity function. 

89 If None, uses the default domain from preferences. 

90 

91 Returns: 

92 Chebfun: A Chebfun object representing the identity function on the specified domain. 

93 

94 Examples: 

95 >>> import numpy as np 

96 >>> x = Chebfun.initidentity([-1, 1]) 

97 >>> float(x(0.5)) 

98 0.5 

99 >>> np.allclose(x([0, 0.5, 1]), [0, 0.5, 1]) 

100 True 

101 """ 

102 return cls(generate_funs(domain, Bndfun.initidentity)) 

103 

104 @classmethod 

105 def initconst(cls, c: Any, domain: Any = None) -> Chebfun: 

106 """Initialize a Chebfun representing a constant function f(x) = c. 

107 

108 Args: 

109 c (float or complex): The constant value. 

110 domain (array-like, optional): Domain on which to define the constant function. 

111 If None, uses the default domain from preferences. 

112 

113 Returns: 

114 Chebfun: A Chebfun object representing the constant function on the specified domain. 

115 

116 Examples: 

117 >>> import numpy as np 

118 >>> f = Chebfun.initconst(3.14, [-1, 1]) 

119 >>> float(f(0)) 

120 3.14 

121 >>> float(f(0.5)) 

122 3.14 

123 >>> np.allclose(f([0, 0.5, 1]), [3.14, 3.14, 3.14]) 

124 True 

125 """ 

126 return cls(generate_funs(domain, Bndfun.initconst, {"c": c})) 

127 

128 @classmethod 

129 def initfun_adaptive( 

130 cls, 

131 f: Callable[..., Any], 

132 domain: Any = None, 

133 *, 

134 sing: str | None = None, 

135 params: Any = None, 

136 ) -> Chebfun: 

137 """Initialize a Chebfun by adaptively sampling a function. 

138 

139 This method determines the appropriate number of points needed to represent 

140 the function to the specified tolerance using an adaptive algorithm. 

141 

142 Args: 

143 f (callable): The function to be approximated. 

144 domain (array-like, optional): Domain on which to define the function. 

145 If None, uses the default domain from preferences. 

146 sing: Optional endpoint-singularity hint, one of ``"left"``, 

147 ``"right"``, or ``"both"``. When set, the appropriate boundary 

148 pieces are built as :class:`~chebpy.singfun.Singfun` instances 

149 using the Adcock-Richardson clustering map; interior pieces 

150 remain :class:`~chebpy.bndfun.Bndfun`. 

151 params: Slit-strip map parameters (a :class:`~chebpy.maps.MapParams`). 

152 Ignored when ``sing`` is ``None``. Default ``None`` (uses 

153 :class:`~chebpy.maps.MapParams` defaults). 

154 

155 Returns: 

156 Chebfun: A Chebfun object representing the function on the specified domain. 

157 

158 Examples: 

159 >>> import numpy as np 

160 >>> f = Chebfun.initfun_adaptive(lambda x: np.sin(x), [-np.pi, np.pi]) 

161 >>> bool(abs(f(0)) < 1e-10) 

162 True 

163 >>> bool(abs(f(np.pi/2) - 1) < 1e-10) 

164 True 

165 """ 

166 if sing is None: 

167 return cls(generate_funs(domain, Bndfun.initfun_adaptive, {"f": f})) 

168 

169 return cls(generate_singular_funs(f, domain, sing=sing, params=params)) 

170 

171 @classmethod 

172 def initfun_fixedlen(cls, f: Callable[..., Any], n: Any, domain: Any = None) -> Chebfun: 

173 """Initialize a Chebfun with a fixed number of points. 

174 

175 This method uses a specified number of points to represent the function, 

176 rather than determining the number adaptively. 

177 

178 Args: 

179 f (callable): The function to be approximated. 

180 n (int or array-like): Number of points to use. If a single value, uses the same 

181 number for each interval. If an array, must have one fewer elements than 

182 the size of the domain. 

183 domain (array-like, optional): Domain on which to define the function. 

184 If None, uses the default domain from preferences. 

185 

186 Returns: 

187 Chebfun: A Chebfun object representing the function on the specified domain. 

188 

189 Raises: 

190 BadFunLengthArgument: If n is an array and its size doesn't match domain.size - 1. 

191 """ 

192 nn = np.array(n) 

193 if nn.size < 2: 

194 funs = generate_funs(domain, Bndfun.initfun_fixedlen, {"f": f, "n": n}) 

195 else: 

196 domain = Domain(domain if domain is not None else prefs.domain) 

197 if not nn.size == domain.size - 1: 

198 raise BadFunLengthArgument 

199 funs = [] 

200 for interval, length in zip(domain.intervals, nn, strict=False): 

201 funs.append(Bndfun.initfun_fixedlen(f, interval, length)) 

202 return cls(funs) 

203 

204 @classmethod 

205 def initfun( 

206 cls, 

207 f: Callable[..., Any], 

208 domain: Any = None, 

209 n: Any = None, 

210 *, 

211 sing: str | None = None, 

212 params: Any = None, 

213 ) -> Chebfun: 

214 """Initialize a Chebfun from a function. 

215 

216 This is a general-purpose constructor that delegates to either initfun_adaptive 

217 or initfun_fixedlen based on whether n is provided. 

218 

219 Args: 

220 f (callable): The function to be approximated. 

221 domain (array-like, optional): Domain on which to define the function. 

222 If None, uses the default domain from preferences. 

223 n (int or array-like, optional): Number of points to use. If None, determines 

224 the number adaptively. If provided, uses a fixed number of points. 

225 sing: Optional endpoint-singularity hint forwarded to 

226 :meth:`initfun_adaptive`. Only valid when ``n is None``. 

227 params: Slit-strip map parameters (a :class:`~chebpy.maps.MapParams`). 

228 Forwarded to :meth:`initfun_adaptive`. 

229 

230 Returns: 

231 Chebfun: A Chebfun object representing the function on the specified domain. 

232 """ 

233 if n is None: 

234 return cls.initfun_adaptive(f, domain, sing=sing, params=params) 

235 if sing is not None: 

236 msg = ( 

237 "fixed-length construction with sing= is not supported in v1; " 

238 "pass n=None for adaptive Singfun construction." 

239 ) 

240 raise NotImplementedError(msg) 

241 return cls.initfun_fixedlen(f, n, domain) 

242 

243 # -------------------- 

244 # operator overloads 

245 # -------------------- 

246 def __add__(self, f: Any) -> Any: 

247 """Add a Chebfun with another Chebfun or a scalar. 

248 

249 Args: 

250 f (Chebfun or scalar): The object to add to this Chebfun. 

251 

252 Returns: 

253 Chebfun: A new Chebfun representing the sum. 

254 """ 

255 return self._apply_binop(f, operator.add) 

256 

257 @self_empty(np.array([])) 

258 @float_argument 

259 def __call__(self, x: Any) -> Any: 

260 """Evaluate the Chebfun at points x. 

261 

262 This method evaluates the Chebfun at the specified points. It handles interior 

263 points, breakpoints, and points outside the domain appropriately. 

264 

265 Args: 

266 x (float or array-like): Points at which to evaluate the Chebfun. 

267 

268 Returns: 

269 float or numpy.ndarray: The value(s) of the Chebfun at the specified point(s). 

270 Returns a scalar if x is a scalar, otherwise an array of the same size as x. 

271 """ 

272 # initialise output 

273 dtype = complex if self.iscomplex else float 

274 out = np.full(x.size, np.nan, dtype=dtype) 

275 

276 # evaluate a fun when x is an interior point 

277 for fun in self: 

278 sa, sb = fun.support[0], fun.support[-1] 

279 idx = np.logical_and(sa < x, x < sb) 

280 out[idx] = fun(x[idx]) 

281 

282 # evaluate the breakpoint data for x at a breakpoint 

283 breakpoints = self.breakpoints 

284 for break_point in breakpoints: 

285 out[x == break_point] = self.breakdata[break_point] 

286 

287 # first and last funs used to evaluate outside of the chebfun domain 

288 lpts, rpts = x < breakpoints[0], x > breakpoints[-1] 

289 out[lpts] = self.funs[0](x[lpts]) 

290 out[rpts] = self.funs[-1](x[rpts]) 

291 return out 

292 

293 def __iter__(self) -> Iterator[Any]: 

294 """Return an iterator over the functions in this Chebfun. 

295 

296 Returns: 

297 iterator: An iterator over the functions (funs) in this Chebfun. 

298 """ 

299 return self.funs.__iter__() 

300 

301 def __len__(self) -> int: 

302 """Return the total number of coefficients across all funs. 

303 

304 Returns: 

305 int: The sum of sizes of all constituent funs. 

306 """ 

307 return sum(f.size for f in self.funs) 

308 

309 def __eq__(self, other: object) -> bool: 

310 """Test for equality between two Chebfun objects. 

311 

312 Two Chebfun objects are considered equal if they have the same domain 

313 and their function values are equal (within tolerance) at a set of test points. 

314 

315 Args: 

316 other (object): The object to compare with this Chebfun. 

317 

318 Returns: 

319 bool: True if the objects are equal, False otherwise. 

320 """ 

321 if not isinstance(other, self.__class__): 

322 return False 

323 

324 # Check if both are empty 

325 if self.isempty and other.isempty: 

326 return True 

327 

328 # Check if domains are equal 

329 if self.domain != other.domain: 

330 return False 

331 

332 # Check function values at test points 

333 xx = np.linspace(self.support[0], self.support[1], 100) 

334 tol = 1e2 * max(self.vscale, other.vscale) * prefs.eps 

335 return bool(np.all(np.abs(self(xx) - other(xx)) <= tol)) 

336 

337 def __mul__(self, f: Any) -> Any: 

338 """Multiply a Chebfun with another Chebfun or a scalar. 

339 

340 Args: 

341 f (Chebfun or scalar): The object to multiply with this Chebfun. 

342 

343 Returns: 

344 Chebfun: A new Chebfun representing the product. 

345 """ 

346 return self._apply_binop(f, operator.mul) 

347 

348 def __neg__(self) -> Chebfun: 

349 """Return the negative of this Chebfun. 

350 

351 Returns: 

352 Chebfun: A new Chebfun representing -f(x). 

353 """ 

354 return self.__class__(-self.funs) 

355 

356 def __pos__(self) -> Chebfun: 

357 """Return the positive of this Chebfun (which is the Chebfun itself). 

358 

359 Returns: 

360 Chebfun: This Chebfun object (unchanged). 

361 """ 

362 return self 

363 

364 def __abs__(self) -> Chebfun: 

365 """Return the absolute value of this Chebfun. 

366 

367 Returns: 

368 Chebfun: A new Chebfun representing |f(x)|. 

369 """ 

370 abs_funs = [] 

371 for fun in self.funs: 

372 abs_funs.append(fun.absolute()) 

373 return self.__class__(abs_funs) 

374 

375 def __pow__(self, f: Any) -> Any: 

376 """Raise this Chebfun to a power. 

377 

378 Args: 

379 f (Chebfun or scalar): The exponent to which this Chebfun is raised. 

380 

381 Returns: 

382 Chebfun: A new Chebfun representing self^f. 

383 """ 

384 return self._apply_binop(f, operator.pow) 

385 

386 def __rtruediv__(self, c: Any) -> Chebfun: 

387 """Divide a scalar by this Chebfun. 

388 

389 This method is called when a scalar is divided by a Chebfun, i.e., c / self. 

390 

391 Args: 

392 c (scalar): The scalar numerator. 

393 

394 Returns: 

395 Chebfun: A new Chebfun representing c / self. 

396 

397 Note: 

398 This is executed when truediv(f, self) fails, which is to say whenever c 

399 is not a Chebfun. We proceed on the assumption f is a scalar. 

400 """ 

401 

402 def constfun(cheb: Any, const: Any) -> Any: 

403 return 0.0 * cheb + const 

404 

405 def make_divfun(fun: Any) -> Callable[..., Any]: 

406 return lambda x: constfun(x, c) / fun(x) 

407 

408 newfuns = [fun.initfun_adaptive(make_divfun(fun), fun.interval) for fun in self] 

409 return self.__class__(newfuns) 

410 

411 @self_empty("Chebfun<empty>") 

412 def __repr__(self) -> str: 

413 """Return a string representation of the Chebfun. 

414 

415 This method returns a detailed string representation of the Chebfun, 

416 including information about its domain, intervals, and endpoint values. 

417 

418 Returns: 

419 str: A string representation of the Chebfun. 

420 """ 

421 rowcol = "row" if self.transposed else "column" 

422 numpcs = self.funs.size 

423 plural = "" if numpcs == 1 else "s" 

424 header = f"Chebfun {rowcol} ({numpcs} smooth piece{plural})\n" 

425 toprow = " interval length endpoint values\n" 

426 tmplat = "[{:8.2g},{:8.2g}] {:6} {:8.2g} {:8.2g}\n" 

427 rowdta = "" 

428 for fun in self: 

429 endpts = fun.support 

430 xl, xr = endpts 

431 fl, fr = fun(endpts) 

432 row = tmplat.format(xl, xr, fun.size, fl, fr) 

433 rowdta += row 

434 btmrow = f"vertical scale = {self.vscale:3.2g}" 

435 btmxtr = "" if numpcs == 1 else f" total length = {sum([f.size for f in self])}" 

436 return header + toprow + rowdta + btmrow + btmxtr 

437 

438 def __rsub__(self, f: Any) -> Any: 

439 """Subtract this Chebfun from another object. 

440 

441 This method is called when another object is subtracted by this Chebfun, 

442 i.e., f - self. 

443 

444 Args: 

445 f (Chebfun or scalar): The object from which to subtract this Chebfun. 

446 

447 Returns: 

448 Chebfun: A new Chebfun representing f - self. 

449 """ 

450 return -(self - f) 

451 

452 @cast_arg_to_chebfun 

453 def __rpow__(self, f: Any) -> Any: 

454 """Raise another object to the power of this Chebfun. 

455 

456 This method is called when another object is raised to the power of this Chebfun, 

457 i.e., f ** self. 

458 

459 Args: 

460 f (Chebfun or scalar): The base to be raised to the power of this Chebfun. 

461 

462 Returns: 

463 Chebfun: A new Chebfun representing f ** self. 

464 """ 

465 return f**self 

466 

467 def __truediv__(self, f: Any) -> Any: 

468 """Divide this Chebfun by another object. 

469 

470 Args: 

471 f (Chebfun or scalar): The divisor. 

472 

473 Returns: 

474 Chebfun: A new Chebfun representing self / f. 

475 """ 

476 return self._apply_binop(f, operator.truediv) 

477 

478 __rmul__ = __mul__ 

479 __div__ = __truediv__ 

480 __rdiv__ = __rtruediv__ 

481 __radd__ = __add__ 

482 

483 def __str__(self) -> str: 

484 """Return a human-readable string representation of the Chebfun. 

485 

486 This method returns the same detailed representation as ``__repr__``, 

487 so that ``print(f)`` shows the full summary table. This is consistent 

488 with the behaviour of numpy and pandas objects. 

489 

490 Returns: 

491 str: A detailed string representation of the Chebfun. 

492 """ 

493 return repr(self) 

494 

495 def __sub__(self, f: Any) -> Any: 

496 """Subtract another object from this Chebfun. 

497 

498 Args: 

499 f (Chebfun or scalar): The object to subtract from this Chebfun. 

500 

501 Returns: 

502 Chebfun: A new Chebfun representing self - f. 

503 """ 

504 return self._apply_binop(f, operator.sub) 

505 

506 # ------------------ 

507 # internal helpers 

508 # ------------------ 

509 @self_empty() 

510 def _apply_binop(self, f: Any, op: Callable[..., Any]) -> Any: 

511 """Apply a binary operation between this Chebfun and another object. 

512 

513 This is a funnel method used in the implementation of Chebfun binary 

514 operators. The high-level idea is to first break each chebfun into a 

515 series of pieces corresponding to the union of the domains of each 

516 before applying the supplied binary operator and simplifying. In the 

517 case of the second argument being a scalar we don't need to do the 

518 simplify step, since at the Tech-level these operations are defined 

519 such that there is no change in the number of coefficients. 

520 

521 Args: 

522 f (Chebfun or scalar): The second operand of the binary operation. 

523 op (callable): The binary operation to apply (e.g., operator.add). 

524 

525 Returns: 

526 Chebfun: A new Chebfun resulting from applying the binary operation. 

527 """ 

528 if hasattr(f, "isempty") and f.isempty: 

529 return f 

530 if np.isscalar(f): 

531 chbfn1 = self 

532 chbfn2 = cast(Any, f) * np.ones(self.funs.size) 

533 simplify = False 

534 else: 

535 newdom = self.domain.union(f.domain) 

536 chbfn1 = self._break(newdom) 

537 chbfn2 = f._break(newdom) 

538 simplify = True 

539 newfuns = [] 

540 for fun1, fun2 in zip(chbfn1, chbfn2, strict=False): 

541 newfun = op(fun1, fun2) 

542 if simplify: 

543 newfun = newfun.simplify() 

544 newfuns.append(newfun) 

545 return self.__class__(newfuns) 

546 

547 def _break(self, targetdomain: Domain) -> Chebfun: 

548 """Resample this Chebfun to a new domain. 

549 

550 This method resamples the Chebfun to the supplied Domain object. It is 

551 intended as a private method since one will typically need to have 

552 called either Domain.union(f) or Domain.merge(f) prior to calling this method. 

553 

554 Args: 

555 targetdomain (Domain): The domain to which this Chebfun should be resampled. 

556 

557 Returns: 

558 Chebfun: A new Chebfun resampled to the target domain. 

559 """ 

560 newfuns = [] 

561 subintervals = iter(targetdomain.intervals) 

562 interval = next(subintervals) # next(..) for Python2/3 compatibility 

563 for fun in self: 

564 while interval in fun.interval: 

565 newfun = fun.restrict(interval) 

566 newfuns.append(newfun) 

567 try: 

568 interval = next(subintervals) 

569 except StopIteration: 

570 break 

571 return self.__class__(newfuns) 

572 

573 # ------------ 

574 # properties 

575 # ------------ 

576 @property 

577 def breakpoints(self) -> np.ndarray: 

578 """Get the breakpoints of this Chebfun. 

579 

580 Breakpoints are the points where the Chebfun transitions from one piece to another. 

581 

582 Returns: 

583 numpy.ndarray: Array of breakpoints. 

584 """ 

585 return np.array(list(self.breakdata.keys())) 

586 

587 @property 

588 @self_empty(Domain([])) 

589 def domain(self) -> Domain: 

590 """Get the domain of this Chebfun. 

591 

592 Returns: 

593 Domain: A Domain object corresponding to this Chebfun. 

594 """ 

595 return Domain.from_chebfun(self) 

596 

597 @domain.setter 

598 def domain(self, new_domain: Any) -> None: 

599 """Set the domain of the Chebfun by restricting to the new domain. 

600 

601 Args: 

602 new_domain (array-like): The new domain to which this Chebfun should be restricted. 

603 """ 

604 self.restrict_(new_domain) 

605 

606 @property 

607 @self_empty(Domain([])) 

608 def support(self) -> Any: 

609 """Get the support interval of this Chebfun. 

610 

611 The support is the interval between the first and last breakpoints. 

612 

613 Returns: 

614 numpy.ndarray: Array containing the first and last breakpoints. 

615 """ 

616 return self.domain.support 

617 

618 @property 

619 @self_empty(0.0) 

620 def hscale(self) -> float: 

621 """Get the horizontal scale of this Chebfun. 

622 

623 The horizontal scale is the maximum absolute value of the support interval. 

624 

625 Returns: 

626 float: The horizontal scale. 

627 """ 

628 return float(np.abs(self.support).max()) 

629 

630 @property 

631 @self_empty(False) 

632 def iscomplex(self) -> bool: 

633 """Check if this Chebfun has complex values. 

634 

635 Returns: 

636 bool: True if any of the functions in this Chebfun have complex values, 

637 False otherwise. 

638 """ 

639 return any(fun.iscomplex for fun in self) 

640 

641 @property 

642 @self_empty(False) 

643 def isconst(self) -> bool: 

644 """Check if this Chebfun represents a constant function. 

645 

646 A Chebfun is constant if all of its pieces are constant with the same value. 

647 

648 Returns: 

649 bool: True if this Chebfun represents a constant function, False otherwise. 

650 

651 Note: 

652 TODO: find an abstract way of referencing funs[0].coeffs[0] 

653 """ 

654 c = self.funs[0].coeffs[0] 

655 return all(fun.isconst and fun.coeffs[0] == c for fun in self) 

656 

657 @property 

658 def isempty(self) -> bool: 

659 """Check if this Chebfun is empty. 

660 

661 An empty Chebfun contains no functions. 

662 

663 Returns: 

664 bool: True if this Chebfun is empty, False otherwise. 

665 """ 

666 return self.funs.size == 0 

667 

668 @property 

669 @self_empty(0.0) 

670 def vscale(self) -> Any: 

671 """Get the vertical scale of this Chebfun. 

672 

673 The vertical scale is the maximum of the vertical scales of all pieces. 

674 

675 Returns: 

676 float: The vertical scale. 

677 """ 

678 return np.max([fun.vscale for fun in self]) 

679 

680 @property 

681 @self_empty() 

682 def x(self) -> Chebfun: 

683 """Get the identity function on the support of this Chebfun. 

684 

685 This property returns a new Chebfun representing the identity function f(x) = x 

686 defined on the same support as this Chebfun. 

687 

688 Returns: 

689 Chebfun: A Chebfun representing the identity function on the support of this Chebfun. 

690 """ 

691 return self.__class__.initidentity(self.support) 

692 

693 # ----------- 

694 # utilities 

695 # ---------- 

696 

697 def imag(self) -> Chebfun: 

698 """Get the imaginary part of this Chebfun. 

699 

700 Returns: 

701 Chebfun: A new Chebfun representing the imaginary part of this Chebfun. 

702 If this Chebfun is real-valued, returns a zero Chebfun. 

703 """ 

704 if self.iscomplex: 

705 return self.__class__([fun.imag() for fun in self]) 

706 else: 

707 return self.initconst(0, domain=self.domain) 

708 

709 def real(self) -> Chebfun: 

710 """Get the real part of this Chebfun. 

711 

712 Returns: 

713 Chebfun: A new Chebfun representing the real part of this Chebfun. 

714 If this Chebfun is already real-valued, returns this Chebfun. 

715 """ 

716 if self.iscomplex: 

717 return self.__class__([fun.real() for fun in self]) 

718 else: 

719 return self 

720 

721 def copy(self) -> Chebfun: 

722 """Create a deep copy of this Chebfun. 

723 

724 Returns: 

725 Chebfun: A new Chebfun that is a deep copy of this Chebfun. 

726 """ 

727 return self.__class__([fun.copy() for fun in self]) 

728 

729 @self_empty() 

730 def _restrict(self, subinterval: Any) -> Chebfun: 

731 """Restrict a Chebfun to a subinterval, without simplifying. 

732 

733 This is an internal method that restricts the Chebfun to a subinterval 

734 without performing simplification. 

735 

736 Args: 

737 subinterval (array-like): The subinterval to which this Chebfun should be restricted. 

738 

739 Returns: 

740 Chebfun: A new Chebfun restricted to the specified subinterval, without simplification. 

741 """ 

742 newdom = self.domain.restrict(Domain(subinterval)) 

743 return self._break(newdom) 

744 

745 def restrict(self, subinterval: Any) -> Any: 

746 """Restrict a Chebfun to a subinterval. 

747 

748 This method creates a new Chebfun that is restricted to the specified subinterval 

749 and simplifies the result. 

750 

751 Args: 

752 subinterval (array-like): The subinterval to which this Chebfun should be restricted. 

753 

754 Returns: 

755 Chebfun: A new Chebfun restricted to the specified subinterval. 

756 """ 

757 return self._restrict(subinterval).simplify() 

758 

759 @self_empty() 

760 def restrict_(self, subinterval: Any) -> Chebfun: 

761 """Restrict a Chebfun to a subinterval, modifying the object in place. 

762 

763 This method modifies the current Chebfun by restricting it to the specified 

764 subinterval and simplifying the result. 

765 

766 Args: 

767 subinterval (array-like): The subinterval to which this Chebfun should be restricted. 

768 

769 Returns: 

770 Chebfun: The modified Chebfun (self). 

771 """ 

772 restricted = self._restrict(subinterval).simplify() 

773 self.funs = restricted.funs 

774 self.breakdata = compute_breakdata(self.funs) 

775 return self 

776 

777 @cache 

778 @self_empty(np.array([])) 

779 def roots(self, merge: Any = None) -> np.ndarray: 

780 """Compute the roots of a Chebfun. 

781 

782 This method finds the values x for which f(x) = 0, by computing the roots 

783 of each piece of the Chebfun and combining them. 

784 

785 Args: 

786 merge (bool, optional): Whether to merge roots at breakpoints. If None, 

787 uses the value from preferences. Defaults to None. 

788 

789 Returns: 

790 numpy.ndarray: Array of roots sorted in ascending order. 

791 

792 Examples: 

793 >>> import numpy as np 

794 >>> f = Chebfun.initfun_adaptive(lambda x: x**2 - 1, [-2, 2]) 

795 >>> roots = f.roots() 

796 >>> len(roots) 

797 2 

798 >>> np.allclose(sorted(roots), [-1, 1]) 

799 True 

800 """ 

801 merge = merge if merge is not None else prefs.mergeroots 

802 allrts = [] 

803 prvrts = np.array([]) 

804 htol = 1e2 * self.hscale * prefs.eps 

805 for fun in self: 

806 rts = fun.roots() 

807 # ignore first root if equal to the last root of previous fun 

808 # TODO: there could be multiple roots at breakpoints 

809 if prvrts.size > 0 and rts.size > 0 and merge and abs(prvrts[-1] - rts[0]) <= htol: 

810 rts = rts[1:] 

811 allrts.append(rts) 

812 prvrts = rts 

813 return np.concatenate(list(allrts)) 

814 

815 @self_empty() 

816 def simplify(self) -> Chebfun: 

817 """Simplify each fun in the chebfun.""" 

818 return self.__class__([fun.simplify() for fun in self]) 

819 

820 def translate(self, c: Any) -> Chebfun: 

821 """Translate a chebfun by c, i.e., return f(x-c).""" 

822 return self.__class__([x.translate(c) for x in self]) 

823 

824 # ---------- 

825 # calculus 

826 # ---------- 

827 def cumsum(self) -> Chebfun: 

828 """Compute the indefinite integral (antiderivative) of the Chebfun. 

829 

830 This method computes the indefinite integral of the Chebfun, with the 

831 constant of integration chosen so that the indefinite integral evaluates 

832 to 0 at the left endpoint of the domain. For piecewise functions, constants 

833 are added to ensure continuity across the pieces. 

834 

835 Returns: 

836 Chebfun: A new Chebfun representing the indefinite integral of this Chebfun. 

837 

838 Examples: 

839 >>> import numpy as np 

840 >>> f = Chebfun.initconst(1.0, [-1, 1]) 

841 >>> F = f.cumsum() 

842 >>> bool(abs(F(-1)) < 1e-10) 

843 True 

844 >>> bool(abs(F(1) - 2.0) < 1e-10) 

845 True 

846 """ 

847 newfuns = [] 

848 prevfun = None 

849 for fun in self: 

850 integral = fun.cumsum() 

851 if prevfun: 

852 # enforce continuity by adding the function value 

853 # at the right endpoint of the previous fun 

854 _, fb = prevfun.endvalues 

855 integral = integral + fb 

856 newfuns.append(integral) 

857 prevfun = integral 

858 return self.__class__(newfuns) 

859 

860 def diff(self, n: int = 1) -> Chebfun: 

861 """Compute the derivative of the Chebfun. 

862 

863 This method calculates the nth derivative of the Chebfun with respect to x. 

864 It creates a new Chebfun where each piece is the derivative of the 

865 corresponding piece in the original Chebfun. 

866 

867 Args: 

868 n: Order of differentiation (default: 1). Must be non-negative integer. 

869 

870 Returns: 

871 Chebfun: A new Chebfun representing the nth derivative of this Chebfun. 

872 

873 Examples: 

874 >>> from chebpy import chebfun 

875 >>> f = chebfun(lambda x: x**3) 

876 >>> df1 = f.diff() # first derivative: 3*x**2 

877 >>> df2 = f.diff(2) # second derivative: 6*x 

878 >>> df3 = f.diff(3) # third derivative: 6 

879 >>> bool(abs(df1(0.5) - 0.75) < 1e-10) 

880 True 

881 >>> bool(abs(df2(0.5) - 3.0) < 1e-10) 

882 True 

883 >>> bool(abs(df3(0.5) - 6.0) < 1e-10) 

884 True 

885 """ 

886 if not isinstance(n, int): 

887 raise TypeError(n) 

888 if n == 0: 

889 return self 

890 if n < 0: 

891 raise ValueError(n) 

892 

893 result = self 

894 for _ in range(n): 

895 dfuns = np.array([fun.diff() for fun in result]) 

896 result = self.__class__(dfuns) 

897 return result 

898 

899 def conv(self, g: Chebfun) -> Chebfun: 

900 """Compute the convolution of this Chebfun with g. 

901 

902 Computes h(x) = (f ★ g)(x) = ∫ f(t) g(x-t) dt, where domain(f) is 

903 [a, b] and domain(g) is [c, d]. The result is a piecewise Chebfun on 

904 [a + c, b + d] whose breakpoints are the pairwise sums of the 

905 breakpoints of f and g. 

906 

907 Both f and g may be piecewise (contain an arbitrary number of funs). 

908 

909 When both inputs are single-piece with equal-width domains, the fast 

910 Hale-Townsend Legendre convolution algorithm is used. Otherwise, each 

911 output sub-interval is constructed adaptively using Gauss-Legendre 

912 quadrature. 

913 

914 The algorithm is based on: 

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

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

917 

918 Args: 

919 g (Chebfun): A Chebfun (single-piece or piecewise). 

920 

921 Returns: 

922 Chebfun: A piecewise Chebfun on [a + c, b + d] representing 

923 (f ★ g). 

924 

925 Examples: 

926 >>> import numpy as np 

927 >>> from chebpy import chebfun 

928 >>> f = chebfun(lambda x: np.ones_like(x), [-1, 1]) 

929 >>> h = f.conv(f) 

930 >>> bool(abs(h(0.0) - 2.0) < 1e-10) 

931 True 

932 >>> bool(abs(h(-1.0) - 1.0) < 1e-10) 

933 True 

934 >>> bool(abs(h(1.0) - 1.0) < 1e-10) 

935 True 

936 """ 

937 return convolve(self, g) 

938 

939 def sum(self) -> Any: 

940 """Compute the definite integral of the Chebfun over its domain. 

941 

942 This method calculates the definite integral of the Chebfun over its 

943 entire domain of definition by summing the definite integrals of each 

944 piece. 

945 

946 Returns: 

947 float or complex: The definite integral of the Chebfun over its domain. 

948 

949 Examples: 

950 >>> import numpy as np 

951 >>> f = Chebfun.initfun_adaptive(lambda x: x**2, [-1, 1]) 

952 >>> bool(abs(f.sum() - 2.0/3.0) < 1e-10) 

953 True 

954 >>> g = Chebfun.initconst(1.0, [-1, 1]) 

955 >>> bool(abs(g.sum() - 2.0) < 1e-10) 

956 True 

957 """ 

958 return np.sum([fun.sum() for fun in self]) 

959 

960 def dot(self, f: Any) -> Any: 

961 """Compute the dot product of this Chebfun with another function. 

962 

963 This method calculates the inner product (dot product) of this Chebfun 

964 with another function f by multiplying them pointwise and then integrating 

965 the result over the domain. 

966 

967 Args: 

968 f (Chebfun or scalar): The function or scalar to compute the dot product with. 

969 If not a Chebfun, it will be converted to one. 

970 

971 Returns: 

972 float or complex: The dot product of this Chebfun with f. 

973 """ 

974 return (self * f).sum() 

975 

976 def norm(self, p: Any = 2) -> Any: 

977 """Compute the Lp norm of the Chebfun over its domain. 

978 

979 This method calculates the Lp norm of the Chebfun. The L2 norm is the 

980 default and is computed as sqrt(integral(|f|^2)). For p=inf, returns 

981 the maximum absolute value by checking critical points (extrema). 

982 

983 Args: 

984 p (int or float): The norm type. Supported values are 1, 2, positive 

985 integers/floats, or np.inf. Defaults to 2 (L2 norm). 

986 

987 Returns: 

988 float: The Lp norm of the Chebfun. 

989 

990 Examples: 

991 >>> from chebpy import chebfun 

992 >>> import numpy as np 

993 >>> f = chebfun(lambda x: x**2, [-1, 1]) 

994 >>> np.allclose(f.norm(), 0.6324555320336759) # L2 norm 

995 True 

996 >>> np.allclose(f.norm(np.inf), 1.0) # Maximum absolute value 

997 True 

998 """ 

999 if p == 2: 

1000 # L2 norm: sqrt(integral(|f|^2)) 

1001 return np.sqrt(self.dot(self)) 

1002 elif p == np.inf: 

1003 # L-infinity norm: max|f(x)| 

1004 df = self.diff() 

1005 critical_pts = df.roots() 

1006 # Add endpoints 

1007 endpoints = np.array([self.domain[0], self.domain[-1]]) 

1008 # Combine all test points 

1009 test_pts = np.concatenate([critical_pts, endpoints]) 

1010 # Evaluate and find max 

1011 vals = np.abs(self(test_pts)) 

1012 return np.max(vals) 

1013 elif p == 1: 

1014 # L1 norm: integral(|f|) 

1015 return self.absolute().sum() 

1016 elif p > 0: 

1017 # General Lp norm: (integral(|f|^p))^(1/p) 

1018 f_abs = self.absolute() 

1019 f_pow_p = f_abs**p 

1020 integral = f_pow_p.sum() 

1021 return integral ** (1.0 / p) 

1022 else: 

1023 raise ValueError(p) 

1024 

1025 # ---------- 

1026 # utilities 

1027 # ---------- 

1028 @self_empty() 

1029 def absolute(self) -> Chebfun: 

1030 """Absolute value of a Chebfun.""" 

1031 return _pointwise.absolute(self) 

1032 

1033 abs = absolute 

1034 

1035 @self_empty() 

1036 def sign(self) -> Chebfun: 

1037 """Sign function of a Chebfun. 

1038 

1039 Computes the piecewise sign of a Chebfun by finding its roots 

1040 and splitting the domain at those points, then creating constant 

1041 pieces with the appropriate sign values. 

1042 

1043 Returns: 

1044 Chebfun: A new Chebfun representing sign(f(x)). 

1045 """ 

1046 return _pointwise.sign(self) 

1047 

1048 @self_empty() 

1049 def ceil(self) -> Chebfun: 

1050 """Ceiling function of a Chebfun. 

1051 

1052 Computes the piecewise ceiling of a Chebfun by finding where 

1053 the function crosses integer values and splitting the domain 

1054 at those points, then creating constant pieces with the 

1055 appropriate ceiling values. 

1056 

1057 Returns: 

1058 Chebfun: A new Chebfun representing ceil(f(x)). 

1059 """ 

1060 return _pointwise.ceil(self) 

1061 

1062 @self_empty() 

1063 def floor(self) -> Chebfun: 

1064 """Floor function of a Chebfun. 

1065 

1066 Computes the piecewise floor of a Chebfun by finding where 

1067 the function crosses integer values and splitting the domain 

1068 at those points, then creating constant pieces with the 

1069 appropriate floor values. 

1070 

1071 Returns: 

1072 Chebfun: A new Chebfun representing floor(f(x)). 

1073 """ 

1074 return _pointwise.floor(self) 

1075 

1076 @self_empty() 

1077 @cast_arg_to_chebfun 

1078 def maximum(self, other: Any) -> Any: 

1079 """Pointwise maximum of self and another chebfun.""" 

1080 return _pointwise.maximum_minimum(self, other, operator.ge) 

1081 

1082 @self_empty() 

1083 @cast_arg_to_chebfun 

1084 def minimum(self, other: Any) -> Any: 

1085 """Pointwise minimum of self and another chebfun.""" 

1086 return _pointwise.maximum_minimum(self, other, operator.lt) 

1087 

1088 # ---------- 

1089 # plotting 

1090 # ---------- 

1091 def plot(self, ax: Axes | None = None, **kwds: Any) -> Any: 

1092 """Plot the Chebfun over its domain. 

1093 

1094 This method plots the Chebfun over its domain using matplotlib. 

1095 For complex-valued Chebfuns, it plots the real part against the imaginary part. 

1096 

1097 For Chebfuns with ``±inf`` endpoints (containing :class:`CompactFun` 

1098 pieces), each unbounded endpoint is replaced for plotting purposes 

1099 with the corresponding ``plot_support`` endpoint of the outermost 

1100 :class:`CompactFun` piece, so the decay-to-zero region is visible. 

1101 

1102 Args: 

1103 ax (matplotlib.axes.Axes, optional): The axes on which to plot. If None, 

1104 a new axes will be created. Defaults to None. 

1105 **kwds: Additional keyword arguments to pass to matplotlib's plot function. 

1106 

1107 Returns: 

1108 matplotlib.axes.Axes: The axes on which the plot was created. 

1109 """ 

1110 return plot_chebfun(self, ax=ax, **kwds) 

1111 

1112 def plotcoeffs(self, ax: Axes | None = None, **kwds: Any) -> Axes: 

1113 """Plot the coefficients of the Chebfun on a semilogy scale. 

1114 

1115 This method plots the absolute values of the coefficients for each piece 

1116 of the Chebfun on a semilogy scale, which is useful for visualizing the 

1117 decay of coefficients in the Chebyshev series. 

1118 

1119 Args: 

1120 ax (matplotlib.axes.Axes, optional): The axes on which to plot. If None, 

1121 a new axes will be created. Defaults to None. 

1122 **kwds: Additional keyword arguments to pass to matplotlib's semilogy function. 

1123 

1124 Returns: 

1125 matplotlib.axes.Axes: The axes on which the plot was created. 

1126 """ 

1127 return cast(Axes, plotcoeffs_chebfun(self, ax=ax, **kwds)) 

1128 

1129 

1130# --------- 

1131# ufuncs 

1132# --------- 

1133register_ufuncs(Chebfun)