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

329 statements  

« prev     ^ index     » next       coverage.py v7.15.2, created at 2026-07-22 21:50 +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 ._ufuncs import register_ufuncs 

22from .bndfun import Bndfun 

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

24from .exceptions import BadFunLengthArgument 

25from .plotting import plot_chebfun, plotcoeffs_chebfun 

26from .settings import _preferences as prefs 

27from .utilities import Domain, check_funs, compute_breakdata, generate_funs 

28 

29 

30class Chebfun: 

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

32 

33 The Chebfun class represents functions using piecewise polynomial approximations 

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

35 with these function representations, including: 

36 

37 - Function evaluation at arbitrary points 

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

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

40 - Rootfinding 

41 - Plotting 

42 

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

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

45 approximations, allowing for efficient representation of functions with varying 

46 complexity across different intervals. 

47 

48 Attributes: 

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

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

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

52 """ 

53 

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

55 """Initialize a Chebfun object. 

56 

57 Args: 

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

59 These will be checked and sorted using check_funs. 

60 """ 

61 self.funs = check_funs(funs) 

62 self.breakdata = compute_breakdata(self.funs) 

63 self.transposed = False 

64 

65 @classmethod 

66 def initempty(cls) -> Chebfun: 

67 """Initialize an empty Chebfun. 

68 

69 Returns: 

70 Chebfun: An empty Chebfun object with no functions. 

71 

72 Examples: 

73 >>> f = Chebfun.initempty() 

74 >>> f.isempty 

75 True 

76 """ 

77 return cls([]) 

78 

79 @classmethod 

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

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

82 

83 Args: 

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

85 If None, uses the default domain from preferences. 

86 

87 Returns: 

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

89 

90 Examples: 

91 >>> import numpy as np 

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

93 >>> float(x(0.5)) 

94 0.5 

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

96 True 

97 """ 

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

99 

100 @classmethod 

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

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

103 

104 Args: 

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

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

107 If None, uses the default domain from preferences. 

108 

109 Returns: 

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

111 

112 Examples: 

113 >>> import numpy as np 

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

115 >>> float(f(0)) 

116 3.14 

117 >>> float(f(0.5)) 

118 3.14 

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

120 True 

121 """ 

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

123 

124 @classmethod 

125 def initfun_adaptive( 

126 cls, 

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

128 domain: Any = None, 

129 *, 

130 sing: str | None = None, 

131 params: Any = None, 

132 ) -> Chebfun: 

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

134 

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

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

137 

138 Args: 

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

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

141 If None, uses the default domain from preferences. 

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

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

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

145 using the Adcock-Richardson clustering map; interior pieces 

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

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

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

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

150 

151 Returns: 

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

153 

154 Examples: 

155 >>> import numpy as np 

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

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

158 True 

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

160 True 

161 """ 

162 if sing is None: 

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

164 from ._singular_construction import generate_singular_funs 

165 

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

167 

168 @classmethod 

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

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

171 

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

173 rather than determining the number adaptively. 

174 

175 Args: 

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

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

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

179 the size of the domain. 

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

181 If None, uses the default domain from preferences. 

182 

183 Returns: 

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

185 

186 Raises: 

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

188 """ 

189 nn = np.array(n) 

190 if nn.size < 2: 

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

192 else: 

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

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

195 raise BadFunLengthArgument 

196 funs = [] 

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

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

199 return cls(funs) 

200 

201 @classmethod 

202 def initfun( 

203 cls, 

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

205 domain: Any = None, 

206 n: Any = None, 

207 *, 

208 sing: str | None = None, 

209 params: Any = None, 

210 ) -> Chebfun: 

211 """Initialize a Chebfun from a function. 

212 

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

214 or initfun_fixedlen based on whether n is provided. 

215 

216 Args: 

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

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

219 If None, uses the default domain from preferences. 

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

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

222 sing: Optional endpoint-singularity hint forwarded to 

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

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

225 Forwarded to :meth:`initfun_adaptive`. 

226 

227 Returns: 

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

229 """ 

230 if n is None: 

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

232 if sing is not None: 

233 msg = ( 

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

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

236 ) 

237 raise NotImplementedError(msg) 

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

239 

240 # -------------------- 

241 # operator overloads 

242 # -------------------- 

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

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

245 

246 Args: 

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

248 

249 Returns: 

250 Chebfun: A new Chebfun representing the sum. 

251 """ 

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

253 

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

255 @float_argument 

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

257 """Evaluate the Chebfun at points x. 

258 

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

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

261 

262 Args: 

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

264 

265 Returns: 

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

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

268 """ 

269 # initialise output 

270 dtype = complex if self.iscomplex else float 

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

272 

273 # evaluate a fun when x is an interior point 

274 for fun in self: 

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

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

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

278 

279 # evaluate the breakpoint data for x at a breakpoint 

280 breakpoints = self.breakpoints 

281 for break_point in breakpoints: 

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

283 

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

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

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

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

288 return out 

289 

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

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

292 

293 Returns: 

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

295 """ 

296 return self.funs.__iter__() 

297 

298 def __len__(self) -> int: 

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

300 

301 Returns: 

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

303 """ 

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

305 

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

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

308 

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

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

311 

312 Args: 

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

314 

315 Returns: 

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

317 """ 

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

319 return False 

320 

321 # Check if both are empty 

322 if self.isempty and other.isempty: 

323 return True 

324 

325 # Check if domains are equal 

326 if self.domain != other.domain: 

327 return False 

328 

329 # Check function values at test points 

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

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

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

333 

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

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

336 

337 Args: 

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

339 

340 Returns: 

341 Chebfun: A new Chebfun representing the product. 

342 """ 

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

344 

345 def __neg__(self) -> Chebfun: 

346 """Return the negative of this Chebfun. 

347 

348 Returns: 

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

350 """ 

351 return self.__class__(-self.funs) 

352 

353 def __pos__(self) -> Chebfun: 

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

355 

356 Returns: 

357 Chebfun: This Chebfun object (unchanged). 

358 """ 

359 return self 

360 

361 def __abs__(self) -> Chebfun: 

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

363 

364 Returns: 

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

366 """ 

367 abs_funs = [] 

368 for fun in self.funs: 

369 abs_funs.append(fun.absolute()) 

370 return self.__class__(abs_funs) 

371 

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

373 """Raise this Chebfun to a power. 

374 

375 Args: 

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

377 

378 Returns: 

379 Chebfun: A new Chebfun representing self^f. 

380 """ 

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

382 

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

384 """Divide a scalar by this Chebfun. 

385 

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

387 

388 Args: 

389 c (scalar): The scalar numerator. 

390 

391 Returns: 

392 Chebfun: A new Chebfun representing c / self. 

393 

394 Note: 

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

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

397 """ 

398 

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

400 return 0.0 * cheb + const 

401 

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

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

404 

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

406 return self.__class__(newfuns) 

407 

408 @self_empty("Chebfun<empty>") 

409 def __repr__(self) -> str: 

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

411 

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

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

414 

415 Returns: 

416 str: A string representation of the Chebfun. 

417 """ 

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

419 numpcs = self.funs.size 

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

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

422 toprow = " interval length endpoint values\n" 

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

424 rowdta = "" 

425 for fun in self: 

426 endpts = fun.support 

427 xl, xr = endpts 

428 fl, fr = fun(endpts) 

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

430 rowdta += row 

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

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

433 return header + toprow + rowdta + btmrow + btmxtr 

434 

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

436 """Subtract this Chebfun from another object. 

437 

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

439 i.e., f - self. 

440 

441 Args: 

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

443 

444 Returns: 

445 Chebfun: A new Chebfun representing f - self. 

446 """ 

447 return -(self - f) 

448 

449 @cast_arg_to_chebfun 

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

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

452 

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

454 i.e., f ** self. 

455 

456 Args: 

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

458 

459 Returns: 

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

461 """ 

462 return f**self 

463 

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

465 """Divide this Chebfun by another object. 

466 

467 Args: 

468 f (Chebfun or scalar): The divisor. 

469 

470 Returns: 

471 Chebfun: A new Chebfun representing self / f. 

472 """ 

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

474 

475 __rmul__ = __mul__ 

476 __div__ = __truediv__ 

477 __rdiv__ = __rtruediv__ 

478 __radd__ = __add__ 

479 

480 def __str__(self) -> str: 

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

482 

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

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

485 with the behaviour of numpy and pandas objects. 

486 

487 Returns: 

488 str: A detailed string representation of the Chebfun. 

489 """ 

490 return repr(self) 

491 

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

493 """Subtract another object from this Chebfun. 

494 

495 Args: 

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

497 

498 Returns: 

499 Chebfun: A new Chebfun representing self - f. 

500 """ 

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

502 

503 # ------------------ 

504 # internal helpers 

505 # ------------------ 

506 @self_empty() 

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

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

509 

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

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

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

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

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

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

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

517 

518 Args: 

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

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

521 

522 Returns: 

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

524 """ 

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

526 return f 

527 if np.isscalar(f): 

528 chbfn1 = self 

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

530 simplify = False 

531 else: 

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

533 chbfn1 = self._break(newdom) 

534 chbfn2 = f._break(newdom) 

535 simplify = True 

536 newfuns = [] 

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

538 newfun = op(fun1, fun2) 

539 if simplify: 

540 newfun = newfun.simplify() 

541 newfuns.append(newfun) 

542 return self.__class__(newfuns) 

543 

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

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

546 

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

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

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

550 

551 Args: 

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

553 

554 Returns: 

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

556 """ 

557 newfuns = [] 

558 subintervals = iter(targetdomain.intervals) 

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

560 for fun in self: 

561 while interval in fun.interval: 

562 newfun = fun.restrict(interval) 

563 newfuns.append(newfun) 

564 try: 

565 interval = next(subintervals) 

566 except StopIteration: 

567 break 

568 return self.__class__(newfuns) 

569 

570 # ------------ 

571 # properties 

572 # ------------ 

573 @property 

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

575 """Get the breakpoints of this Chebfun. 

576 

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

578 

579 Returns: 

580 numpy.ndarray: Array of breakpoints. 

581 """ 

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

583 

584 @property 

585 @self_empty(Domain([])) 

586 def domain(self) -> Domain: 

587 """Get the domain of this Chebfun. 

588 

589 Returns: 

590 Domain: A Domain object corresponding to this Chebfun. 

591 """ 

592 return Domain.from_chebfun(self) 

593 

594 @domain.setter 

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

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

597 

598 Args: 

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

600 """ 

601 self.restrict_(new_domain) 

602 

603 @property 

604 @self_empty(Domain([])) 

605 def support(self) -> Any: 

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

607 

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

609 

610 Returns: 

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

612 """ 

613 return self.domain.support 

614 

615 @property 

616 @self_empty(0.0) 

617 def hscale(self) -> float: 

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

619 

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

621 

622 Returns: 

623 float: The horizontal scale. 

624 """ 

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

626 

627 @property 

628 @self_empty(False) 

629 def iscomplex(self) -> bool: 

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

631 

632 Returns: 

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

634 False otherwise. 

635 """ 

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

637 

638 @property 

639 @self_empty(False) 

640 def isconst(self) -> bool: 

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

642 

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

644 

645 Returns: 

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

647 

648 Note: 

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

650 """ 

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

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

653 

654 @property 

655 def isempty(self) -> bool: 

656 """Check if this Chebfun is empty. 

657 

658 An empty Chebfun contains no functions. 

659 

660 Returns: 

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

662 """ 

663 return self.funs.size == 0 

664 

665 @property 

666 @self_empty(0.0) 

667 def vscale(self) -> Any: 

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

669 

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

671 

672 Returns: 

673 float: The vertical scale. 

674 """ 

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

676 

677 @property 

678 @self_empty() 

679 def x(self) -> Chebfun: 

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

681 

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

683 defined on the same support as this Chebfun. 

684 

685 Returns: 

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

687 """ 

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

689 

690 # ----------- 

691 # utilities 

692 # ---------- 

693 

694 def imag(self) -> Chebfun: 

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

696 

697 Returns: 

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

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

700 """ 

701 if self.iscomplex: 

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

703 else: 

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

705 

706 def real(self) -> Chebfun: 

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

708 

709 Returns: 

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

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

712 """ 

713 if self.iscomplex: 

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

715 else: 

716 return self 

717 

718 def copy(self) -> Chebfun: 

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

720 

721 Returns: 

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

723 """ 

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

725 

726 @self_empty() 

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

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

729 

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

731 without performing simplification. 

732 

733 Args: 

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

735 

736 Returns: 

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

738 """ 

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

740 return self._break(newdom) 

741 

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

743 """Restrict a Chebfun to a subinterval. 

744 

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

746 and simplifies the result. 

747 

748 Args: 

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

750 

751 Returns: 

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

753 """ 

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

755 

756 @self_empty() 

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

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

759 

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

761 subinterval and simplifying the result. 

762 

763 Args: 

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

765 

766 Returns: 

767 Chebfun: The modified Chebfun (self). 

768 """ 

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

770 self.funs = restricted.funs 

771 self.breakdata = compute_breakdata(self.funs) 

772 return self 

773 

774 @cache 

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

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

777 """Compute the roots of a Chebfun. 

778 

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

780 of each piece of the Chebfun and combining them. 

781 

782 Args: 

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

784 uses the value from preferences. Defaults to None. 

785 

786 Returns: 

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

788 

789 Examples: 

790 >>> import numpy as np 

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

792 >>> roots = f.roots() 

793 >>> len(roots) 

794 2 

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

796 True 

797 """ 

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

799 allrts = [] 

800 prvrts = np.array([]) 

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

802 for fun in self: 

803 rts = fun.roots() 

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

805 # TODO: there could be multiple roots at breakpoints 

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

807 rts = rts[1:] 

808 allrts.append(rts) 

809 prvrts = rts 

810 return np.concatenate(list(allrts)) 

811 

812 @self_empty() 

813 def simplify(self) -> Chebfun: 

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

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

816 

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

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

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

820 

821 # ---------- 

822 # calculus 

823 # ---------- 

824 def cumsum(self) -> Chebfun: 

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

826 

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

828 constant of integration chosen so that the indefinite integral evaluates 

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

830 are added to ensure continuity across the pieces. 

831 

832 Returns: 

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

834 

835 Examples: 

836 >>> import numpy as np 

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

838 >>> F = f.cumsum() 

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

840 True 

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

842 True 

843 """ 

844 newfuns = [] 

845 prevfun = None 

846 for fun in self: 

847 integral = fun.cumsum() 

848 if prevfun: 

849 # enforce continuity by adding the function value 

850 # at the right endpoint of the previous fun 

851 _, fb = prevfun.endvalues 

852 integral = integral + fb 

853 newfuns.append(integral) 

854 prevfun = integral 

855 return self.__class__(newfuns) 

856 

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

858 """Compute the derivative of the Chebfun. 

859 

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

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

862 corresponding piece in the original Chebfun. 

863 

864 Args: 

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

866 

867 Returns: 

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

869 

870 Examples: 

871 >>> from chebpy import chebfun 

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

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

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

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

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

877 True 

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

879 True 

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

881 True 

882 """ 

883 if not isinstance(n, int): 

884 raise TypeError(n) 

885 if n == 0: 

886 return self 

887 if n < 0: 

888 raise ValueError(n) 

889 

890 result = self 

891 for _ in range(n): 

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

893 result = self.__class__(dfuns) 

894 return result 

895 

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

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

898 

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

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

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

902 breakpoints of f and g. 

903 

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

905 

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

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

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

909 quadrature. 

910 

911 The algorithm is based on: 

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

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

914 

915 Args: 

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

917 

918 Returns: 

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

920 (f ★ g). 

921 

922 Examples: 

923 >>> import numpy as np 

924 >>> from chebpy import chebfun 

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

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

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

928 True 

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

930 True 

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

932 True 

933 """ 

934 from ._convolution import convolve 

935 

936 return convolve(self, g) 

937 

938 def sum(self) -> Any: 

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

940 

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

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

943 piece. 

944 

945 Returns: 

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

947 

948 Examples: 

949 >>> import numpy as np 

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

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

952 True 

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

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

955 True 

956 """ 

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

958 

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

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

961 

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

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

964 the result over the domain. 

965 

966 Args: 

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

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

969 

970 Returns: 

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

972 """ 

973 return (self * f).sum() 

974 

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

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

977 

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

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

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

981 

982 Args: 

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

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

985 

986 Returns: 

987 float: The Lp norm of the Chebfun. 

988 

989 Examples: 

990 >>> from chebpy import chebfun 

991 >>> import numpy as np 

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

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

994 True 

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

996 True 

997 """ 

998 if p == 2: 

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

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

1001 elif p == np.inf: 

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

1003 df = self.diff() 

1004 critical_pts = df.roots() 

1005 # Add endpoints 

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

1007 # Combine all test points 

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

1009 # Evaluate and find max 

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

1011 return np.max(vals) 

1012 elif p == 1: 

1013 # L1 norm: integral(|f|) 

1014 return self.absolute().sum() 

1015 elif p > 0: 

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

1017 f_abs = self.absolute() 

1018 f_pow_p = f_abs**p 

1019 integral = f_pow_p.sum() 

1020 return integral ** (1.0 / p) 

1021 else: 

1022 raise ValueError(p) 

1023 

1024 # ---------- 

1025 # utilities 

1026 # ---------- 

1027 @self_empty() 

1028 def absolute(self) -> Chebfun: 

1029 """Absolute value of a Chebfun.""" 

1030 from ._pointwise import absolute 

1031 

1032 return absolute(self) 

1033 

1034 abs = absolute 

1035 

1036 @self_empty() 

1037 def sign(self) -> Chebfun: 

1038 """Sign function of a Chebfun. 

1039 

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

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

1042 pieces with the appropriate sign values. 

1043 

1044 Returns: 

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

1046 """ 

1047 from ._pointwise import sign 

1048 

1049 return sign(self) 

1050 

1051 @self_empty() 

1052 def ceil(self) -> Chebfun: 

1053 """Ceiling function of a Chebfun. 

1054 

1055 Computes the piecewise ceiling of a Chebfun by finding where 

1056 the function crosses integer values and splitting the domain 

1057 at those points, then creating constant pieces with the 

1058 appropriate ceiling values. 

1059 

1060 Returns: 

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

1062 """ 

1063 from ._pointwise import ceil 

1064 

1065 return ceil(self) 

1066 

1067 @self_empty() 

1068 def floor(self) -> Chebfun: 

1069 """Floor function of a Chebfun. 

1070 

1071 Computes the piecewise floor of a Chebfun by finding where 

1072 the function crosses integer values and splitting the domain 

1073 at those points, then creating constant pieces with the 

1074 appropriate floor values. 

1075 

1076 Returns: 

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

1078 """ 

1079 from ._pointwise import floor 

1080 

1081 return floor(self) 

1082 

1083 @self_empty() 

1084 @cast_arg_to_chebfun 

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

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

1087 from ._pointwise import maximum_minimum 

1088 

1089 return maximum_minimum(self, other, operator.ge) 

1090 

1091 @self_empty() 

1092 @cast_arg_to_chebfun 

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

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

1095 from ._pointwise import maximum_minimum 

1096 

1097 return maximum_minimum(self, other, operator.lt) 

1098 

1099 # ---------- 

1100 # plotting 

1101 # ---------- 

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

1103 """Plot the Chebfun over its domain. 

1104 

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

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

1107 

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

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

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

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

1112 

1113 Args: 

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

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

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

1117 

1118 Returns: 

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

1120 """ 

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

1122 

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

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

1125 

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

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

1128 decay of coefficients in the Chebyshev series. 

1129 

1130 Args: 

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

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

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

1134 

1135 Returns: 

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

1137 """ 

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

1139 

1140 

1141# --------- 

1142# ufuncs 

1143# --------- 

1144register_ufuncs(Chebfun)