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
« 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.
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.
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"""
12from __future__ import annotations
14import operator
15from collections.abc import Callable, Iterator
16from typing import Any, cast
18import numpy as np
19from matplotlib.axes import Axes
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
30class Chebfun:
31 """Main class for representing and manipulating functions in ChebPy.
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:
37 - Function evaluation at arbitrary points
38 - Algebraic operations (addition, multiplication, etc.)
39 - Calculus operations (differentiation, integration, etc.)
40 - Rootfinding
41 - Plotting
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.
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 """
54 def __init__(self, funs: Any) -> None:
55 """Initialize a Chebfun object.
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
65 @classmethod
66 def initempty(cls) -> Chebfun:
67 """Initialize an empty Chebfun.
69 Returns:
70 Chebfun: An empty Chebfun object with no functions.
72 Examples:
73 >>> f = Chebfun.initempty()
74 >>> f.isempty
75 True
76 """
77 return cls([])
79 @classmethod
80 def initidentity(cls, domain: Any = None) -> Chebfun:
81 """Initialize a Chebfun representing the identity function f(x) = x.
83 Args:
84 domain (array-like, optional): Domain on which to define the identity function.
85 If None, uses the default domain from preferences.
87 Returns:
88 Chebfun: A Chebfun object representing the identity function on the specified domain.
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))
100 @classmethod
101 def initconst(cls, c: Any, domain: Any = None) -> Chebfun:
102 """Initialize a Chebfun representing a constant function f(x) = c.
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.
109 Returns:
110 Chebfun: A Chebfun object representing the constant function on the specified domain.
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}))
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.
135 This method determines the appropriate number of points needed to represent
136 the function to the specified tolerance using an adaptive algorithm.
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).
151 Returns:
152 Chebfun: A Chebfun object representing the function on the specified domain.
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
166 return cls(generate_singular_funs(f, domain, sing=sing, params=params))
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.
172 This method uses a specified number of points to represent the function,
173 rather than determining the number adaptively.
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.
183 Returns:
184 Chebfun: A Chebfun object representing the function on the specified domain.
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)
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.
213 This is a general-purpose constructor that delegates to either initfun_adaptive
214 or initfun_fixedlen based on whether n is provided.
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`.
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)
240 # --------------------
241 # operator overloads
242 # --------------------
243 def __add__(self, f: Any) -> Any:
244 """Add a Chebfun with another Chebfun or a scalar.
246 Args:
247 f (Chebfun or scalar): The object to add to this Chebfun.
249 Returns:
250 Chebfun: A new Chebfun representing the sum.
251 """
252 return self._apply_binop(f, operator.add)
254 @self_empty(np.array([]))
255 @float_argument
256 def __call__(self, x: Any) -> Any:
257 """Evaluate the Chebfun at points x.
259 This method evaluates the Chebfun at the specified points. It handles interior
260 points, breakpoints, and points outside the domain appropriately.
262 Args:
263 x (float or array-like): Points at which to evaluate the Chebfun.
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)
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])
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]
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
290 def __iter__(self) -> Iterator[Any]:
291 """Return an iterator over the functions in this Chebfun.
293 Returns:
294 iterator: An iterator over the functions (funs) in this Chebfun.
295 """
296 return self.funs.__iter__()
298 def __len__(self) -> int:
299 """Return the total number of coefficients across all funs.
301 Returns:
302 int: The sum of sizes of all constituent funs.
303 """
304 return sum(f.size for f in self.funs)
306 def __eq__(self, other: object) -> bool:
307 """Test for equality between two Chebfun objects.
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.
312 Args:
313 other (object): The object to compare with this Chebfun.
315 Returns:
316 bool: True if the objects are equal, False otherwise.
317 """
318 if not isinstance(other, self.__class__):
319 return False
321 # Check if both are empty
322 if self.isempty and other.isempty:
323 return True
325 # Check if domains are equal
326 if self.domain != other.domain:
327 return False
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))
334 def __mul__(self, f: Any) -> Any:
335 """Multiply a Chebfun with another Chebfun or a scalar.
337 Args:
338 f (Chebfun or scalar): The object to multiply with this Chebfun.
340 Returns:
341 Chebfun: A new Chebfun representing the product.
342 """
343 return self._apply_binop(f, operator.mul)
345 def __neg__(self) -> Chebfun:
346 """Return the negative of this Chebfun.
348 Returns:
349 Chebfun: A new Chebfun representing -f(x).
350 """
351 return self.__class__(-self.funs)
353 def __pos__(self) -> Chebfun:
354 """Return the positive of this Chebfun (which is the Chebfun itself).
356 Returns:
357 Chebfun: This Chebfun object (unchanged).
358 """
359 return self
361 def __abs__(self) -> Chebfun:
362 """Return the absolute value of this Chebfun.
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)
372 def __pow__(self, f: Any) -> Any:
373 """Raise this Chebfun to a power.
375 Args:
376 f (Chebfun or scalar): The exponent to which this Chebfun is raised.
378 Returns:
379 Chebfun: A new Chebfun representing self^f.
380 """
381 return self._apply_binop(f, operator.pow)
383 def __rtruediv__(self, c: Any) -> Chebfun:
384 """Divide a scalar by this Chebfun.
386 This method is called when a scalar is divided by a Chebfun, i.e., c / self.
388 Args:
389 c (scalar): The scalar numerator.
391 Returns:
392 Chebfun: A new Chebfun representing c / self.
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 """
399 def constfun(cheb: Any, const: Any) -> Any:
400 return 0.0 * cheb + const
402 def make_divfun(fun: Any) -> Callable[..., Any]:
403 return lambda x: constfun(x, c) / fun(x)
405 newfuns = [fun.initfun_adaptive(make_divfun(fun), fun.interval) for fun in self]
406 return self.__class__(newfuns)
408 @self_empty("Chebfun<empty>")
409 def __repr__(self) -> str:
410 """Return a string representation of the Chebfun.
412 This method returns a detailed string representation of the Chebfun,
413 including information about its domain, intervals, and endpoint values.
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
435 def __rsub__(self, f: Any) -> Any:
436 """Subtract this Chebfun from another object.
438 This method is called when another object is subtracted by this Chebfun,
439 i.e., f - self.
441 Args:
442 f (Chebfun or scalar): The object from which to subtract this Chebfun.
444 Returns:
445 Chebfun: A new Chebfun representing f - self.
446 """
447 return -(self - f)
449 @cast_arg_to_chebfun
450 def __rpow__(self, f: Any) -> Any:
451 """Raise another object to the power of this Chebfun.
453 This method is called when another object is raised to the power of this Chebfun,
454 i.e., f ** self.
456 Args:
457 f (Chebfun or scalar): The base to be raised to the power of this Chebfun.
459 Returns:
460 Chebfun: A new Chebfun representing f ** self.
461 """
462 return f**self
464 def __truediv__(self, f: Any) -> Any:
465 """Divide this Chebfun by another object.
467 Args:
468 f (Chebfun or scalar): The divisor.
470 Returns:
471 Chebfun: A new Chebfun representing self / f.
472 """
473 return self._apply_binop(f, operator.truediv)
475 __rmul__ = __mul__
476 __div__ = __truediv__
477 __rdiv__ = __rtruediv__
478 __radd__ = __add__
480 def __str__(self) -> str:
481 """Return a human-readable string representation of the Chebfun.
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.
487 Returns:
488 str: A detailed string representation of the Chebfun.
489 """
490 return repr(self)
492 def __sub__(self, f: Any) -> Any:
493 """Subtract another object from this Chebfun.
495 Args:
496 f (Chebfun or scalar): The object to subtract from this Chebfun.
498 Returns:
499 Chebfun: A new Chebfun representing self - f.
500 """
501 return self._apply_binop(f, operator.sub)
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.
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.
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).
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)
544 def _break(self, targetdomain: Domain) -> Chebfun:
545 """Resample this Chebfun to a new domain.
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.
551 Args:
552 targetdomain (Domain): The domain to which this Chebfun should be resampled.
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)
570 # ------------
571 # properties
572 # ------------
573 @property
574 def breakpoints(self) -> np.ndarray:
575 """Get the breakpoints of this Chebfun.
577 Breakpoints are the points where the Chebfun transitions from one piece to another.
579 Returns:
580 numpy.ndarray: Array of breakpoints.
581 """
582 return np.array(list(self.breakdata.keys()))
584 @property
585 @self_empty(Domain([]))
586 def domain(self) -> Domain:
587 """Get the domain of this Chebfun.
589 Returns:
590 Domain: A Domain object corresponding to this Chebfun.
591 """
592 return Domain.from_chebfun(self)
594 @domain.setter
595 def domain(self, new_domain: Any) -> None:
596 """Set the domain of the Chebfun by restricting to the new domain.
598 Args:
599 new_domain (array-like): The new domain to which this Chebfun should be restricted.
600 """
601 self.restrict_(new_domain)
603 @property
604 @self_empty(Domain([]))
605 def support(self) -> Any:
606 """Get the support interval of this Chebfun.
608 The support is the interval between the first and last breakpoints.
610 Returns:
611 numpy.ndarray: Array containing the first and last breakpoints.
612 """
613 return self.domain.support
615 @property
616 @self_empty(0.0)
617 def hscale(self) -> float:
618 """Get the horizontal scale of this Chebfun.
620 The horizontal scale is the maximum absolute value of the support interval.
622 Returns:
623 float: The horizontal scale.
624 """
625 return float(np.abs(self.support).max())
627 @property
628 @self_empty(False)
629 def iscomplex(self) -> bool:
630 """Check if this Chebfun has complex values.
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)
638 @property
639 @self_empty(False)
640 def isconst(self) -> bool:
641 """Check if this Chebfun represents a constant function.
643 A Chebfun is constant if all of its pieces are constant with the same value.
645 Returns:
646 bool: True if this Chebfun represents a constant function, False otherwise.
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)
654 @property
655 def isempty(self) -> bool:
656 """Check if this Chebfun is empty.
658 An empty Chebfun contains no functions.
660 Returns:
661 bool: True if this Chebfun is empty, False otherwise.
662 """
663 return self.funs.size == 0
665 @property
666 @self_empty(0.0)
667 def vscale(self) -> Any:
668 """Get the vertical scale of this Chebfun.
670 The vertical scale is the maximum of the vertical scales of all pieces.
672 Returns:
673 float: The vertical scale.
674 """
675 return np.max([fun.vscale for fun in self])
677 @property
678 @self_empty()
679 def x(self) -> Chebfun:
680 """Get the identity function on the support of this Chebfun.
682 This property returns a new Chebfun representing the identity function f(x) = x
683 defined on the same support as this Chebfun.
685 Returns:
686 Chebfun: A Chebfun representing the identity function on the support of this Chebfun.
687 """
688 return self.__class__.initidentity(self.support)
690 # -----------
691 # utilities
692 # ----------
694 def imag(self) -> Chebfun:
695 """Get the imaginary part of this Chebfun.
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)
706 def real(self) -> Chebfun:
707 """Get the real part of this Chebfun.
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
718 def copy(self) -> Chebfun:
719 """Create a deep copy of this Chebfun.
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])
726 @self_empty()
727 def _restrict(self, subinterval: Any) -> Chebfun:
728 """Restrict a Chebfun to a subinterval, without simplifying.
730 This is an internal method that restricts the Chebfun to a subinterval
731 without performing simplification.
733 Args:
734 subinterval (array-like): The subinterval to which this Chebfun should be restricted.
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)
742 def restrict(self, subinterval: Any) -> Any:
743 """Restrict a Chebfun to a subinterval.
745 This method creates a new Chebfun that is restricted to the specified subinterval
746 and simplifies the result.
748 Args:
749 subinterval (array-like): The subinterval to which this Chebfun should be restricted.
751 Returns:
752 Chebfun: A new Chebfun restricted to the specified subinterval.
753 """
754 return self._restrict(subinterval).simplify()
756 @self_empty()
757 def restrict_(self, subinterval: Any) -> Chebfun:
758 """Restrict a Chebfun to a subinterval, modifying the object in place.
760 This method modifies the current Chebfun by restricting it to the specified
761 subinterval and simplifying the result.
763 Args:
764 subinterval (array-like): The subinterval to which this Chebfun should be restricted.
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
774 @cache
775 @self_empty(np.array([]))
776 def roots(self, merge: Any = None) -> np.ndarray:
777 """Compute the roots of a Chebfun.
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.
782 Args:
783 merge (bool, optional): Whether to merge roots at breakpoints. If None,
784 uses the value from preferences. Defaults to None.
786 Returns:
787 numpy.ndarray: Array of roots sorted in ascending order.
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))
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])
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])
821 # ----------
822 # calculus
823 # ----------
824 def cumsum(self) -> Chebfun:
825 """Compute the indefinite integral (antiderivative) of the Chebfun.
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.
832 Returns:
833 Chebfun: A new Chebfun representing the indefinite integral of this Chebfun.
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)
857 def diff(self, n: int = 1) -> Chebfun:
858 """Compute the derivative of the Chebfun.
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.
864 Args:
865 n: Order of differentiation (default: 1). Must be non-negative integer.
867 Returns:
868 Chebfun: A new Chebfun representing the nth derivative of this Chebfun.
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)
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
896 def conv(self, g: Chebfun) -> Chebfun:
897 """Compute the convolution of this Chebfun with g.
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.
904 Both f and g may be piecewise (contain an arbitrary number of funs).
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.
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.
915 Args:
916 g (Chebfun): A Chebfun (single-piece or piecewise).
918 Returns:
919 Chebfun: A piecewise Chebfun on [a + c, b + d] representing
920 (f ★ g).
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
936 return convolve(self, g)
938 def sum(self) -> Any:
939 """Compute the definite integral of the Chebfun over its domain.
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.
945 Returns:
946 float or complex: The definite integral of the Chebfun over its domain.
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])
959 def dot(self, f: Any) -> Any:
960 """Compute the dot product of this Chebfun with another function.
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.
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.
970 Returns:
971 float or complex: The dot product of this Chebfun with f.
972 """
973 return (self * f).sum()
975 def norm(self, p: Any = 2) -> Any:
976 """Compute the Lp norm of the Chebfun over its domain.
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).
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).
986 Returns:
987 float: The Lp norm of the Chebfun.
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)
1024 # ----------
1025 # utilities
1026 # ----------
1027 @self_empty()
1028 def absolute(self) -> Chebfun:
1029 """Absolute value of a Chebfun."""
1030 from ._pointwise import absolute
1032 return absolute(self)
1034 abs = absolute
1036 @self_empty()
1037 def sign(self) -> Chebfun:
1038 """Sign function of a Chebfun.
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.
1044 Returns:
1045 Chebfun: A new Chebfun representing sign(f(x)).
1046 """
1047 from ._pointwise import sign
1049 return sign(self)
1051 @self_empty()
1052 def ceil(self) -> Chebfun:
1053 """Ceiling function of a Chebfun.
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.
1060 Returns:
1061 Chebfun: A new Chebfun representing ceil(f(x)).
1062 """
1063 from ._pointwise import ceil
1065 return ceil(self)
1067 @self_empty()
1068 def floor(self) -> Chebfun:
1069 """Floor function of a Chebfun.
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.
1076 Returns:
1077 Chebfun: A new Chebfun representing floor(f(x)).
1078 """
1079 from ._pointwise import floor
1081 return floor(self)
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
1089 return maximum_minimum(self, other, operator.ge)
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
1097 return maximum_minimum(self, other, operator.lt)
1099 # ----------
1100 # plotting
1101 # ----------
1102 def plot(self, ax: Axes | None = None, **kwds: Any) -> Any:
1103 """Plot the Chebfun over its domain.
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.
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.
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.
1118 Returns:
1119 matplotlib.axes.Axes: The axes on which the plot was created.
1120 """
1121 return plot_chebfun(self, ax=ax, **kwds)
1123 def plotcoeffs(self, ax: Axes | None = None, **kwds: Any) -> Axes:
1124 """Plot the coefficients of the Chebfun on a semilogy scale.
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.
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.
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))
1141# ---------
1142# ufuncs
1143# ---------
1144register_ufuncs(Chebfun)