Coverage for /usr/lib/python3/dist-packages/sympy/series/fourier.py: 23%
277 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Fourier Series"""
3from sympy.core.numbers import (oo, pi)
4from sympy.core.symbol import Wild
5from sympy.core.expr import Expr
6from sympy.core.add import Add
7from sympy.core.containers import Tuple
8from sympy.core.singleton import S
9from sympy.core.symbol import Dummy, Symbol
10from sympy.core.sympify import sympify
11from sympy.functions.elementary.trigonometric import sin, cos, sinc
12from sympy.series.series_class import SeriesBase
13from sympy.series.sequences import SeqFormula
14from sympy.sets.sets import Interval
15from sympy.utilities.iterables import is_sequence
18def fourier_cos_seq(func, limits, n):
19 """Returns the cos sequence in a Fourier series"""
20 from sympy.integrals import integrate
21 x, L = limits[0], limits[2] - limits[1]
22 cos_term = cos(2*n*pi*x / L)
23 formula = 2 * cos_term * integrate(func * cos_term, limits) / L
24 a0 = formula.subs(n, S.Zero) / 2
25 return a0, SeqFormula(2 * cos_term * integrate(func * cos_term, limits)
26 / L, (n, 1, oo))
29def fourier_sin_seq(func, limits, n):
30 """Returns the sin sequence in a Fourier series"""
31 from sympy.integrals import integrate
32 x, L = limits[0], limits[2] - limits[1]
33 sin_term = sin(2*n*pi*x / L)
34 return SeqFormula(2 * sin_term * integrate(func * sin_term, limits)
35 / L, (n, 1, oo))
38def _process_limits(func, limits):
39 """
40 Limits should be of the form (x, start, stop).
41 x should be a symbol. Both start and stop should be bounded.
43 Explanation
44 ===========
46 * If x is not given, x is determined from func.
47 * If limits is None. Limit of the form (x, -pi, pi) is returned.
49 Examples
50 ========
52 >>> from sympy.series.fourier import _process_limits as pari
53 >>> from sympy.abc import x
54 >>> pari(x**2, (x, -2, 2))
55 (x, -2, 2)
56 >>> pari(x**2, (-2, 2))
57 (x, -2, 2)
58 >>> pari(x**2, None)
59 (x, -pi, pi)
60 """
61 def _find_x(func):
62 free = func.free_symbols
63 if len(free) == 1:
64 return free.pop()
65 elif not free:
66 return Dummy('k')
67 else:
68 raise ValueError(
69 " specify dummy variables for %s. If the function contains"
70 " more than one free symbol, a dummy variable should be"
71 " supplied explicitly e.g. FourierSeries(m*n**2, (n, -pi, pi))"
72 % func)
74 x, start, stop = None, None, None
75 if limits is None:
76 x, start, stop = _find_x(func), -pi, pi
77 if is_sequence(limits, Tuple):
78 if len(limits) == 3:
79 x, start, stop = limits
80 elif len(limits) == 2:
81 x = _find_x(func)
82 start, stop = limits
84 if not isinstance(x, Symbol) or start is None or stop is None:
85 raise ValueError('Invalid limits given: %s' % str(limits))
87 unbounded = [S.NegativeInfinity, S.Infinity]
88 if start in unbounded or stop in unbounded:
89 raise ValueError("Both the start and end value should be bounded")
91 return sympify((x, start, stop))
94def finite_check(f, x, L):
96 def check_fx(exprs, x):
97 return x not in exprs.free_symbols
99 def check_sincos(_expr, x, L):
100 if isinstance(_expr, (sin, cos)):
101 sincos_args = _expr.args[0]
103 if sincos_args.match(a*(pi/L)*x + b) is not None:
104 return True
105 else:
106 return False
108 from sympy.simplify.fu import TR2, TR1, sincos_to_sum
109 _expr = sincos_to_sum(TR2(TR1(f)))
110 add_coeff = _expr.as_coeff_add()
112 a = Wild('a', properties=[lambda k: k.is_Integer, lambda k: k != S.Zero, ])
113 b = Wild('b', properties=[lambda k: x not in k.free_symbols, ])
115 for s in add_coeff[1]:
116 mul_coeffs = s.as_coeff_mul()[1]
117 for t in mul_coeffs:
118 if not (check_fx(t, x) or check_sincos(t, x, L)):
119 return False, f
121 return True, _expr
124class FourierSeries(SeriesBase):
125 r"""Represents Fourier sine/cosine series.
127 Explanation
128 ===========
130 This class only represents a fourier series.
131 No computation is performed.
133 For how to compute Fourier series, see the :func:`fourier_series`
134 docstring.
136 See Also
137 ========
139 sympy.series.fourier.fourier_series
140 """
141 def __new__(cls, *args):
142 args = map(sympify, args)
143 return Expr.__new__(cls, *args)
145 @property
146 def function(self):
147 return self.args[0]
149 @property
150 def x(self):
151 return self.args[1][0]
153 @property
154 def period(self):
155 return (self.args[1][1], self.args[1][2])
157 @property
158 def a0(self):
159 return self.args[2][0]
161 @property
162 def an(self):
163 return self.args[2][1]
165 @property
166 def bn(self):
167 return self.args[2][2]
169 @property
170 def interval(self):
171 return Interval(0, oo)
173 @property
174 def start(self):
175 return self.interval.inf
177 @property
178 def stop(self):
179 return self.interval.sup
181 @property
182 def length(self):
183 return oo
185 @property
186 def L(self):
187 return abs(self.period[1] - self.period[0]) / 2
189 def _eval_subs(self, old, new):
190 x = self.x
191 if old.has(x):
192 return self
194 def truncate(self, n=3):
195 """
196 Return the first n nonzero terms of the series.
198 If ``n`` is None return an iterator.
200 Parameters
201 ==========
203 n : int or None
204 Amount of non-zero terms in approximation or None.
206 Returns
207 =======
209 Expr or iterator :
210 Approximation of function expanded into Fourier series.
212 Examples
213 ========
215 >>> from sympy import fourier_series, pi
216 >>> from sympy.abc import x
217 >>> s = fourier_series(x, (x, -pi, pi))
218 >>> s.truncate(4)
219 2*sin(x) - sin(2*x) + 2*sin(3*x)/3 - sin(4*x)/2
221 See Also
222 ========
224 sympy.series.fourier.FourierSeries.sigma_approximation
225 """
226 if n is None:
227 return iter(self)
229 terms = []
230 for t in self:
231 if len(terms) == n:
232 break
233 if t is not S.Zero:
234 terms.append(t)
236 return Add(*terms)
238 def sigma_approximation(self, n=3):
239 r"""
240 Return :math:`\sigma`-approximation of Fourier series with respect
241 to order n.
243 Explanation
244 ===========
246 Sigma approximation adjusts a Fourier summation to eliminate the Gibbs
247 phenomenon which would otherwise occur at discontinuities.
248 A sigma-approximated summation for a Fourier series of a T-periodical
249 function can be written as
251 .. math::
252 s(\theta) = \frac{1}{2} a_0 + \sum _{k=1}^{m-1}
253 \operatorname{sinc} \Bigl( \frac{k}{m} \Bigr) \cdot
254 \left[ a_k \cos \Bigl( \frac{2\pi k}{T} \theta \Bigr)
255 + b_k \sin \Bigl( \frac{2\pi k}{T} \theta \Bigr) \right],
257 where :math:`a_0, a_k, b_k, k=1,\ldots,{m-1}` are standard Fourier
258 series coefficients and
259 :math:`\operatorname{sinc} \Bigl( \frac{k}{m} \Bigr)` is a Lanczos
260 :math:`\sigma` factor (expressed in terms of normalized
261 :math:`\operatorname{sinc}` function).
263 Parameters
264 ==========
266 n : int
267 Highest order of the terms taken into account in approximation.
269 Returns
270 =======
272 Expr :
273 Sigma approximation of function expanded into Fourier series.
275 Examples
276 ========
278 >>> from sympy import fourier_series, pi
279 >>> from sympy.abc import x
280 >>> s = fourier_series(x, (x, -pi, pi))
281 >>> s.sigma_approximation(4)
282 2*sin(x)*sinc(pi/4) - 2*sin(2*x)/pi + 2*sin(3*x)*sinc(3*pi/4)/3
284 See Also
285 ========
287 sympy.series.fourier.FourierSeries.truncate
289 Notes
290 =====
292 The behaviour of
293 :meth:`~sympy.series.fourier.FourierSeries.sigma_approximation`
294 is different from :meth:`~sympy.series.fourier.FourierSeries.truncate`
295 - it takes all nonzero terms of degree smaller than n, rather than
296 first n nonzero ones.
298 References
299 ==========
301 .. [1] https://en.wikipedia.org/wiki/Gibbs_phenomenon
302 .. [2] https://en.wikipedia.org/wiki/Sigma_approximation
303 """
304 terms = [sinc(pi * i / n) * t for i, t in enumerate(self[:n])
305 if t is not S.Zero]
306 return Add(*terms)
308 def shift(self, s):
309 """
310 Shift the function by a term independent of x.
312 Explanation
313 ===========
315 f(x) -> f(x) + s
317 This is fast, if Fourier series of f(x) is already
318 computed.
320 Examples
321 ========
323 >>> from sympy import fourier_series, pi
324 >>> from sympy.abc import x
325 >>> s = fourier_series(x**2, (x, -pi, pi))
326 >>> s.shift(1).truncate()
327 -4*cos(x) + cos(2*x) + 1 + pi**2/3
328 """
329 s, x = sympify(s), self.x
331 if x in s.free_symbols:
332 raise ValueError("'%s' should be independent of %s" % (s, x))
334 a0 = self.a0 + s
335 sfunc = self.function + s
337 return self.func(sfunc, self.args[1], (a0, self.an, self.bn))
339 def shiftx(self, s):
340 """
341 Shift x by a term independent of x.
343 Explanation
344 ===========
346 f(x) -> f(x + s)
348 This is fast, if Fourier series of f(x) is already
349 computed.
351 Examples
352 ========
354 >>> from sympy import fourier_series, pi
355 >>> from sympy.abc import x
356 >>> s = fourier_series(x**2, (x, -pi, pi))
357 >>> s.shiftx(1).truncate()
358 -4*cos(x + 1) + cos(2*x + 2) + pi**2/3
359 """
360 s, x = sympify(s), self.x
362 if x in s.free_symbols:
363 raise ValueError("'%s' should be independent of %s" % (s, x))
365 an = self.an.subs(x, x + s)
366 bn = self.bn.subs(x, x + s)
367 sfunc = self.function.subs(x, x + s)
369 return self.func(sfunc, self.args[1], (self.a0, an, bn))
371 def scale(self, s):
372 """
373 Scale the function by a term independent of x.
375 Explanation
376 ===========
378 f(x) -> s * f(x)
380 This is fast, if Fourier series of f(x) is already
381 computed.
383 Examples
384 ========
386 >>> from sympy import fourier_series, pi
387 >>> from sympy.abc import x
388 >>> s = fourier_series(x**2, (x, -pi, pi))
389 >>> s.scale(2).truncate()
390 -8*cos(x) + 2*cos(2*x) + 2*pi**2/3
391 """
392 s, x = sympify(s), self.x
394 if x in s.free_symbols:
395 raise ValueError("'%s' should be independent of %s" % (s, x))
397 an = self.an.coeff_mul(s)
398 bn = self.bn.coeff_mul(s)
399 a0 = self.a0 * s
400 sfunc = self.args[0] * s
402 return self.func(sfunc, self.args[1], (a0, an, bn))
404 def scalex(self, s):
405 """
406 Scale x by a term independent of x.
408 Explanation
409 ===========
411 f(x) -> f(s*x)
413 This is fast, if Fourier series of f(x) is already
414 computed.
416 Examples
417 ========
419 >>> from sympy import fourier_series, pi
420 >>> from sympy.abc import x
421 >>> s = fourier_series(x**2, (x, -pi, pi))
422 >>> s.scalex(2).truncate()
423 -4*cos(2*x) + cos(4*x) + pi**2/3
424 """
425 s, x = sympify(s), self.x
427 if x in s.free_symbols:
428 raise ValueError("'%s' should be independent of %s" % (s, x))
430 an = self.an.subs(x, x * s)
431 bn = self.bn.subs(x, x * s)
432 sfunc = self.function.subs(x, x * s)
434 return self.func(sfunc, self.args[1], (self.a0, an, bn))
436 def _eval_as_leading_term(self, x, logx=None, cdir=0):
437 for t in self:
438 if t is not S.Zero:
439 return t
441 def _eval_term(self, pt):
442 if pt == 0:
443 return self.a0
444 return self.an.coeff(pt) + self.bn.coeff(pt)
446 def __neg__(self):
447 return self.scale(-1)
449 def __add__(self, other):
450 if isinstance(other, FourierSeries):
451 if self.period != other.period:
452 raise ValueError("Both the series should have same periods")
454 x, y = self.x, other.x
455 function = self.function + other.function.subs(y, x)
457 if self.x not in function.free_symbols:
458 return function
460 an = self.an + other.an
461 bn = self.bn + other.bn
462 a0 = self.a0 + other.a0
464 return self.func(function, self.args[1], (a0, an, bn))
466 return Add(self, other)
468 def __sub__(self, other):
469 return self.__add__(-other)
472class FiniteFourierSeries(FourierSeries):
473 r"""Represents Finite Fourier sine/cosine series.
475 For how to compute Fourier series, see the :func:`fourier_series`
476 docstring.
478 Parameters
479 ==========
481 f : Expr
482 Expression for finding fourier_series
484 limits : ( x, start, stop)
485 x is the independent variable for the expression f
486 (start, stop) is the period of the fourier series
488 exprs: (a0, an, bn) or Expr
489 a0 is the constant term a0 of the fourier series
490 an is a dictionary of coefficients of cos terms
491 an[k] = coefficient of cos(pi*(k/L)*x)
492 bn is a dictionary of coefficients of sin terms
493 bn[k] = coefficient of sin(pi*(k/L)*x)
495 or exprs can be an expression to be converted to fourier form
497 Methods
498 =======
500 This class is an extension of FourierSeries class.
501 Please refer to sympy.series.fourier.FourierSeries for
502 further information.
504 See Also
505 ========
507 sympy.series.fourier.FourierSeries
508 sympy.series.fourier.fourier_series
509 """
511 def __new__(cls, f, limits, exprs):
512 f = sympify(f)
513 limits = sympify(limits)
514 exprs = sympify(exprs)
516 if not (isinstance(exprs, Tuple) and len(exprs) == 3): # exprs is not of form (a0, an, bn)
517 # Converts the expression to fourier form
518 c, e = exprs.as_coeff_add()
519 from sympy.simplify.fu import TR10
520 rexpr = c + Add(*[TR10(i) for i in e])
521 a0, exp_ls = rexpr.expand(trig=False, power_base=False, power_exp=False, log=False).as_coeff_add()
523 x = limits[0]
524 L = abs(limits[2] - limits[1]) / 2
526 a = Wild('a', properties=[lambda k: k.is_Integer, lambda k: k is not S.Zero, ])
527 b = Wild('b', properties=[lambda k: x not in k.free_symbols, ])
529 an = {}
530 bn = {}
532 # separates the coefficients of sin and cos terms in dictionaries an, and bn
533 for p in exp_ls:
534 t = p.match(b * cos(a * (pi / L) * x))
535 q = p.match(b * sin(a * (pi / L) * x))
536 if t:
537 an[t[a]] = t[b] + an.get(t[a], S.Zero)
538 elif q:
539 bn[q[a]] = q[b] + bn.get(q[a], S.Zero)
540 else:
541 a0 += p
543 exprs = Tuple(a0, an, bn)
545 return Expr.__new__(cls, f, limits, exprs)
547 @property
548 def interval(self):
549 _length = 1 if self.a0 else 0
550 _length += max(set(self.an.keys()).union(set(self.bn.keys()))) + 1
551 return Interval(0, _length)
553 @property
554 def length(self):
555 return self.stop - self.start
557 def shiftx(self, s):
558 s, x = sympify(s), self.x
560 if x in s.free_symbols:
561 raise ValueError("'%s' should be independent of %s" % (s, x))
563 _expr = self.truncate().subs(x, x + s)
564 sfunc = self.function.subs(x, x + s)
566 return self.func(sfunc, self.args[1], _expr)
568 def scale(self, s):
569 s, x = sympify(s), self.x
571 if x in s.free_symbols:
572 raise ValueError("'%s' should be independent of %s" % (s, x))
574 _expr = self.truncate() * s
575 sfunc = self.function * s
577 return self.func(sfunc, self.args[1], _expr)
579 def scalex(self, s):
580 s, x = sympify(s), self.x
582 if x in s.free_symbols:
583 raise ValueError("'%s' should be independent of %s" % (s, x))
585 _expr = self.truncate().subs(x, x * s)
586 sfunc = self.function.subs(x, x * s)
588 return self.func(sfunc, self.args[1], _expr)
590 def _eval_term(self, pt):
591 if pt == 0:
592 return self.a0
594 _term = self.an.get(pt, S.Zero) * cos(pt * (pi / self.L) * self.x) \
595 + self.bn.get(pt, S.Zero) * sin(pt * (pi / self.L) * self.x)
596 return _term
598 def __add__(self, other):
599 if isinstance(other, FourierSeries):
600 return other.__add__(fourier_series(self.function, self.args[1],\
601 finite=False))
602 elif isinstance(other, FiniteFourierSeries):
603 if self.period != other.period:
604 raise ValueError("Both the series should have same periods")
606 x, y = self.x, other.x
607 function = self.function + other.function.subs(y, x)
609 if self.x not in function.free_symbols:
610 return function
612 return fourier_series(function, limits=self.args[1])
615def fourier_series(f, limits=None, finite=True):
616 r"""Computes the Fourier trigonometric series expansion.
618 Explanation
619 ===========
621 Fourier trigonometric series of $f(x)$ over the interval $(a, b)$
622 is defined as:
624 .. math::
625 \frac{a_0}{2} + \sum_{n=1}^{\infty}
626 (a_n \cos(\frac{2n \pi x}{L}) + b_n \sin(\frac{2n \pi x}{L}))
628 where the coefficients are:
630 .. math::
631 L = b - a
633 .. math::
634 a_0 = \frac{2}{L} \int_{a}^{b}{f(x) dx}
636 .. math::
637 a_n = \frac{2}{L} \int_{a}^{b}{f(x) \cos(\frac{2n \pi x}{L}) dx}
639 .. math::
640 b_n = \frac{2}{L} \int_{a}^{b}{f(x) \sin(\frac{2n \pi x}{L}) dx}
642 The condition whether the function $f(x)$ given should be periodic
643 or not is more than necessary, because it is sufficient to consider
644 the series to be converging to $f(x)$ only in the given interval,
645 not throughout the whole real line.
647 This also brings a lot of ease for the computation because
648 you do not have to make $f(x)$ artificially periodic by
649 wrapping it with piecewise, modulo operations,
650 but you can shape the function to look like the desired periodic
651 function only in the interval $(a, b)$, and the computed series will
652 automatically become the series of the periodic version of $f(x)$.
654 This property is illustrated in the examples section below.
656 Parameters
657 ==========
659 limits : (sym, start, end), optional
660 *sym* denotes the symbol the series is computed with respect to.
662 *start* and *end* denotes the start and the end of the interval
663 where the fourier series converges to the given function.
665 Default range is specified as $-\pi$ and $\pi$.
667 Returns
668 =======
670 FourierSeries
671 A symbolic object representing the Fourier trigonometric series.
673 Examples
674 ========
676 Computing the Fourier series of $f(x) = x^2$:
678 >>> from sympy import fourier_series, pi
679 >>> from sympy.abc import x
680 >>> f = x**2
681 >>> s = fourier_series(f, (x, -pi, pi))
682 >>> s1 = s.truncate(n=3)
683 >>> s1
684 -4*cos(x) + cos(2*x) + pi**2/3
686 Shifting of the Fourier series:
688 >>> s.shift(1).truncate()
689 -4*cos(x) + cos(2*x) + 1 + pi**2/3
690 >>> s.shiftx(1).truncate()
691 -4*cos(x + 1) + cos(2*x + 2) + pi**2/3
693 Scaling of the Fourier series:
695 >>> s.scale(2).truncate()
696 -8*cos(x) + 2*cos(2*x) + 2*pi**2/3
697 >>> s.scalex(2).truncate()
698 -4*cos(2*x) + cos(4*x) + pi**2/3
700 Computing the Fourier series of $f(x) = x$:
702 This illustrates how truncating to the higher order gives better
703 convergence.
705 .. plot::
706 :context: reset
707 :format: doctest
708 :include-source: True
710 >>> from sympy import fourier_series, pi, plot
711 >>> from sympy.abc import x
712 >>> f = x
713 >>> s = fourier_series(f, (x, -pi, pi))
714 >>> s1 = s.truncate(n = 3)
715 >>> s2 = s.truncate(n = 5)
716 >>> s3 = s.truncate(n = 7)
717 >>> p = plot(f, s1, s2, s3, (x, -pi, pi), show=False, legend=True)
719 >>> p[0].line_color = (0, 0, 0)
720 >>> p[0].label = 'x'
721 >>> p[1].line_color = (0.7, 0.7, 0.7)
722 >>> p[1].label = 'n=3'
723 >>> p[2].line_color = (0.5, 0.5, 0.5)
724 >>> p[2].label = 'n=5'
725 >>> p[3].line_color = (0.3, 0.3, 0.3)
726 >>> p[3].label = 'n=7'
728 >>> p.show()
730 This illustrates how the series converges to different sawtooth
731 waves if the different ranges are specified.
733 .. plot::
734 :context: close-figs
735 :format: doctest
736 :include-source: True
738 >>> s1 = fourier_series(x, (x, -1, 1)).truncate(10)
739 >>> s2 = fourier_series(x, (x, -pi, pi)).truncate(10)
740 >>> s3 = fourier_series(x, (x, 0, 1)).truncate(10)
741 >>> p = plot(x, s1, s2, s3, (x, -5, 5), show=False, legend=True)
743 >>> p[0].line_color = (0, 0, 0)
744 >>> p[0].label = 'x'
745 >>> p[1].line_color = (0.7, 0.7, 0.7)
746 >>> p[1].label = '[-1, 1]'
747 >>> p[2].line_color = (0.5, 0.5, 0.5)
748 >>> p[2].label = '[-pi, pi]'
749 >>> p[3].line_color = (0.3, 0.3, 0.3)
750 >>> p[3].label = '[0, 1]'
752 >>> p.show()
754 Notes
755 =====
757 Computing Fourier series can be slow
758 due to the integration required in computing
759 an, bn.
761 It is faster to compute Fourier series of a function
762 by using shifting and scaling on an already
763 computed Fourier series rather than computing
764 again.
766 e.g. If the Fourier series of ``x**2`` is known
767 the Fourier series of ``x**2 - 1`` can be found by shifting by ``-1``.
769 See Also
770 ========
772 sympy.series.fourier.FourierSeries
774 References
775 ==========
777 .. [1] https://mathworld.wolfram.com/FourierSeries.html
778 """
779 f = sympify(f)
781 limits = _process_limits(f, limits)
782 x = limits[0]
784 if x not in f.free_symbols:
785 return f
787 if finite:
788 L = abs(limits[2] - limits[1]) / 2
789 is_finite, res_f = finite_check(f, x, L)
790 if is_finite:
791 return FiniteFourierSeries(f, limits, res_f)
793 n = Dummy('n')
794 center = (limits[1] + limits[2]) / 2
795 if center.is_zero:
796 neg_f = f.subs(x, -x)
797 if f == neg_f:
798 a0, an = fourier_cos_seq(f, limits, n)
799 bn = SeqFormula(0, (1, oo))
800 return FourierSeries(f, limits, (a0, an, bn))
801 elif f == -neg_f:
802 a0 = S.Zero
803 an = SeqFormula(0, (1, oo))
804 bn = fourier_sin_seq(f, limits, n)
805 return FourierSeries(f, limits, (a0, an, bn))
806 a0, an = fourier_cos_seq(f, limits, n)
807 bn = fourier_sin_seq(f, limits, n)
808 return FourierSeries(f, limits, (a0, an, bn))