Coverage for src / chebpy / chebfun.py: 98%
347 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-20 05:55 +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"""
12import operator
14import matplotlib.pyplot as plt
15import numpy as np
17from .bndfun import Bndfun
18from .decorators import cache, cast_arg_to_chebfun, float_argument, self_empty
19from .exceptions import BadFunLengthArgument, SupportMismatch
20from .plotting import plotfun
21from .settings import _preferences as prefs
22from .utilities import Domain, check_funs, compute_breakdata, generate_funs
25class Chebfun:
26 """Main class for representing and manipulating functions in ChebPy.
28 The Chebfun class represents functions using piecewise polynomial approximations
29 on arbitrary intervals. It provides a comprehensive set of operations for working
30 with these function representations, including:
32 - Function evaluation at arbitrary points
33 - Algebraic operations (addition, multiplication, etc.)
34 - Calculus operations (differentiation, integration, etc.)
35 - Rootfinding
36 - Plotting
38 Chebfun objects can be created from callable functions, constant values, or
39 directly from function pieces. The class supports both adaptive and fixed-length
40 approximations, allowing for efficient representation of functions with varying
41 complexity across different intervals.
43 Attributes:
44 funs (numpy.ndarray): Array of function pieces that make up the Chebfun.
45 breakdata (OrderedDict): Mapping of breakpoints to function values.
46 transposed (bool): Flag indicating if the Chebfun is transposed.
47 """
49 def __init__(self, funs):
50 """Initialize a Chebfun object.
52 Args:
53 funs (list): List of function objects to be included in the Chebfun.
54 These will be checked and sorted using check_funs.
55 """
56 self.funs = check_funs(funs)
57 self.breakdata = compute_breakdata(self.funs)
58 self.transposed = False
60 @classmethod
61 def initempty(cls):
62 """Initialize an empty Chebfun.
64 Returns:
65 Chebfun: An empty Chebfun object with no functions.
67 Examples:
68 >>> f = Chebfun.initempty()
69 >>> f.isempty
70 True
71 """
72 return cls([])
74 @classmethod
75 def initidentity(cls, domain=None):
76 """Initialize a Chebfun representing the identity function f(x) = x.
78 Args:
79 domain (array-like, optional): Domain on which to define the identity function.
80 If None, uses the default domain from preferences.
82 Returns:
83 Chebfun: A Chebfun object representing the identity function on the specified domain.
85 Examples:
86 >>> import numpy as np
87 >>> x = Chebfun.initidentity([-1, 1])
88 >>> float(x(0.5))
89 0.5
90 >>> np.allclose(x([0, 0.5, 1]), [0, 0.5, 1])
91 True
92 """
93 return cls(generate_funs(domain, Bndfun.initidentity))
95 @classmethod
96 def initconst(cls, c, domain=None):
97 """Initialize a Chebfun representing a constant function f(x) = c.
99 Args:
100 c (float or complex): The constant value.
101 domain (array-like, optional): Domain on which to define the constant function.
102 If None, uses the default domain from preferences.
104 Returns:
105 Chebfun: A Chebfun object representing the constant function on the specified domain.
107 Examples:
108 >>> import numpy as np
109 >>> f = Chebfun.initconst(3.14, [-1, 1])
110 >>> float(f(0))
111 3.14
112 >>> float(f(0.5))
113 3.14
114 >>> np.allclose(f([0, 0.5, 1]), [3.14, 3.14, 3.14])
115 True
116 """
117 return cls(generate_funs(domain, Bndfun.initconst, {"c": c}))
119 @classmethod
120 def initfun_adaptive(cls, f, domain=None):
121 """Initialize a Chebfun by adaptively sampling a function.
123 This method determines the appropriate number of points needed to represent
124 the function to the specified tolerance using an adaptive algorithm.
126 Args:
127 f (callable): The function to be approximated.
128 domain (array-like, optional): Domain on which to define the function.
129 If None, uses the default domain from preferences.
131 Returns:
132 Chebfun: A Chebfun object representing the function on the specified domain.
134 Examples:
135 >>> import numpy as np
136 >>> f = Chebfun.initfun_adaptive(lambda x: np.sin(x), [-np.pi, np.pi])
137 >>> bool(abs(f(0)) < 1e-10)
138 True
139 >>> bool(abs(f(np.pi/2) - 1) < 1e-10)
140 True
141 """
142 return cls(generate_funs(domain, Bndfun.initfun_adaptive, {"f": f}))
144 @classmethod
145 def initfun_fixedlen(cls, f, n, domain=None):
146 """Initialize a Chebfun with a fixed number of points.
148 This method uses a specified number of points to represent the function,
149 rather than determining the number adaptively.
151 Args:
152 f (callable): The function to be approximated.
153 n (int or array-like): Number of points to use. If a single value, uses the same
154 number for each interval. If an array, must have one fewer elements than
155 the size of the domain.
156 domain (array-like, optional): Domain on which to define the function.
157 If None, uses the default domain from preferences.
159 Returns:
160 Chebfun: A Chebfun object representing the function on the specified domain.
162 Raises:
163 BadFunLengthArgument: If n is an array and its size doesn't match domain.size - 1.
164 """
165 nn = np.array(n)
166 if nn.size < 2:
167 funs = generate_funs(domain, Bndfun.initfun_fixedlen, {"f": f, "n": n})
168 else:
169 domain = Domain(domain if domain is not None else prefs.domain)
170 if not nn.size == domain.size - 1:
171 raise BadFunLengthArgument
172 funs = []
173 for interval, length in zip(domain.intervals, nn):
174 funs.append(Bndfun.initfun_fixedlen(f, interval, length))
175 return cls(funs)
177 @classmethod
178 def initfun(cls, f, domain=None, n=None):
179 """Initialize a Chebfun from a function.
181 This is a general-purpose constructor that delegates to either initfun_adaptive
182 or initfun_fixedlen based on whether n is provided.
184 Args:
185 f (callable): The function to be approximated.
186 domain (array-like, optional): Domain on which to define the function.
187 If None, uses the default domain from preferences.
188 n (int or array-like, optional): Number of points to use. If None, determines
189 the number adaptively. If provided, uses a fixed number of points.
191 Returns:
192 Chebfun: A Chebfun object representing the function on the specified domain.
193 """
194 if n is None:
195 return cls.initfun_adaptive(f, domain)
196 else:
197 return cls.initfun_fixedlen(f, n, domain)
199 # --------------------
200 # operator overloads
201 # --------------------
202 def __add__(self, f):
203 """Add a Chebfun with another Chebfun or a scalar.
205 Args:
206 f (Chebfun or scalar): The object to add to this Chebfun.
208 Returns:
209 Chebfun: A new Chebfun representing the sum.
210 """
211 return self._apply_binop(f, operator.add)
213 @self_empty(np.array([]))
214 @float_argument
215 def __call__(self, x):
216 """Evaluate the Chebfun at points x.
218 This method evaluates the Chebfun at the specified points. It handles interior
219 points, breakpoints, and points outside the domain appropriately.
221 Args:
222 x (float or array-like): Points at which to evaluate the Chebfun.
224 Returns:
225 float or numpy.ndarray: The value(s) of the Chebfun at the specified point(s).
226 Returns a scalar if x is a scalar, otherwise an array of the same size as x.
227 """
228 # initialise output
229 dtype = complex if self.iscomplex else float
230 out = np.full(x.size, np.nan, dtype=dtype)
232 # evaluate a fun when x is an interior point
233 for fun in self:
234 idx = fun.interval.isinterior(x)
235 out[idx] = fun(x[idx])
237 # evaluate the breakpoint data for x at a breakpoint
238 breakpoints = self.breakpoints
239 for break_point in breakpoints:
240 out[x == break_point] = self.breakdata[break_point]
242 # first and last funs used to evaluate outside of the chebfun domain
243 lpts, rpts = x < breakpoints[0], x > breakpoints[-1]
244 out[lpts] = self.funs[0](x[lpts])
245 out[rpts] = self.funs[-1](x[rpts])
246 return out
248 def __iter__(self):
249 """Return an iterator over the functions in this Chebfun.
251 Returns:
252 iterator: An iterator over the functions (funs) in this Chebfun.
253 """
254 return self.funs.__iter__()
256 def __len__(self):
257 """Return the total number of coefficients across all funs.
259 Returns:
260 int: The sum of sizes of all constituent funs.
261 """
262 return sum(f.size for f in self.funs)
264 def __eq__(self, other):
265 """Test for equality between two Chebfun objects.
267 Two Chebfun objects are considered equal if they have the same domain
268 and their function values are equal (within tolerance) at a set of test points.
270 Args:
271 other (object): The object to compare with this Chebfun.
273 Returns:
274 bool: True if the objects are equal, False otherwise.
275 """
276 if not isinstance(other, self.__class__):
277 return False
279 # Check if both are empty
280 if self.isempty and other.isempty:
281 return True
283 # Check if domains are equal
284 if self.domain != other.domain:
285 return False
287 # Check function values at test points
288 xx = np.linspace(self.support[0], self.support[1], 100)
289 tol = 1e2 * max(self.vscale, other.vscale) * prefs.eps
290 return np.all(np.abs(self(xx) - other(xx)) <= tol)
292 def __mul__(self, f):
293 """Multiply a Chebfun with another Chebfun or a scalar.
295 Args:
296 f (Chebfun or scalar): The object to multiply with this Chebfun.
298 Returns:
299 Chebfun: A new Chebfun representing the product.
300 """
301 return self._apply_binop(f, operator.mul)
303 def __neg__(self):
304 """Return the negative of this Chebfun.
306 Returns:
307 Chebfun: A new Chebfun representing -f(x).
308 """
309 return self.__class__(-self.funs)
311 def __pos__(self):
312 """Return the positive of this Chebfun (which is the Chebfun itself).
314 Returns:
315 Chebfun: This Chebfun object (unchanged).
316 """
317 return self
319 def __abs__(self):
320 """Return the absolute value of this Chebfun.
322 Returns:
323 Chebfun: A new Chebfun representing |f(x)|.
324 """
325 abs_funs = []
326 for fun in self.funs:
327 abs_funs.append(fun.absolute())
328 return self.__class__(abs_funs)
330 def __pow__(self, f):
331 """Raise this Chebfun to a power.
333 Args:
334 f (Chebfun or scalar): The exponent to which this Chebfun is raised.
336 Returns:
337 Chebfun: A new Chebfun representing self^f.
338 """
339 return self._apply_binop(f, operator.pow)
341 def __rtruediv__(self, c):
342 """Divide a scalar by this Chebfun.
344 This method is called when a scalar is divided by a Chebfun, i.e., c / self.
346 Args:
347 c (scalar): The scalar numerator.
349 Returns:
350 Chebfun: A new Chebfun representing c / self.
352 Note:
353 This is executed when truediv(f, self) fails, which is to say whenever c
354 is not a Chebfun. We proceed on the assumption f is a scalar.
355 """
357 def constfun(cheb, const):
358 return 0.0 * cheb + const
360 newfuns = [fun.initfun_adaptive(lambda x: constfun(x, c) / fun(x), fun.interval) for fun in self]
361 return self.__class__(newfuns)
363 @self_empty("Chebfun<empty>")
364 def __repr__(self):
365 """Return a string representation of the Chebfun.
367 This method returns a detailed string representation of the Chebfun,
368 including information about its domain, intervals, and endpoint values.
370 Returns:
371 str: A string representation of the Chebfun.
372 """
373 rowcol = "row" if self.transposed else "column"
374 numpcs = self.funs.size
375 plural = "" if numpcs == 1 else "s"
376 header = f"Chebfun {rowcol} ({numpcs} smooth piece{plural})\n"
377 domain_info = f"domain: {self.support}\n"
378 toprow = " interval length endpoint values\n"
379 tmplat = "[{:8.2g},{:8.2g}] {:6} {:8.2g} {:8.2g}\n"
380 rowdta = ""
381 for fun in self:
382 endpts = fun.support
383 xl, xr = endpts
384 fl, fr = fun(endpts)
385 row = tmplat.format(xl, xr, fun.size, fl, fr)
386 rowdta += row
387 btmrow = f"vertical scale = {self.vscale:3.2g}"
388 btmxtr = "" if numpcs == 1 else f" total length = {sum([f.size for f in self])}"
389 return header + domain_info + toprow + rowdta + btmrow + btmxtr
391 def __rsub__(self, f):
392 """Subtract this Chebfun from another object.
394 This method is called when another object is subtracted by this Chebfun,
395 i.e., f - self.
397 Args:
398 f (Chebfun or scalar): The object from which to subtract this Chebfun.
400 Returns:
401 Chebfun: A new Chebfun representing f - self.
402 """
403 return -(self - f)
405 @cast_arg_to_chebfun
406 def __rpow__(self, f):
407 """Raise another object to the power of this Chebfun.
409 This method is called when another object is raised to the power of this Chebfun,
410 i.e., f ** self.
412 Args:
413 f (Chebfun or scalar): The base to be raised to the power of this Chebfun.
415 Returns:
416 Chebfun: A new Chebfun representing f ** self.
417 """
418 return f**self
420 def __truediv__(self, f):
421 """Divide this Chebfun by another object.
423 Args:
424 f (Chebfun or scalar): The divisor.
426 Returns:
427 Chebfun: A new Chebfun representing self / f.
428 """
429 return self._apply_binop(f, operator.truediv)
431 __rmul__ = __mul__
432 __div__ = __truediv__
433 __rdiv__ = __rtruediv__
434 __radd__ = __add__
436 def __str__(self):
437 """Return a concise string representation of the Chebfun.
439 This method returns a brief string representation of the Chebfun,
440 showing its orientation, number of pieces, total size, and domain.
442 Returns:
443 str: A concise string representation of the Chebfun.
444 """
445 rowcol = "row" if self.transposed else "col"
446 domain_str = f"domain {self.support}" if not self.isempty else "empty"
447 out = f"<Chebfun-{rowcol},{self.funs.size},{sum([f.size for f in self])}, {domain_str}>\n"
448 return out
450 def __sub__(self, f):
451 """Subtract another object from this Chebfun.
453 Args:
454 f (Chebfun or scalar): The object to subtract from this Chebfun.
456 Returns:
457 Chebfun: A new Chebfun representing self - f.
458 """
459 return self._apply_binop(f, operator.sub)
461 # ------------------
462 # internal helpers
463 # ------------------
464 @self_empty()
465 def _apply_binop(self, f, op):
466 """Apply a binary operation between this Chebfun and another object.
468 This is a funnel method used in the implementation of Chebfun binary
469 operators. The high-level idea is to first break each chebfun into a
470 series of pieces corresponding to the union of the domains of each
471 before applying the supplied binary operator and simplifying. In the
472 case of the second argument being a scalar we don't need to do the
473 simplify step, since at the Tech-level these operations are defined
474 such that there is no change in the number of coefficients.
476 Args:
477 f (Chebfun or scalar): The second operand of the binary operation.
478 op (callable): The binary operation to apply (e.g., operator.add).
480 Returns:
481 Chebfun: A new Chebfun resulting from applying the binary operation.
482 """
483 if hasattr(f, "isempty") and f.isempty:
484 return f
485 if np.isscalar(f):
486 chbfn1 = self
487 chbfn2 = f * np.ones(self.funs.size)
488 simplify = False
489 else:
490 newdom = self.domain.union(f.domain)
491 chbfn1 = self._break(newdom)
492 chbfn2 = f._break(newdom)
493 simplify = True
494 newfuns = []
495 for fun1, fun2 in zip(chbfn1, chbfn2):
496 newfun = op(fun1, fun2)
497 if simplify:
498 newfun = newfun.simplify()
499 newfuns.append(newfun)
500 return self.__class__(newfuns)
502 def _break(self, targetdomain):
503 """Resample this Chebfun to a new domain.
505 This method resamples the Chebfun to the supplied Domain object. It is
506 intended as a private method since one will typically need to have
507 called either Domain.union(f) or Domain.merge(f) prior to calling this method.
509 Args:
510 targetdomain (Domain): The domain to which this Chebfun should be resampled.
512 Returns:
513 Chebfun: A new Chebfun resampled to the target domain.
514 """
515 newfuns = []
516 subintervals = targetdomain.intervals
517 interval = next(subintervals) # next(..) for Python2/3 compatibility
518 for fun in self:
519 while interval in fun.interval:
520 newfun = fun.restrict(interval)
521 newfuns.append(newfun)
522 try:
523 interval = next(subintervals)
524 except StopIteration:
525 break
526 return self.__class__(newfuns)
528 # ------------
529 # properties
530 # ------------
531 @property
532 def breakpoints(self):
533 """Get the breakpoints of this Chebfun.
535 Breakpoints are the points where the Chebfun transitions from one piece to another.
537 Returns:
538 numpy.ndarray: Array of breakpoints.
539 """
540 return np.array([x for x in self.breakdata.keys()])
542 @property
543 @self_empty(np.array([]))
544 def domain(self):
545 """Get the domain of this Chebfun.
547 Returns:
548 Domain: A Domain object corresponding to this Chebfun.
549 """
550 return Domain.from_chebfun(self)
552 @domain.setter
553 def domain(self, new_domain):
554 """Set the domain of the Chebfun by restricting to the new domain.
556 Args:
557 new_domain (array-like): The new domain to which this Chebfun should be restricted.
558 """
559 self.restrict_(new_domain)
561 @property
562 @self_empty(Domain([]))
563 def support(self):
564 """Get the support interval of this Chebfun.
566 The support is the interval between the first and last breakpoints.
568 Returns:
569 numpy.ndarray: Array containing the first and last breakpoints.
570 """
571 return self.domain.support
573 @property
574 @self_empty(0.0)
575 def hscale(self):
576 """Get the horizontal scale of this Chebfun.
578 The horizontal scale is the maximum absolute value of the support interval.
580 Returns:
581 float: The horizontal scale.
582 """
583 return float(np.abs(self.support).max())
585 @property
586 @self_empty(False)
587 def iscomplex(self):
588 """Check if this Chebfun has complex values.
590 Returns:
591 bool: True if any of the functions in this Chebfun have complex values,
592 False otherwise.
593 """
594 return any(fun.iscomplex for fun in self)
596 @property
597 @self_empty(False)
598 def isconst(self):
599 """Check if this Chebfun represents a constant function.
601 A Chebfun is constant if all of its pieces are constant with the same value.
603 Returns:
604 bool: True if this Chebfun represents a constant function, False otherwise.
606 Note:
607 TODO: find an abstract way of referencing funs[0].coeffs[0]
608 """
609 c = self.funs[0].coeffs[0]
610 return all(fun.isconst and fun.coeffs[0] == c for fun in self)
612 @property
613 def isempty(self):
614 """Check if this Chebfun is empty.
616 An empty Chebfun contains no functions.
618 Returns:
619 bool: True if this Chebfun is empty, False otherwise.
620 """
621 return self.funs.size == 0
623 @property
624 @self_empty(0.0)
625 def vscale(self):
626 """Get the vertical scale of this Chebfun.
628 The vertical scale is the maximum of the vertical scales of all pieces.
630 Returns:
631 float: The vertical scale.
632 """
633 return np.max([fun.vscale for fun in self])
635 @property
636 @self_empty()
637 def x(self):
638 """Get the identity function on the support of this Chebfun.
640 This property returns a new Chebfun representing the identity function f(x) = x
641 defined on the same support as this Chebfun.
643 Returns:
644 Chebfun: A Chebfun representing the identity function on the support of this Chebfun.
645 """
646 return self.__class__.initidentity(self.support)
648 # -----------
649 # utilities
650 # ----------
652 def imag(self):
653 """Get the imaginary part of this Chebfun.
655 Returns:
656 Chebfun: A new Chebfun representing the imaginary part of this Chebfun.
657 If this Chebfun is real-valued, returns a zero Chebfun.
658 """
659 if self.iscomplex:
660 return self.__class__([fun.imag() for fun in self])
661 else:
662 return self.initconst(0, domain=self.domain)
664 def real(self):
665 """Get the real part of this Chebfun.
667 Returns:
668 Chebfun: A new Chebfun representing the real part of this Chebfun.
669 If this Chebfun is already real-valued, returns this Chebfun.
670 """
671 if self.iscomplex:
672 return self.__class__([fun.real() for fun in self])
673 else:
674 return self
676 def copy(self):
677 """Create a deep copy of this Chebfun.
679 Returns:
680 Chebfun: A new Chebfun that is a deep copy of this Chebfun.
681 """
682 return self.__class__([fun.copy() for fun in self])
684 @self_empty()
685 def _restrict(self, subinterval):
686 """Restrict a Chebfun to a subinterval, without simplifying.
688 This is an internal method that restricts the Chebfun to a subinterval
689 without performing simplification.
691 Args:
692 subinterval (array-like): The subinterval to which this Chebfun should be restricted.
694 Returns:
695 Chebfun: A new Chebfun restricted to the specified subinterval, without simplification.
696 """
697 newdom = self.domain.restrict(Domain(subinterval))
698 return self._break(newdom)
700 def restrict(self, subinterval):
701 """Restrict a Chebfun to a subinterval.
703 This method creates a new Chebfun that is restricted to the specified subinterval
704 and simplifies the result.
706 Args:
707 subinterval (array-like): The subinterval to which this Chebfun should be restricted.
709 Returns:
710 Chebfun: A new Chebfun restricted to the specified subinterval.
711 """
712 return self._restrict(subinterval).simplify()
714 @self_empty()
715 def restrict_(self, subinterval):
716 """Restrict a Chebfun to a subinterval, modifying the object in place.
718 This method modifies the current Chebfun by restricting it to the specified
719 subinterval and simplifying the result.
721 Args:
722 subinterval (array-like): The subinterval to which this Chebfun should be restricted.
724 Returns:
725 Chebfun: The modified Chebfun (self).
726 """
727 restricted = self._restrict(subinterval).simplify()
728 self.funs = restricted.funs
729 self.breakdata = compute_breakdata(self.funs)
730 return self
732 @cache
733 @self_empty(np.array([]))
734 def roots(self, merge=None):
735 """Compute the roots of a Chebfun.
737 This method finds the values x for which f(x) = 0, by computing the roots
738 of each piece of the Chebfun and combining them.
740 Args:
741 merge (bool, optional): Whether to merge roots at breakpoints. If None,
742 uses the value from preferences. Defaults to None.
744 Returns:
745 numpy.ndarray: Array of roots sorted in ascending order.
747 Examples:
748 >>> import numpy as np
749 >>> f = Chebfun.initfun_adaptive(lambda x: x**2 - 1, [-2, 2])
750 >>> roots = f.roots()
751 >>> len(roots)
752 2
753 >>> np.allclose(sorted(roots), [-1, 1])
754 True
755 """
756 merge = merge if merge is not None else prefs.mergeroots
757 allrts = []
758 prvrts = np.array([])
759 htol = 1e2 * self.hscale * prefs.eps
760 for fun in self:
761 rts = fun.roots()
762 # ignore first root if equal to the last root of previous fun
763 # TODO: there could be multiple roots at breakpoints
764 if prvrts.size > 0 and rts.size > 0:
765 if merge and abs(prvrts[-1] - rts[0]) <= htol:
766 rts = rts[1:]
767 allrts.append(rts)
768 prvrts = rts
769 return np.concatenate([x for x in allrts])
771 @self_empty()
772 def simplify(self):
773 """Simplify each fun in the chebfun."""
774 return self.__class__([fun.simplify() for fun in self])
776 def translate(self, c):
777 """Translate a chebfun by c, i.e., return f(x-c)."""
778 return self.__class__([x.translate(c) for x in self])
780 # ----------
781 # calculus
782 # ----------
783 def cumsum(self):
784 """Compute the indefinite integral (antiderivative) of the Chebfun.
786 This method computes the indefinite integral of the Chebfun, with the
787 constant of integration chosen so that the indefinite integral evaluates
788 to 0 at the left endpoint of the domain. For piecewise functions, constants
789 are added to ensure continuity across the pieces.
791 Returns:
792 Chebfun: A new Chebfun representing the indefinite integral of this Chebfun.
794 Examples:
795 >>> import numpy as np
796 >>> f = Chebfun.initconst(1.0, [-1, 1])
797 >>> F = f.cumsum()
798 >>> bool(abs(F(-1)) < 1e-10)
799 True
800 >>> bool(abs(F(1) - 2.0) < 1e-10)
801 True
802 """
803 newfuns = []
804 prevfun = None
805 for fun in self:
806 integral = fun.cumsum()
807 if prevfun:
808 # enforce continuity by adding the function value
809 # at the right endpoint of the previous fun
810 _, fb = prevfun.endvalues
811 integral = integral + fb
812 newfuns.append(integral)
813 prevfun = integral
814 return self.__class__(newfuns)
816 def diff(self, n=1):
817 """Compute the derivative of the Chebfun.
819 This method calculates the nth derivative of the Chebfun with respect to x.
820 It creates a new Chebfun where each piece is the derivative of the
821 corresponding piece in the original Chebfun.
823 Args:
824 n: Order of differentiation (default: 1). Must be non-negative integer.
826 Returns:
827 Chebfun: A new Chebfun representing the nth derivative of this Chebfun.
829 Examples:
830 >>> from chebpy import chebfun
831 >>> f = chebfun(lambda x: x**3)
832 >>> df1 = f.diff() # first derivative: 3*x**2
833 >>> df2 = f.diff(2) # second derivative: 6*x
834 >>> df3 = f.diff(3) # third derivative: 6
835 >>> bool(abs(df1(0.5) - 0.75) < 1e-10)
836 True
837 >>> bool(abs(df2(0.5) - 3.0) < 1e-10)
838 True
839 >>> bool(abs(df3(0.5) - 6.0) < 1e-10)
840 True
841 """
842 if not isinstance(n, int):
843 raise TypeError("Derivative order must be an integer")
844 if n == 0:
845 return self
846 if n < 0:
847 raise ValueError("Derivative order must be non-negative")
849 result = self
850 for _ in range(n):
851 dfuns = np.array([fun.diff() for fun in result])
852 result = self.__class__(dfuns)
853 return result
855 def sum(self):
856 """Compute the definite integral of the Chebfun over its domain.
858 This method calculates the definite integral of the Chebfun over its
859 entire domain of definition by summing the definite integrals of each
860 piece.
862 Returns:
863 float or complex: The definite integral of the Chebfun over its domain.
865 Examples:
866 >>> import numpy as np
867 >>> f = Chebfun.initfun_adaptive(lambda x: x**2, [-1, 1])
868 >>> bool(abs(f.sum() - 2.0/3.0) < 1e-10)
869 True
870 >>> g = Chebfun.initconst(1.0, [-1, 1])
871 >>> bool(abs(g.sum() - 2.0) < 1e-10)
872 True
873 """
874 return np.sum([fun.sum() for fun in self])
876 def dot(self, f):
877 """Compute the dot product of this Chebfun with another function.
879 This method calculates the inner product (dot product) of this Chebfun
880 with another function f by multiplying them pointwise and then integrating
881 the result over the domain.
883 Args:
884 f (Chebfun or scalar): The function or scalar to compute the dot product with.
885 If not a Chebfun, it will be converted to one.
887 Returns:
888 float or complex: The dot product of this Chebfun with f.
889 """
890 return (self * f).sum()
892 def norm(self, p=2):
893 """Compute the Lp norm of the Chebfun over its domain.
895 This method calculates the Lp norm of the Chebfun. The L2 norm is the
896 default and is computed as sqrt(integral(|f|^2)). For p=inf, returns
897 the maximum absolute value by checking critical points (extrema).
899 Args:
900 p (int or float): The norm type. Supported values are 1, 2, positive
901 integers/floats, or np.inf. Defaults to 2 (L2 norm).
903 Returns:
904 float: The Lp norm of the Chebfun.
906 Examples:
907 >>> from chebpy import chebfun
908 >>> import numpy as np
909 >>> f = chebfun(lambda x: x**2, [-1, 1])
910 >>> np.allclose(f.norm(), 0.6324555320336759) # L2 norm
911 True
912 >>> np.allclose(f.norm(np.inf), 1.0) # Maximum absolute value
913 True
914 """
915 if p == 2:
916 # L2 norm: sqrt(integral(|f|^2))
917 return np.sqrt(self.dot(self))
918 elif p == np.inf:
919 # L-infinity norm: max|f(x)|
920 df = self.diff()
921 critical_pts = df.roots()
922 # Add endpoints
923 endpoints = np.array([self.domain[0], self.domain[-1]])
924 # Combine all test points
925 test_pts = np.concatenate([critical_pts, endpoints])
926 # Evaluate and find max
927 vals = np.abs(self(test_pts))
928 return np.max(vals)
929 elif p == 1:
930 # L1 norm: integral(|f|)
931 return self.absolute().sum()
932 elif p > 0:
933 # General Lp norm: (integral(|f|^p))^(1/p)
934 f_abs = self.absolute()
935 f_pow_p = f_abs**p
936 integral = f_pow_p.sum()
937 return integral ** (1.0 / p)
938 else:
939 raise ValueError(f"norm(p={p}): p must be positive or np.inf")
941 # ----------
942 # utilities
943 # ----------
944 @self_empty()
945 def absolute(self):
946 """Absolute value of a Chebfun."""
947 newdom = self.domain.merge(self.roots())
948 funs = [x.absolute() for x in self._break(newdom)]
949 return self.__class__(funs)
951 abs = absolute
953 @self_empty()
954 @cast_arg_to_chebfun
955 def maximum(self, other):
956 """Pointwise maximum of self and another chebfun."""
957 return self._maximum_minimum(other, operator.ge)
959 @self_empty()
960 @cast_arg_to_chebfun
961 def minimum(self, other):
962 """Pointwise mimimum of self and another chebfun."""
963 return self._maximum_minimum(other, operator.lt)
965 def _maximum_minimum(self, other, comparator):
966 """Method for computing the pointwise maximum/minimum of two Chebfuns.
968 This internal method implements the algorithm for computing the pointwise
969 maximum or minimum of two Chebfun objects, based on the provided comparator.
970 It is used by the maximum() and minimum() methods.
972 Args:
973 other (Chebfun): Another Chebfun to compare with this one.
974 comparator (callable): A function that compares two values and returns
975 a boolean. For maximum, this is operator.ge (>=), and for minimum,
976 this is operator.lt (<).
978 Returns:
979 Chebfun: A new Chebfun representing the pointwise maximum or minimum.
980 """
981 # Handle empty Chebfuns
982 if self.isempty or other.isempty:
983 return self.__class__.initempty()
985 # Find the intersection of domains
986 try:
987 # Try to use union if supports match
988 newdom = self.domain.union(other.domain)
989 except SupportMismatch:
990 # If supports don't match, find the intersection
991 a_min, a_max = self.support
992 b_min, b_max = other.support
994 # Calculate intersection
995 c_min = max(a_min, b_min)
996 c_max = min(a_max, b_max)
998 # If there's no intersection, return empty
999 if c_min >= c_max:
1000 return self.__class__.initempty()
1002 # Restrict both functions to the intersection
1003 self_restricted = self.restrict([c_min, c_max])
1004 other_restricted = other.restrict([c_min, c_max])
1006 # Recursively call with the restricted functions
1007 return self_restricted._maximum_minimum(other_restricted, comparator)
1009 # Continue with the original algorithm
1010 roots = (self - other).roots()
1011 newdom = newdom.merge(roots)
1012 switch = newdom.support.merge(roots)
1014 # Handle the case where switch is empty
1015 if switch.size == 0: # pragma: no cover
1016 return self.__class__.initempty()
1018 keys = 0.5 * ((-1) ** np.arange(switch.size - 1) + 1)
1019 if switch.size > 0 and comparator(other(switch[0]), self(switch[0])):
1020 keys = 1 - keys
1021 funs = np.array([])
1022 for interval, use_self in zip(switch.intervals, keys):
1023 subdom = newdom.restrict(interval)
1024 if use_self:
1025 subfun = self.restrict(subdom)
1026 else:
1027 subfun = other.restrict(subdom)
1028 funs = np.append(funs, subfun.funs)
1029 return self.__class__(funs)
1031 # ----------
1032 # plotting
1033 # ----------
1034 def plot(self, ax=None, **kwds):
1035 """Plot the Chebfun over its domain.
1037 This method plots the Chebfun over its domain using matplotlib.
1038 For complex-valued Chebfuns, it plots the real part against the imaginary part.
1040 Args:
1041 ax (matplotlib.axes.Axes, optional): The axes on which to plot. If None,
1042 a new axes will be created. Defaults to None.
1043 **kwds: Additional keyword arguments to pass to matplotlib's plot function.
1045 Returns:
1046 matplotlib.axes.Axes: The axes on which the plot was created.
1047 """
1048 return plotfun(self, self.support, ax=ax, **kwds)
1050 def plotcoeffs(self, ax=None, **kwds):
1051 """Plot the coefficients of the Chebfun on a semilogy scale.
1053 This method plots the absolute values of the coefficients for each piece
1054 of the Chebfun on a semilogy scale, which is useful for visualizing the
1055 decay of coefficients in the Chebyshev series.
1057 Args:
1058 ax (matplotlib.axes.Axes, optional): The axes on which to plot. If None,
1059 a new axes will be created. Defaults to None.
1060 **kwds: Additional keyword arguments to pass to matplotlib's semilogy function.
1062 Returns:
1063 matplotlib.axes.Axes: The axes on which the plot was created.
1064 """
1065 ax = ax or plt.gca()
1066 for fun in self:
1067 fun.plotcoeffs(ax=ax, **kwds)
1068 return ax
1071# ---------
1072# ufuncs
1073# ---------
1074def add_ufunc(op):
1075 """Add a NumPy universal function method to the Chebfun class.
1077 This function creates a method that applies a NumPy universal function (ufunc)
1078 to each piece of a Chebfun and returns a new Chebfun representing the result.
1080 Args:
1081 op (callable): The NumPy universal function to apply.
1083 Note:
1084 The created method will have the same name as the NumPy function
1085 and will take no arguments other than self.
1086 """
1088 @self_empty()
1089 def method(self):
1090 """Apply a NumPy universal function to this Chebfun.
1092 This method applies a NumPy universal function (ufunc) to each piece
1093 of this Chebfun and returns a new Chebfun representing the result.
1095 Args:
1096 self (Chebfun): The Chebfun object to which the function is applied.
1098 Returns:
1099 Chebfun: A new Chebfun representing op(f(x)).
1100 """
1101 return self.__class__([op(fun) for fun in self])
1103 name = op.__name__
1104 method.__name__ = name
1105 method.__doc__ = method.__doc__
1106 setattr(Chebfun, name, method)
1109ufuncs = (
1110 np.arccos,
1111 np.arccosh,
1112 np.arcsin,
1113 np.arcsinh,
1114 np.arctan,
1115 np.arctanh,
1116 np.cos,
1117 np.cosh,
1118 np.exp,
1119 np.exp2,
1120 np.expm1,
1121 np.log,
1122 np.log2,
1123 np.log10,
1124 np.log1p,
1125 np.sinh,
1126 np.sin,
1127 np.tan,
1128 np.tanh,
1129 np.sqrt,
1130)
1132for op in ufuncs:
1133 add_ufunc(op)