Coverage for src / chebpy / chebfun.py: 99%
558 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-05 07:22 +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
18import matplotlib.pyplot as plt
19import numpy as np
20from matplotlib.axes import Axes
22from .algorithms import _conv_legendre, cheb2leg, leg2cheb
23from .bndfun import Bndfun
24from .chebtech import Chebtech
25from .decorators import cache, cast_arg_to_chebfun, float_argument, self_empty
26from .exceptions import BadFunLengthArgument, SupportMismatch
27from .plotting import plotfun
28from .settings import _preferences as prefs
29from .trigtech import Trigtech
30from .utilities import Domain, Interval, check_funs, compute_breakdata, generate_funs
33def _generate_singular_funs(
34 f: Callable[..., Any],
35 domain: Any,
36 *,
37 sing: str,
38 params: Any,
39) -> list[Any]:
40 """Build per-piece funs for a Chebfun with endpoint singularities.
42 The leftmost / rightmost pieces (depending on ``sing``) are built as
43 :class:`~chebpy.singfun.Singfun` instances using the Adcock-Richardson
44 clustering map; all interior pieces are ordinary :class:`Bndfun`.
46 Args:
47 f: Callable evaluating the function in logical coordinates.
48 domain: Breakpoint sequence; outermost endpoints must be finite.
49 sing: One of ``"left"``, ``"right"``, ``"both"``.
50 params: A :class:`~chebpy.maps.MapParams` instance carrying
51 ``(L, alpha)`` for the slit-strip clustering map. ``None`` means
52 use :class:`~chebpy.maps.MapParams` defaults.
54 Returns:
55 list: Per-piece funs ready to feed to :class:`Chebfun`.
57 Raises:
58 ValueError: If ``sing`` is not one of the recognised values.
59 """
60 # Local import: chebfun and singfun are siblings under classicfun and
61 # would otherwise risk a cyclic top-level import.
62 from .maps import MapParams
63 from .singfun import Singfun
65 if sing not in ("left", "right", "both"):
66 msg = f"sing must be 'left', 'right', or 'both'; got {sing!r}"
67 raise ValueError(msg)
68 if params is None:
69 params = MapParams()
70 dom = Domain(domain if domain is not None else prefs.domain)
71 intervals = list(dom.intervals)
72 n_pieces = len(intervals)
73 funs: list[Any] = []
74 for i, interval in enumerate(intervals):
75 is_first = i == 0
76 is_last = i == n_pieces - 1
77 piece_sing: str | None = None
78 if sing == "left" and is_first:
79 piece_sing = "left"
80 elif sing == "right" and is_last:
81 piece_sing = "right"
82 elif sing == "both":
83 if is_first and is_last:
84 piece_sing = "both"
85 elif is_first:
86 piece_sing = "left"
87 elif is_last:
88 piece_sing = "right"
89 if piece_sing is None:
90 funs.append(Bndfun.initfun_adaptive(f, interval))
91 else:
92 funs.append(Singfun.initfun_adaptive(f, interval, sing=piece_sing, params=params))
93 return funs
96class Chebfun:
97 """Main class for representing and manipulating functions in ChebPy.
99 The Chebfun class represents functions using piecewise polynomial approximations
100 on arbitrary intervals. It provides a comprehensive set of operations for working
101 with these function representations, including:
103 - Function evaluation at arbitrary points
104 - Algebraic operations (addition, multiplication, etc.)
105 - Calculus operations (differentiation, integration, etc.)
106 - Rootfinding
107 - Plotting
109 Chebfun objects can be created from callable functions, constant values, or
110 directly from function pieces. The class supports both adaptive and fixed-length
111 approximations, allowing for efficient representation of functions with varying
112 complexity across different intervals.
114 Attributes:
115 funs (numpy.ndarray): Array of function pieces that make up the Chebfun.
116 breakdata (OrderedDict): Mapping of breakpoints to function values.
117 transposed (bool): Flag indicating if the Chebfun is transposed.
118 """
120 def __init__(self, funs: Any) -> None:
121 """Initialize a Chebfun object.
123 Args:
124 funs (list): List of function objects to be included in the Chebfun.
125 These will be checked and sorted using check_funs.
126 """
127 self.funs = check_funs(funs)
128 self.breakdata = compute_breakdata(self.funs)
129 self.transposed = False
131 @classmethod
132 def initempty(cls) -> Chebfun:
133 """Initialize an empty Chebfun.
135 Returns:
136 Chebfun: An empty Chebfun object with no functions.
138 Examples:
139 >>> f = Chebfun.initempty()
140 >>> f.isempty
141 True
142 """
143 return cls([])
145 @classmethod
146 def initidentity(cls, domain: Any = None) -> Chebfun:
147 """Initialize a Chebfun representing the identity function f(x) = x.
149 Args:
150 domain (array-like, optional): Domain on which to define the identity function.
151 If None, uses the default domain from preferences.
153 Returns:
154 Chebfun: A Chebfun object representing the identity function on the specified domain.
156 Examples:
157 >>> import numpy as np
158 >>> x = Chebfun.initidentity([-1, 1])
159 >>> float(x(0.5))
160 0.5
161 >>> np.allclose(x([0, 0.5, 1]), [0, 0.5, 1])
162 True
163 """
164 return cls(generate_funs(domain, Bndfun.initidentity))
166 @classmethod
167 def initconst(cls, c: Any, domain: Any = None) -> Chebfun:
168 """Initialize a Chebfun representing a constant function f(x) = c.
170 Args:
171 c (float or complex): The constant value.
172 domain (array-like, optional): Domain on which to define the constant function.
173 If None, uses the default domain from preferences.
175 Returns:
176 Chebfun: A Chebfun object representing the constant function on the specified domain.
178 Examples:
179 >>> import numpy as np
180 >>> f = Chebfun.initconst(3.14, [-1, 1])
181 >>> float(f(0))
182 3.14
183 >>> float(f(0.5))
184 3.14
185 >>> np.allclose(f([0, 0.5, 1]), [3.14, 3.14, 3.14])
186 True
187 """
188 return cls(generate_funs(domain, Bndfun.initconst, {"c": c}))
190 @classmethod
191 def initfun_adaptive(
192 cls,
193 f: Callable[..., Any],
194 domain: Any = None,
195 *,
196 sing: str | None = None,
197 params: Any = None,
198 ) -> Chebfun:
199 """Initialize a Chebfun by adaptively sampling a function.
201 This method determines the appropriate number of points needed to represent
202 the function to the specified tolerance using an adaptive algorithm.
204 Args:
205 f (callable): The function to be approximated.
206 domain (array-like, optional): Domain on which to define the function.
207 If None, uses the default domain from preferences.
208 sing: Optional endpoint-singularity hint, one of ``"left"``,
209 ``"right"``, or ``"both"``. When set, the appropriate boundary
210 pieces are built as :class:`~chebpy.singfun.Singfun` instances
211 using the Adcock-Richardson clustering map; interior pieces
212 remain :class:`~chebpy.bndfun.Bndfun`.
213 params: Slit-strip map parameters (a :class:`~chebpy.maps.MapParams`).
214 Ignored when ``sing`` is ``None``. Default ``None`` (uses
215 :class:`~chebpy.maps.MapParams` defaults).
217 Returns:
218 Chebfun: A Chebfun object representing the function on the specified domain.
220 Examples:
221 >>> import numpy as np
222 >>> f = Chebfun.initfun_adaptive(lambda x: np.sin(x), [-np.pi, np.pi])
223 >>> bool(abs(f(0)) < 1e-10)
224 True
225 >>> bool(abs(f(np.pi/2) - 1) < 1e-10)
226 True
227 """
228 if sing is None:
229 return cls(generate_funs(domain, Bndfun.initfun_adaptive, {"f": f}))
230 return cls(_generate_singular_funs(f, domain, sing=sing, params=params))
232 @classmethod
233 def initfun_fixedlen(cls, f: Callable[..., Any], n: Any, domain: Any = None) -> Chebfun:
234 """Initialize a Chebfun with a fixed number of points.
236 This method uses a specified number of points to represent the function,
237 rather than determining the number adaptively.
239 Args:
240 f (callable): The function to be approximated.
241 n (int or array-like): Number of points to use. If a single value, uses the same
242 number for each interval. If an array, must have one fewer elements than
243 the size of the domain.
244 domain (array-like, optional): Domain on which to define the function.
245 If None, uses the default domain from preferences.
247 Returns:
248 Chebfun: A Chebfun object representing the function on the specified domain.
250 Raises:
251 BadFunLengthArgument: If n is an array and its size doesn't match domain.size - 1.
252 """
253 nn = np.array(n)
254 if nn.size < 2:
255 funs = generate_funs(domain, Bndfun.initfun_fixedlen, {"f": f, "n": n})
256 else:
257 domain = Domain(domain if domain is not None else prefs.domain)
258 if not nn.size == domain.size - 1:
259 raise BadFunLengthArgument
260 funs = []
261 for interval, length in zip(domain.intervals, nn, strict=False):
262 funs.append(Bndfun.initfun_fixedlen(f, interval, length))
263 return cls(funs)
265 @classmethod
266 def initfun(
267 cls,
268 f: Callable[..., Any],
269 domain: Any = None,
270 n: Any = None,
271 *,
272 sing: str | None = None,
273 params: Any = None,
274 ) -> Chebfun:
275 """Initialize a Chebfun from a function.
277 This is a general-purpose constructor that delegates to either initfun_adaptive
278 or initfun_fixedlen based on whether n is provided.
280 Args:
281 f (callable): The function to be approximated.
282 domain (array-like, optional): Domain on which to define the function.
283 If None, uses the default domain from preferences.
284 n (int or array-like, optional): Number of points to use. If None, determines
285 the number adaptively. If provided, uses a fixed number of points.
286 sing: Optional endpoint-singularity hint forwarded to
287 :meth:`initfun_adaptive`. Only valid when ``n is None``.
288 params: Slit-strip map parameters (a :class:`~chebpy.maps.MapParams`).
289 Forwarded to :meth:`initfun_adaptive`.
291 Returns:
292 Chebfun: A Chebfun object representing the function on the specified domain.
293 """
294 if n is None:
295 return cls.initfun_adaptive(f, domain, sing=sing, params=params)
296 if sing is not None:
297 msg = (
298 "fixed-length construction with sing= is not supported in v1; "
299 "pass n=None for adaptive Singfun construction."
300 )
301 raise NotImplementedError(msg)
302 return cls.initfun_fixedlen(f, n, domain)
304 # --------------------
305 # operator overloads
306 # --------------------
307 def __add__(self, f: Any) -> Any:
308 """Add a Chebfun with another Chebfun or a scalar.
310 Args:
311 f (Chebfun or scalar): The object to add to this Chebfun.
313 Returns:
314 Chebfun: A new Chebfun representing the sum.
315 """
316 return self._apply_binop(f, operator.add)
318 @self_empty(np.array([]))
319 @float_argument
320 def __call__(self, x: Any) -> Any:
321 """Evaluate the Chebfun at points x.
323 This method evaluates the Chebfun at the specified points. It handles interior
324 points, breakpoints, and points outside the domain appropriately.
326 Args:
327 x (float or array-like): Points at which to evaluate the Chebfun.
329 Returns:
330 float or numpy.ndarray: The value(s) of the Chebfun at the specified point(s).
331 Returns a scalar if x is a scalar, otherwise an array of the same size as x.
332 """
333 # initialise output
334 dtype = complex if self.iscomplex else float
335 out = np.full(x.size, np.nan, dtype=dtype)
337 # evaluate a fun when x is an interior point
338 for fun in self:
339 sa, sb = fun.support[0], fun.support[-1]
340 idx = np.logical_and(sa < x, x < sb)
341 out[idx] = fun(x[idx])
343 # evaluate the breakpoint data for x at a breakpoint
344 breakpoints = self.breakpoints
345 for break_point in breakpoints:
346 out[x == break_point] = self.breakdata[break_point]
348 # first and last funs used to evaluate outside of the chebfun domain
349 lpts, rpts = x < breakpoints[0], x > breakpoints[-1]
350 out[lpts] = self.funs[0](x[lpts])
351 out[rpts] = self.funs[-1](x[rpts])
352 return out
354 def __iter__(self) -> Iterator[Any]:
355 """Return an iterator over the functions in this Chebfun.
357 Returns:
358 iterator: An iterator over the functions (funs) in this Chebfun.
359 """
360 return self.funs.__iter__()
362 def __len__(self) -> int:
363 """Return the total number of coefficients across all funs.
365 Returns:
366 int: The sum of sizes of all constituent funs.
367 """
368 return sum(f.size for f in self.funs)
370 def __eq__(self, other: object) -> bool:
371 """Test for equality between two Chebfun objects.
373 Two Chebfun objects are considered equal if they have the same domain
374 and their function values are equal (within tolerance) at a set of test points.
376 Args:
377 other (object): The object to compare with this Chebfun.
379 Returns:
380 bool: True if the objects are equal, False otherwise.
381 """
382 if not isinstance(other, self.__class__):
383 return False
385 # Check if both are empty
386 if self.isempty and other.isempty:
387 return True
389 # Check if domains are equal
390 if self.domain != other.domain:
391 return False
393 # Check function values at test points
394 xx = np.linspace(self.support[0], self.support[1], 100)
395 tol = 1e2 * max(self.vscale, other.vscale) * prefs.eps
396 return bool(np.all(np.abs(self(xx) - other(xx)) <= tol))
398 def __mul__(self, f: Any) -> Any:
399 """Multiply a Chebfun with another Chebfun or a scalar.
401 Args:
402 f (Chebfun or scalar): The object to multiply with this Chebfun.
404 Returns:
405 Chebfun: A new Chebfun representing the product.
406 """
407 return self._apply_binop(f, operator.mul)
409 def __neg__(self) -> Chebfun:
410 """Return the negative of this Chebfun.
412 Returns:
413 Chebfun: A new Chebfun representing -f(x).
414 """
415 return self.__class__(-self.funs)
417 def __pos__(self) -> Chebfun:
418 """Return the positive of this Chebfun (which is the Chebfun itself).
420 Returns:
421 Chebfun: This Chebfun object (unchanged).
422 """
423 return self
425 def __abs__(self) -> Chebfun:
426 """Return the absolute value of this Chebfun.
428 Returns:
429 Chebfun: A new Chebfun representing |f(x)|.
430 """
431 abs_funs = []
432 for fun in self.funs:
433 abs_funs.append(fun.absolute())
434 return self.__class__(abs_funs)
436 def __pow__(self, f: Any) -> Any:
437 """Raise this Chebfun to a power.
439 Args:
440 f (Chebfun or scalar): The exponent to which this Chebfun is raised.
442 Returns:
443 Chebfun: A new Chebfun representing self^f.
444 """
445 return self._apply_binop(f, operator.pow)
447 def __rtruediv__(self, c: Any) -> Chebfun:
448 """Divide a scalar by this Chebfun.
450 This method is called when a scalar is divided by a Chebfun, i.e., c / self.
452 Args:
453 c (scalar): The scalar numerator.
455 Returns:
456 Chebfun: A new Chebfun representing c / self.
458 Note:
459 This is executed when truediv(f, self) fails, which is to say whenever c
460 is not a Chebfun. We proceed on the assumption f is a scalar.
461 """
463 def constfun(cheb: Any, const: Any) -> Any:
464 return 0.0 * cheb + const
466 def make_divfun(fun: Any) -> Callable[..., Any]:
467 return lambda x: constfun(x, c) / fun(x)
469 newfuns = [fun.initfun_adaptive(make_divfun(fun), fun.interval) for fun in self]
470 return self.__class__(newfuns)
472 @self_empty("Chebfun<empty>")
473 def __repr__(self) -> str:
474 """Return a string representation of the Chebfun.
476 This method returns a detailed string representation of the Chebfun,
477 including information about its domain, intervals, and endpoint values.
479 Returns:
480 str: A string representation of the Chebfun.
481 """
482 rowcol = "row" if self.transposed else "column"
483 numpcs = self.funs.size
484 plural = "" if numpcs == 1 else "s"
485 header = f"Chebfun {rowcol} ({numpcs} smooth piece{plural})\n"
486 toprow = " interval length endpoint values\n"
487 tmplat = "[{:8.2g},{:8.2g}] {:6} {:8.2g} {:8.2g}\n"
488 rowdta = ""
489 for fun in self:
490 endpts = fun.support
491 xl, xr = endpts
492 fl, fr = fun(endpts)
493 row = tmplat.format(xl, xr, fun.size, fl, fr)
494 rowdta += row
495 btmrow = f"vertical scale = {self.vscale:3.2g}"
496 btmxtr = "" if numpcs == 1 else f" total length = {sum([f.size for f in self])}"
497 return header + toprow + rowdta + btmrow + btmxtr
499 def __rsub__(self, f: Any) -> Any:
500 """Subtract this Chebfun from another object.
502 This method is called when another object is subtracted by this Chebfun,
503 i.e., f - self.
505 Args:
506 f (Chebfun or scalar): The object from which to subtract this Chebfun.
508 Returns:
509 Chebfun: A new Chebfun representing f - self.
510 """
511 return -(self - f)
513 @cast_arg_to_chebfun
514 def __rpow__(self, f: Any) -> Any:
515 """Raise another object to the power of this Chebfun.
517 This method is called when another object is raised to the power of this Chebfun,
518 i.e., f ** self.
520 Args:
521 f (Chebfun or scalar): The base to be raised to the power of this Chebfun.
523 Returns:
524 Chebfun: A new Chebfun representing f ** self.
525 """
526 return f**self
528 def __truediv__(self, f: Any) -> Any:
529 """Divide this Chebfun by another object.
531 Args:
532 f (Chebfun or scalar): The divisor.
534 Returns:
535 Chebfun: A new Chebfun representing self / f.
536 """
537 return self._apply_binop(f, operator.truediv)
539 __rmul__ = __mul__
540 __div__ = __truediv__
541 __rdiv__ = __rtruediv__
542 __radd__ = __add__
544 def __str__(self) -> str:
545 """Return a human-readable string representation of the Chebfun.
547 This method returns the same detailed representation as ``__repr__``,
548 so that ``print(f)`` shows the full summary table. This is consistent
549 with the behaviour of numpy and pandas objects.
551 Returns:
552 str: A detailed string representation of the Chebfun.
553 """
554 return repr(self)
556 def __sub__(self, f: Any) -> Any:
557 """Subtract another object from this Chebfun.
559 Args:
560 f (Chebfun or scalar): The object to subtract from this Chebfun.
562 Returns:
563 Chebfun: A new Chebfun representing self - f.
564 """
565 return self._apply_binop(f, operator.sub)
567 # ------------------
568 # internal helpers
569 # ------------------
570 @self_empty()
571 def _apply_binop(self, f: Any, op: Callable[..., Any]) -> Any:
572 """Apply a binary operation between this Chebfun and another object.
574 This is a funnel method used in the implementation of Chebfun binary
575 operators. The high-level idea is to first break each chebfun into a
576 series of pieces corresponding to the union of the domains of each
577 before applying the supplied binary operator and simplifying. In the
578 case of the second argument being a scalar we don't need to do the
579 simplify step, since at the Tech-level these operations are defined
580 such that there is no change in the number of coefficients.
582 Args:
583 f (Chebfun or scalar): The second operand of the binary operation.
584 op (callable): The binary operation to apply (e.g., operator.add).
586 Returns:
587 Chebfun: A new Chebfun resulting from applying the binary operation.
588 """
589 if hasattr(f, "isempty") and f.isempty:
590 return f
591 if np.isscalar(f):
592 chbfn1 = self
593 chbfn2 = f * np.ones(self.funs.size)
594 simplify = False
595 else:
596 newdom = self.domain.union(f.domain)
597 chbfn1 = self._break(newdom)
598 chbfn2 = f._break(newdom)
599 simplify = True
600 newfuns = []
601 for fun1, fun2 in zip(chbfn1, chbfn2, strict=False):
602 newfun = op(fun1, fun2)
603 if simplify:
604 newfun = newfun.simplify()
605 newfuns.append(newfun)
606 return self.__class__(newfuns)
608 def _break(self, targetdomain: Domain) -> Chebfun:
609 """Resample this Chebfun to a new domain.
611 This method resamples the Chebfun to the supplied Domain object. It is
612 intended as a private method since one will typically need to have
613 called either Domain.union(f) or Domain.merge(f) prior to calling this method.
615 Args:
616 targetdomain (Domain): The domain to which this Chebfun should be resampled.
618 Returns:
619 Chebfun: A new Chebfun resampled to the target domain.
620 """
621 newfuns = []
622 subintervals = iter(targetdomain.intervals)
623 interval = next(subintervals) # next(..) for Python2/3 compatibility
624 for fun in self:
625 while interval in fun.interval:
626 newfun = fun.restrict(interval)
627 newfuns.append(newfun)
628 try:
629 interval = next(subintervals)
630 except StopIteration:
631 break
632 return self.__class__(newfuns)
634 # ------------
635 # properties
636 # ------------
637 @property
638 def breakpoints(self) -> np.ndarray:
639 """Get the breakpoints of this Chebfun.
641 Breakpoints are the points where the Chebfun transitions from one piece to another.
643 Returns:
644 numpy.ndarray: Array of breakpoints.
645 """
646 return np.array(list(self.breakdata.keys()))
648 @property
649 @self_empty(Domain([]))
650 def domain(self) -> Domain:
651 """Get the domain of this Chebfun.
653 Returns:
654 Domain: A Domain object corresponding to this Chebfun.
655 """
656 return Domain.from_chebfun(self)
658 @domain.setter
659 def domain(self, new_domain: Any) -> None:
660 """Set the domain of the Chebfun by restricting to the new domain.
662 Args:
663 new_domain (array-like): The new domain to which this Chebfun should be restricted.
664 """
665 self.restrict_(new_domain)
667 @property
668 @self_empty(Domain([]))
669 def support(self) -> Any:
670 """Get the support interval of this Chebfun.
672 The support is the interval between the first and last breakpoints.
674 Returns:
675 numpy.ndarray: Array containing the first and last breakpoints.
676 """
677 return self.domain.support
679 @property
680 @self_empty(0.0)
681 def hscale(self) -> float:
682 """Get the horizontal scale of this Chebfun.
684 The horizontal scale is the maximum absolute value of the support interval.
686 Returns:
687 float: The horizontal scale.
688 """
689 return float(np.abs(self.support).max())
691 @property
692 @self_empty(False)
693 def iscomplex(self) -> bool:
694 """Check if this Chebfun has complex values.
696 Returns:
697 bool: True if any of the functions in this Chebfun have complex values,
698 False otherwise.
699 """
700 return any(fun.iscomplex for fun in self)
702 @property
703 @self_empty(False)
704 def isconst(self) -> bool:
705 """Check if this Chebfun represents a constant function.
707 A Chebfun is constant if all of its pieces are constant with the same value.
709 Returns:
710 bool: True if this Chebfun represents a constant function, False otherwise.
712 Note:
713 TODO: find an abstract way of referencing funs[0].coeffs[0]
714 """
715 c = self.funs[0].coeffs[0]
716 return all(fun.isconst and fun.coeffs[0] == c for fun in self)
718 @property
719 def isempty(self) -> bool:
720 """Check if this Chebfun is empty.
722 An empty Chebfun contains no functions.
724 Returns:
725 bool: True if this Chebfun is empty, False otherwise.
726 """
727 return self.funs.size == 0
729 @property
730 @self_empty(0.0)
731 def vscale(self) -> Any:
732 """Get the vertical scale of this Chebfun.
734 The vertical scale is the maximum of the vertical scales of all pieces.
736 Returns:
737 float: The vertical scale.
738 """
739 return np.max([fun.vscale for fun in self])
741 @property
742 @self_empty()
743 def x(self) -> Chebfun:
744 """Get the identity function on the support of this Chebfun.
746 This property returns a new Chebfun representing the identity function f(x) = x
747 defined on the same support as this Chebfun.
749 Returns:
750 Chebfun: A Chebfun representing the identity function on the support of this Chebfun.
751 """
752 return self.__class__.initidentity(self.support)
754 # -----------
755 # utilities
756 # ----------
758 def imag(self) -> Chebfun:
759 """Get the imaginary part of this Chebfun.
761 Returns:
762 Chebfun: A new Chebfun representing the imaginary part of this Chebfun.
763 If this Chebfun is real-valued, returns a zero Chebfun.
764 """
765 if self.iscomplex:
766 return self.__class__([fun.imag() for fun in self])
767 else:
768 return self.initconst(0, domain=self.domain)
770 def real(self) -> Chebfun:
771 """Get the real part of this Chebfun.
773 Returns:
774 Chebfun: A new Chebfun representing the real part of this Chebfun.
775 If this Chebfun is already real-valued, returns this Chebfun.
776 """
777 if self.iscomplex:
778 return self.__class__([fun.real() for fun in self])
779 else:
780 return self
782 def copy(self) -> Chebfun:
783 """Create a deep copy of this Chebfun.
785 Returns:
786 Chebfun: A new Chebfun that is a deep copy of this Chebfun.
787 """
788 return self.__class__([fun.copy() for fun in self])
790 @self_empty()
791 def _restrict(self, subinterval: Any) -> Chebfun:
792 """Restrict a Chebfun to a subinterval, without simplifying.
794 This is an internal method that restricts the Chebfun to a subinterval
795 without performing simplification.
797 Args:
798 subinterval (array-like): The subinterval to which this Chebfun should be restricted.
800 Returns:
801 Chebfun: A new Chebfun restricted to the specified subinterval, without simplification.
802 """
803 newdom = self.domain.restrict(Domain(subinterval))
804 return self._break(newdom)
806 def restrict(self, subinterval: Any) -> Any:
807 """Restrict a Chebfun to a subinterval.
809 This method creates a new Chebfun that is restricted to the specified subinterval
810 and simplifies the result.
812 Args:
813 subinterval (array-like): The subinterval to which this Chebfun should be restricted.
815 Returns:
816 Chebfun: A new Chebfun restricted to the specified subinterval.
817 """
818 return self._restrict(subinterval).simplify()
820 @self_empty()
821 def restrict_(self, subinterval: Any) -> Chebfun:
822 """Restrict a Chebfun to a subinterval, modifying the object in place.
824 This method modifies the current Chebfun by restricting it to the specified
825 subinterval and simplifying the result.
827 Args:
828 subinterval (array-like): The subinterval to which this Chebfun should be restricted.
830 Returns:
831 Chebfun: The modified Chebfun (self).
832 """
833 restricted = self._restrict(subinterval).simplify()
834 self.funs = restricted.funs
835 self.breakdata = compute_breakdata(self.funs)
836 return self
838 @cache
839 @self_empty(np.array([]))
840 def roots(self, merge: Any = None) -> np.ndarray:
841 """Compute the roots of a Chebfun.
843 This method finds the values x for which f(x) = 0, by computing the roots
844 of each piece of the Chebfun and combining them.
846 Args:
847 merge (bool, optional): Whether to merge roots at breakpoints. If None,
848 uses the value from preferences. Defaults to None.
850 Returns:
851 numpy.ndarray: Array of roots sorted in ascending order.
853 Examples:
854 >>> import numpy as np
855 >>> f = Chebfun.initfun_adaptive(lambda x: x**2 - 1, [-2, 2])
856 >>> roots = f.roots()
857 >>> len(roots)
858 2
859 >>> np.allclose(sorted(roots), [-1, 1])
860 True
861 """
862 merge = merge if merge is not None else prefs.mergeroots
863 allrts = []
864 prvrts = np.array([])
865 htol = 1e2 * self.hscale * prefs.eps
866 for fun in self:
867 rts = fun.roots()
868 # ignore first root if equal to the last root of previous fun
869 # TODO: there could be multiple roots at breakpoints
870 if prvrts.size > 0 and rts.size > 0 and merge and abs(prvrts[-1] - rts[0]) <= htol:
871 rts = rts[1:]
872 allrts.append(rts)
873 prvrts = rts
874 return np.concatenate(list(allrts))
876 @self_empty()
877 def simplify(self) -> Chebfun:
878 """Simplify each fun in the chebfun."""
879 return self.__class__([fun.simplify() for fun in self])
881 def translate(self, c: Any) -> Chebfun:
882 """Translate a chebfun by c, i.e., return f(x-c)."""
883 return self.__class__([x.translate(c) for x in self])
885 # ----------
886 # calculus
887 # ----------
888 def cumsum(self) -> Chebfun:
889 """Compute the indefinite integral (antiderivative) of the Chebfun.
891 This method computes the indefinite integral of the Chebfun, with the
892 constant of integration chosen so that the indefinite integral evaluates
893 to 0 at the left endpoint of the domain. For piecewise functions, constants
894 are added to ensure continuity across the pieces.
896 Returns:
897 Chebfun: A new Chebfun representing the indefinite integral of this Chebfun.
899 Examples:
900 >>> import numpy as np
901 >>> f = Chebfun.initconst(1.0, [-1, 1])
902 >>> F = f.cumsum()
903 >>> bool(abs(F(-1)) < 1e-10)
904 True
905 >>> bool(abs(F(1) - 2.0) < 1e-10)
906 True
907 """
908 newfuns = []
909 prevfun = None
910 for fun in self:
911 integral = fun.cumsum()
912 if prevfun:
913 # enforce continuity by adding the function value
914 # at the right endpoint of the previous fun
915 _, fb = prevfun.endvalues
916 integral = integral + fb
917 newfuns.append(integral)
918 prevfun = integral
919 return self.__class__(newfuns)
921 def diff(self, n: int = 1) -> Chebfun:
922 """Compute the derivative of the Chebfun.
924 This method calculates the nth derivative of the Chebfun with respect to x.
925 It creates a new Chebfun where each piece is the derivative of the
926 corresponding piece in the original Chebfun.
928 Args:
929 n: Order of differentiation (default: 1). Must be non-negative integer.
931 Returns:
932 Chebfun: A new Chebfun representing the nth derivative of this Chebfun.
934 Examples:
935 >>> from chebpy import chebfun
936 >>> f = chebfun(lambda x: x**3)
937 >>> df1 = f.diff() # first derivative: 3*x**2
938 >>> df2 = f.diff(2) # second derivative: 6*x
939 >>> df3 = f.diff(3) # third derivative: 6
940 >>> bool(abs(df1(0.5) - 0.75) < 1e-10)
941 True
942 >>> bool(abs(df2(0.5) - 3.0) < 1e-10)
943 True
944 >>> bool(abs(df3(0.5) - 6.0) < 1e-10)
945 True
946 """
947 if not isinstance(n, int):
948 raise TypeError(n)
949 if n == 0:
950 return self
951 if n < 0:
952 raise ValueError(n)
954 result = self
955 for _ in range(n):
956 dfuns = np.array([fun.diff() for fun in result])
957 result = self.__class__(dfuns)
958 return result
960 def conv(self, g: Chebfun) -> Chebfun:
961 """Compute the convolution of this Chebfun with g.
963 Computes h(x) = (f ★ g)(x) = ∫ f(t) g(x-t) dt, where domain(f) is
964 [a, b] and domain(g) is [c, d]. The result is a piecewise Chebfun on
965 [a + c, b + d] whose breakpoints are the pairwise sums of the
966 breakpoints of f and g.
968 Both f and g may be piecewise (contain an arbitrary number of funs).
970 When both inputs are single-piece with equal-width domains, the fast
971 Hale-Townsend Legendre convolution algorithm is used. Otherwise, each
972 output sub-interval is constructed adaptively using Gauss-Legendre
973 quadrature.
975 The algorithm is based on:
976 N. Hale and A. Townsend, "An algorithm for the convolution of
977 Legendre series", SIAM J. Sci. Comput., 36(3), A1207-A1220, 2014.
979 Args:
980 g (Chebfun): A Chebfun (single-piece or piecewise).
982 Returns:
983 Chebfun: A piecewise Chebfun on [a + c, b + d] representing
984 (f ★ g).
986 Examples:
987 >>> import numpy as np
988 >>> from chebpy import chebfun
989 >>> f = chebfun(lambda x: np.ones_like(x), [-1, 1])
990 >>> h = f.conv(f)
991 >>> bool(abs(h(0.0) - 2.0) < 1e-10)
992 True
993 >>> bool(abs(h(-1.0) - 1.0) < 1e-10)
994 True
995 >>> bool(abs(h(1.0) - 1.0) < 1e-10)
996 True
997 """
998 if self.isempty or g.isempty:
999 return self.__class__.initempty()
1001 # Trigtech inputs are not supported: the underlying algorithms assume
1002 # Chebyshev coefficients and would silently produce incorrect results
1003 # if applied to Fourier coefficients. See chebfun's @trigtech/conv.m
1004 # for the same restriction; periodic functions need a circular
1005 # convolution (circconv) instead, which has different semantics
1006 # (fixed period, no domain expansion).
1007 if any(isinstance(fun.onefun, Trigtech) for fun in self.funs) or any(
1008 isinstance(fun.onefun, Trigtech) for fun in g.funs
1009 ):
1010 raise NotImplementedError(
1011 "conv() is not supported for trigfun (Trigtech-backed) inputs. "
1012 "Aperiodic convolution and periodic (circular) convolution are "
1013 "distinct operations; a dedicated circconv() for trigfuns is "
1014 "not yet implemented."
1015 )
1017 # Singfun inputs are not supported: the Hale-Townsend Legendre algorithm
1018 # and the Gauss-Legendre fallback both assume an affine map between the
1019 # logical and reference variables, which the Adcock-Richardson clustering
1020 # map breaks.
1021 from .singfun import Singfun
1023 if any(isinstance(fun, Singfun) for fun in self.funs) or any(isinstance(fun, Singfun) for fun in g.funs):
1024 raise NotImplementedError(
1025 "conv() is not supported for Chebfuns containing Singfun pieces "
1026 "(functions with endpoint singularities represented by a non-affine "
1027 "clustering map)."
1028 )
1030 # Fast path: both single-piece with equal-width finite domains.
1031 # CompactFun pieces (with possibly infinite logical support) always
1032 # take the general piecewise path so the output is wrapped correctly.
1033 from .compactfun import CompactFun
1034 from .exceptions import DivergentIntegralError
1036 # Convolution of a function with non-zero asymptotic limits diverges
1037 # on an unbounded interval, so refuse early with a clear error
1038 # pointing the user at the algebraic-closure escape hatch.
1039 for label, h in (("self", self), ("other", g)):
1040 for piece in h.funs:
1041 if isinstance(piece, CompactFun) and (piece.tail_left != 0.0 or piece.tail_right != 0.0):
1042 raise DivergentIntegralError( # noqa: TRY003
1043 f"Convolution requires both operands to decay to zero at "
1044 f"±inf; got tail_left={piece.tail_left}, "
1045 f"tail_right={piece.tail_right} for {label}. Consider "
1046 f"subtracting a matched sigmoid first so the residual "
1047 f"has zero tails, then convolving the residual."
1048 )
1050 any_compact = any(isinstance(fun, CompactFun) for fun in self.funs) or any(
1051 isinstance(fun, CompactFun) for fun in g.funs
1052 )
1053 if not any_compact and self.funs.size == 1 and g.funs.size == 1:
1054 f_fun, g_fun = self.funs[0], g.funs[0]
1055 f_w = float(f_fun.support[1]) - float(f_fun.support[0])
1056 g_w = float(g_fun.support[1]) - float(g_fun.support[0])
1057 if np.isclose(f_w, g_w):
1058 return self._conv_equal_width_pair(f_fun, g_fun)
1060 # General piecewise convolution
1061 return self._conv_piecewise(g)
1063 def _conv_equal_width_pair(self, f_fun: Any, g_fun: Any) -> Chebfun:
1064 """Convolve two single Bndfuns of equal width using the fast algorithm.
1066 Uses the Hale-Townsend Legendre convolution. The two funs may be on
1067 different intervals as long as they have the same width.
1068 """
1069 a = float(f_fun.support[0])
1070 b = float(f_fun.support[1])
1071 c = float(g_fun.support[0])
1072 d = float(g_fun.support[1])
1074 h = (b - a) / 2.0 # half-width (same for both funs)
1076 leg_f = cheb2leg(f_fun.coeffs)
1077 leg_g = cheb2leg(g_fun.coeffs)
1079 gamma_left, gamma_right = _conv_legendre(leg_f, leg_g)
1081 gamma_left = h * gamma_left
1082 gamma_right = h * gamma_right
1084 cheb_left = leg2cheb(gamma_left)
1085 cheb_right = leg2cheb(gamma_right)
1087 mid = (a + b + c + d) / 2.0
1088 left_interval = Interval(a + c, mid)
1089 right_interval = Interval(mid, b + d)
1091 left_fun = Bndfun(Chebtech(cheb_left), left_interval)
1092 right_fun = Bndfun(Chebtech(cheb_right), right_interval)
1094 return self.__class__([left_fun, right_fun])
1096 def _conv_piecewise(self, g: Chebfun) -> Chebfun:
1097 """General piecewise convolution via Gauss-Legendre quadrature.
1099 The breakpoints of the result are the sorted, unique pairwise sums of
1100 the breakpoints of self and g. On each sub-interval the convolution
1101 integral is smooth, so we construct it adaptively. When either input
1102 contains :class:`CompactFun` pieces, the corresponding ``±inf``
1103 breakpoints are replaced with the numerical-support bounds for the
1104 purposes of integration; the outermost output pieces are then wrapped
1105 as :class:`CompactFun` so the result preserves the unbounded logical
1106 support.
1107 """
1108 from .compactfun import CompactFun
1110 def _effective_breakpoints(h: Chebfun) -> np.ndarray:
1111 """Return breakpoints with ±inf replaced by numerical_support bounds."""
1112 bps = np.array(h.breakpoints, dtype=float)
1113 if not np.isfinite(bps[0]) and isinstance(h.funs[0], CompactFun):
1114 bps[0] = float(h.funs[0].numerical_support[0])
1115 if not np.isfinite(bps[-1]) and isinstance(h.funs[-1], CompactFun):
1116 bps[-1] = float(h.funs[-1].numerical_support[1])
1117 return bps
1119 f_logical_breaks = np.array(self.breakpoints, dtype=float)
1120 g_logical_breaks = np.array(g.breakpoints, dtype=float)
1121 left_inf = (not np.isfinite(f_logical_breaks[0])) or (not np.isfinite(g_logical_breaks[0]))
1122 right_inf = (not np.isfinite(f_logical_breaks[-1])) or (not np.isfinite(g_logical_breaks[-1]))
1124 f_breaks = _effective_breakpoints(self)
1125 g_breaks = _effective_breakpoints(g)
1126 f_a, f_b = float(f_breaks[0]), float(f_breaks[-1])
1127 g_c, g_d = float(g_breaks[0]), float(g_breaks[-1])
1129 # Output breakpoints: all pairwise sums, uniquified and coalesced
1130 out_breaks = np.unique(np.add.outer(f_breaks, g_breaks).ravel())
1131 hscl = max(abs(out_breaks[0]), abs(out_breaks[-1]), 1.0)
1132 tol = 10.0 * np.finfo(float).eps * hscl
1133 mask = np.concatenate(([True], np.diff(out_breaks) > tol))
1134 out_breaks = out_breaks[mask]
1136 # Quadrature order: sufficient for exact integration of polynomial
1137 # integrand on each smooth sub-interval
1138 max_deg = max(fun.size for fun in self.funs) + max(fun.size for fun in g.funs)
1139 n_quad = max(int(np.ceil((max_deg + 1) / 2)), 16)
1140 quad_nodes, quad_weights = np.polynomial.legendre.leggauss(n_quad)
1142 # Pre-convert breakpoints to plain float lists for the inner loop
1143 f_bps = [float(bp) for bp in f_breaks]
1144 g_bps = [float(bp) for bp in g_breaks]
1146 def conv_eval(x: np.ndarray) -> np.ndarray:
1147 """Evaluate (self ★ g)(x) via Gauss-Legendre quadrature."""
1148 x = np.atleast_1d(np.asarray(x, dtype=float))
1149 result = np.zeros(x.shape)
1150 for idx in range(x.size):
1151 xi = x[idx]
1152 t_lo = max(f_a, xi - g_d)
1153 t_hi = min(f_b, xi - g_c)
1154 if t_hi <= t_lo:
1155 continue
1156 # Break integration at breakpoints of f and shifted breakpoints
1157 # of g so the integrand is polynomial on each sub-interval.
1158 inner = [t_lo, t_hi]
1159 for bp in f_bps:
1160 if t_lo < bp < t_hi:
1161 inner.append(bp)
1162 for bp in g_bps:
1163 shifted = xi - bp
1164 if t_lo < shifted < t_hi:
1165 inner.append(shifted)
1166 inner = sorted(set(inner))
1168 total = 0.0
1169 for j in range(len(inner) - 1):
1170 a_int, b_int = inner[j], inner[j + 1]
1171 hw = (b_int - a_int) / 2.0
1172 mid = (a_int + b_int) / 2.0
1173 nodes = hw * quad_nodes + mid
1174 wts = hw * quad_weights
1175 total += np.dot(wts, self(nodes) * g(xi - nodes))
1176 result[idx] = total
1177 return result
1179 # Build a fun on each output sub-interval. Outermost pieces are
1180 # wrapped as CompactFun when the corresponding logical edge is ±inf;
1181 # interior pieces are always finite Bndfuns.
1182 n_pieces = len(out_breaks) - 1
1183 funs_list = []
1184 for i in range(n_pieces):
1185 a_storage = float(out_breaks[i])
1186 b_storage = float(out_breaks[i + 1])
1187 interval = Interval(a_storage, b_storage)
1188 bnd = Bndfun.initfun_adaptive(conv_eval, interval)
1189 is_first = i == 0
1190 is_last = i == n_pieces - 1
1191 wrap_left = is_first and left_inf
1192 wrap_right = is_last and right_inf
1193 if wrap_left or wrap_right:
1194 a_logical = -np.inf if wrap_left else a_storage
1195 b_logical = np.inf if wrap_right else b_storage
1196 funs_list.append(CompactFun(bnd.onefun, interval, logical_interval=(a_logical, b_logical)))
1197 else:
1198 funs_list.append(bnd)
1200 return self.__class__(funs_list)
1202 def sum(self) -> Any:
1203 """Compute the definite integral of the Chebfun over its domain.
1205 This method calculates the definite integral of the Chebfun over its
1206 entire domain of definition by summing the definite integrals of each
1207 piece.
1209 Returns:
1210 float or complex: The definite integral of the Chebfun over its domain.
1212 Examples:
1213 >>> import numpy as np
1214 >>> f = Chebfun.initfun_adaptive(lambda x: x**2, [-1, 1])
1215 >>> bool(abs(f.sum() - 2.0/3.0) < 1e-10)
1216 True
1217 >>> g = Chebfun.initconst(1.0, [-1, 1])
1218 >>> bool(abs(g.sum() - 2.0) < 1e-10)
1219 True
1220 """
1221 return np.sum([fun.sum() for fun in self])
1223 def dot(self, f: Any) -> Any:
1224 """Compute the dot product of this Chebfun with another function.
1226 This method calculates the inner product (dot product) of this Chebfun
1227 with another function f by multiplying them pointwise and then integrating
1228 the result over the domain.
1230 Args:
1231 f (Chebfun or scalar): The function or scalar to compute the dot product with.
1232 If not a Chebfun, it will be converted to one.
1234 Returns:
1235 float or complex: The dot product of this Chebfun with f.
1236 """
1237 return (self * f).sum()
1239 def norm(self, p: Any = 2) -> Any:
1240 """Compute the Lp norm of the Chebfun over its domain.
1242 This method calculates the Lp norm of the Chebfun. The L2 norm is the
1243 default and is computed as sqrt(integral(|f|^2)). For p=inf, returns
1244 the maximum absolute value by checking critical points (extrema).
1246 Args:
1247 p (int or float): The norm type. Supported values are 1, 2, positive
1248 integers/floats, or np.inf. Defaults to 2 (L2 norm).
1250 Returns:
1251 float: The Lp norm of the Chebfun.
1253 Examples:
1254 >>> from chebpy import chebfun
1255 >>> import numpy as np
1256 >>> f = chebfun(lambda x: x**2, [-1, 1])
1257 >>> np.allclose(f.norm(), 0.6324555320336759) # L2 norm
1258 True
1259 >>> np.allclose(f.norm(np.inf), 1.0) # Maximum absolute value
1260 True
1261 """
1262 if p == 2:
1263 # L2 norm: sqrt(integral(|f|^2))
1264 return np.sqrt(self.dot(self))
1265 elif p == np.inf:
1266 # L-infinity norm: max|f(x)|
1267 df = self.diff()
1268 critical_pts = df.roots()
1269 # Add endpoints
1270 endpoints = np.array([self.domain[0], self.domain[-1]])
1271 # Combine all test points
1272 test_pts = np.concatenate([critical_pts, endpoints])
1273 # Evaluate and find max
1274 vals = np.abs(self(test_pts))
1275 return np.max(vals)
1276 elif p == 1:
1277 # L1 norm: integral(|f|)
1278 return self.absolute().sum()
1279 elif p > 0:
1280 # General Lp norm: (integral(|f|^p))^(1/p)
1281 f_abs = self.absolute()
1282 f_pow_p = f_abs**p
1283 integral = f_pow_p.sum()
1284 return integral ** (1.0 / p)
1285 else:
1286 raise ValueError(p)
1288 # ----------
1289 # utilities
1290 # ----------
1291 @self_empty()
1292 def absolute(self) -> Chebfun:
1293 """Absolute value of a Chebfun."""
1294 newdom = self.domain.merge(self.roots())
1295 funs = [x.absolute() for x in self._break(newdom)]
1296 return self.__class__(funs)
1298 abs = absolute
1300 @self_empty()
1301 def sign(self) -> Chebfun:
1302 """Sign function of a Chebfun.
1304 Computes the piecewise sign of a Chebfun by finding its roots
1305 and splitting the domain at those points, then creating constant
1306 pieces with the appropriate sign values.
1308 Returns:
1309 Chebfun: A new Chebfun representing sign(f(x)).
1310 """
1311 roots = self.roots()
1312 newdom = self.domain.merge(roots)
1313 funs = []
1314 for fun in self._break(newdom):
1315 mid = fun.support[0] + 0.5 * (fun.support[-1] - fun.support[0])
1316 s = float(np.sign(float(self(mid))))
1317 funs.append(Bndfun.initconst(s, fun.interval))
1318 result = self.__class__(funs)
1319 # Set breakdata: at roots sign is 0, elsewhere use sign of function
1320 htol = max(1e2 * self.hscale * prefs.eps, prefs.eps)
1321 for bp in result.breakpoints:
1322 if roots.size > 0 and np.any(np.abs(bp - roots) <= htol):
1323 result.breakdata[bp] = 0.0
1324 else:
1325 result.breakdata[bp] = float(np.sign(float(self(bp))))
1326 return result
1328 @self_empty()
1329 def ceil(self) -> Chebfun:
1330 """Ceiling function of a Chebfun.
1332 Computes the piecewise ceiling of a Chebfun by finding where
1333 the function crosses integer values and splitting the domain
1334 at those points, then creating constant pieces with the
1335 appropriate ceiling values.
1337 Returns:
1338 Chebfun: A new Chebfun representing ceil(f(x)).
1339 """
1340 crossings = self._integer_crossings()
1341 newdom = self.domain.merge(crossings)
1342 funs = []
1343 for fun in self._break(newdom):
1344 mid = fun.support[0] + 0.5 * (fun.support[-1] - fun.support[0])
1345 c = float(np.ceil(float(self(mid))))
1346 funs.append(Bndfun.initconst(c, fun.interval))
1347 result = self.__class__(funs)
1348 for bp in result.breakpoints:
1349 result.breakdata[bp] = float(np.ceil(float(self(bp))))
1350 return result
1352 @self_empty()
1353 def floor(self) -> Chebfun:
1354 """Floor function of a Chebfun.
1356 Computes the piecewise floor of a Chebfun by finding where
1357 the function crosses integer values and splitting the domain
1358 at those points, then creating constant pieces with the
1359 appropriate floor values.
1361 Returns:
1362 Chebfun: A new Chebfun representing floor(f(x)).
1363 """
1364 crossings = self._integer_crossings()
1365 newdom = self.domain.merge(crossings)
1366 funs = []
1367 for fun in self._break(newdom):
1368 mid = fun.support[0] + 0.5 * (fun.support[-1] - fun.support[0])
1369 c = float(np.floor(float(self(mid))))
1370 funs.append(Bndfun.initconst(c, fun.interval))
1371 result = self.__class__(funs)
1372 for bp in result.breakpoints:
1373 result.breakdata[bp] = float(np.floor(float(self(bp))))
1374 return result
1376 def _integer_crossings(self) -> np.ndarray:
1377 """Find where this Chebfun crosses integer values.
1379 This helper method identifies all points in the domain where the
1380 Chebfun value equals an integer, by finding roots of (self - n)
1381 for each integer n in the range of the function.
1383 Returns:
1384 numpy.ndarray: Array of x-values where the function crosses integers.
1385 """
1386 all_values = np.concatenate([fun.values() for fun in self])
1387 lo = int(np.floor(np.min(all_values)))
1388 hi = int(np.ceil(np.max(all_values)))
1389 crossings = []
1390 for n in range(lo, hi + 1):
1391 shifted = self - n
1392 crossings.extend(shifted.roots().tolist())
1393 return np.array(crossings)
1395 @self_empty()
1396 @cast_arg_to_chebfun
1397 def maximum(self, other: Any) -> Any:
1398 """Pointwise maximum of self and another chebfun."""
1399 return self._maximum_minimum(other, operator.ge)
1401 @self_empty()
1402 @cast_arg_to_chebfun
1403 def minimum(self, other: Any) -> Any:
1404 """Pointwise mimimum of self and another chebfun."""
1405 return self._maximum_minimum(other, operator.lt)
1407 def _maximum_minimum(self, other: Chebfun, comparator: Callable[..., bool]) -> Any:
1408 """Method for computing the pointwise maximum/minimum of two Chebfuns.
1410 This internal method implements the algorithm for computing the pointwise
1411 maximum or minimum of two Chebfun objects, based on the provided comparator.
1412 It is used by the maximum() and minimum() methods.
1414 Args:
1415 other (Chebfun): Another Chebfun to compare with this one.
1416 comparator (callable): A function that compares two values and returns
1417 a boolean. For maximum, this is operator.ge (>=), and for minimum,
1418 this is operator.lt (<).
1420 Returns:
1421 Chebfun: A new Chebfun representing the pointwise maximum or minimum.
1422 """
1423 # Handle empty Chebfuns
1424 if self.isempty or other.isempty:
1425 return self.__class__.initempty()
1427 # Find the intersection of domains
1428 try:
1429 # Try to use union if supports match
1430 newdom = self.domain.union(other.domain)
1431 except SupportMismatch:
1432 # If supports don't match, find the intersection
1433 a_min, a_max = self.support
1434 b_min, b_max = other.support
1436 # Calculate intersection
1437 c_min = max(a_min, b_min)
1438 c_max = min(a_max, b_max)
1440 # If there's no intersection, return empty
1441 if c_min >= c_max:
1442 return self.__class__.initempty()
1444 # Restrict both functions to the intersection
1445 self_restricted = self.restrict([c_min, c_max])
1446 other_restricted = other.restrict([c_min, c_max])
1448 # Recursively call with the restricted functions
1449 return self_restricted._maximum_minimum(other_restricted, comparator)
1451 # Continue with the original algorithm
1452 roots = (self - other).roots()
1453 newdom = newdom.merge(roots)
1454 switch = newdom.support.merge(roots)
1456 # Handle the case where switch is empty
1457 if switch.size == 0: # pragma: no cover
1458 return self.__class__.initempty()
1460 keys = 0.5 * ((-1) ** np.arange(switch.size - 1) + 1)
1461 if switch.size > 0 and comparator(other(switch[0]), self(switch[0])):
1462 keys = 1 - keys
1463 funs = np.array([])
1464 for interval, use_self in zip(switch.intervals, keys, strict=False):
1465 subdom = newdom.restrict(interval)
1466 subfun = self.restrict(subdom) if use_self else other.restrict(subdom)
1467 funs = np.append(funs, subfun.funs)
1468 return self.__class__(funs)
1470 # ----------
1471 # plotting
1472 # ----------
1473 def plot(self, ax: Axes | None = None, **kwds: Any) -> Any:
1474 """Plot the Chebfun over its domain.
1476 This method plots the Chebfun over its domain using matplotlib.
1477 For complex-valued Chebfuns, it plots the real part against the imaginary part.
1479 For Chebfuns with ``±inf`` endpoints (containing :class:`CompactFun`
1480 pieces), each unbounded endpoint is replaced for plotting purposes
1481 with the corresponding ``plot_support`` endpoint of the outermost
1482 :class:`CompactFun` piece, so the decay-to-zero region is visible.
1484 Args:
1485 ax (matplotlib.axes.Axes, optional): The axes on which to plot. If None,
1486 a new axes will be created. Defaults to None.
1487 **kwds: Additional keyword arguments to pass to matplotlib's plot function.
1489 Returns:
1490 matplotlib.axes.Axes: The axes on which the plot was created.
1491 """
1492 a, b = float(self.support[0]), float(self.support[-1])
1493 if not np.isfinite(a):
1494 left = self.funs[0]
1495 a = float(left.plot_support[0]) if hasattr(left, "plot_support") else a
1496 if not np.isfinite(b):
1497 right = self.funs[-1]
1498 b = float(right.plot_support[1]) if hasattr(right, "plot_support") else b
1499 support = kwds.pop("support", (a, b))
1500 return plotfun(self, support, ax=ax, **kwds)
1502 def plotcoeffs(self, ax: Axes | None = None, **kwds: Any) -> Axes:
1503 """Plot the coefficients of the Chebfun on a semilogy scale.
1505 This method plots the absolute values of the coefficients for each piece
1506 of the Chebfun on a semilogy scale, which is useful for visualizing the
1507 decay of coefficients in the Chebyshev series.
1509 Args:
1510 ax (matplotlib.axes.Axes, optional): The axes on which to plot. If None,
1511 a new axes will be created. Defaults to None.
1512 **kwds: Additional keyword arguments to pass to matplotlib's semilogy function.
1514 Returns:
1515 matplotlib.axes.Axes: The axes on which the plot was created.
1516 """
1517 ax = ax or plt.gca()
1518 for fun in self:
1519 fun.plotcoeffs(ax=ax, **kwds)
1520 return ax
1523# ---------
1524# ufuncs
1525# ---------
1526def add_ufunc(op: Callable[..., Any]) -> None:
1527 """Add a NumPy universal function method to the Chebfun class.
1529 This function creates a method that applies a NumPy universal function (ufunc)
1530 to each piece of a Chebfun and returns a new Chebfun representing the result.
1532 Args:
1533 op (callable): The NumPy universal function to apply.
1535 Note:
1536 The created method will have the same name as the NumPy function
1537 and will take no arguments other than self.
1538 """
1540 @self_empty()
1541 def method(self: Chebfun) -> Chebfun:
1542 """Apply a NumPy universal function to this Chebfun.
1544 This method applies a NumPy universal function (ufunc) to each piece
1545 of this Chebfun and returns a new Chebfun representing the result.
1547 Args:
1548 self (Chebfun): The Chebfun object to which the function is applied.
1550 Returns:
1551 Chebfun: A new Chebfun representing op(f(x)).
1552 """
1553 return self.__class__([op(fun) for fun in self])
1555 name = op.__name__ # type: ignore[attr-defined]
1556 method.__name__ = name
1557 method.__doc__ = method.__doc__
1558 setattr(Chebfun, name, method)
1561ufuncs = (
1562 np.arccos,
1563 np.arccosh,
1564 np.arcsin,
1565 np.arcsinh,
1566 np.arctan,
1567 np.arctanh,
1568 np.cos,
1569 np.cosh,
1570 np.exp,
1571 np.exp2,
1572 np.expm1,
1573 np.log,
1574 np.log2,
1575 np.log10,
1576 np.log1p,
1577 np.sinh,
1578 np.sin,
1579 np.tan,
1580 np.tanh,
1581 np.sqrt,
1582)
1584for op in ufuncs:
1585 add_ufunc(op)