Search code examples
pythonpython-3.xfftsympypiecewise

Fourier series of Piecewise PYTHON

I tried to find the Fourier Series of

piecewise

With simpy like :

p = Piecewise((sin(t), 0 < t),(sin(t), t < pi), (0 , pi < t), (0, t < 2*pi))
fs = fourier_series(p, (t, 0, 2*pi)).truncate(8)

But it doesn't seem to work. It is stuck in * (looping?). Is there any way to solve that? Perhaps an alternative? Many thanks

Solution

  • I get, with a second or two of delay:

    In [55]: fourier_series(p,(t,0,2*pi))
Out[55]: FourierSeries(Piecewise((sin(t), (t > 0) | (t < pi)), (0, (pi < t) | (t < 2*pi))), (t, 0, 2*pi), (0, SeqFormula(Piecewise((0, Eq(_n, -1) | Eq(_n, 1)), (cos(2*_n*pi)/(_n**2 - 1) - 1/(_n**2 - 1), True))*cos(_n*t)/pi, (_n, 1, oo)), SeqFormula(Piecewise((-pi, Eq(_n, -1)), (pi, Eq(_n, 1)), (sin(2*_n*pi)/(_n**2 - 1), True))*sin(_n*t)/pi, (_n, 1, oo))))

    That's just setting it up.

    _.truncate(8) is taking (too) long. That must be doing the evaluation.

    Does a different truncation work better? I don't see any other controls.

    .truncate(1) returns sin(t). .truncate(2) hangs. Mixing this simple sin(t) with a flat segment must be setting up a difficult case that is analytically difficult. But I'm a bit rusty on this area of math.

    Looking for fourier series with numpy I found:

    How to calculate a Fourier series in Numpy?

    For a FS defined on (0,pi) fs1 = fourier_series(p, (t, 0, pi)):

    In [5]: fs1.truncate(1)
Out[5]: 2/pi
In [6]: fs1.truncate(2)
Out[6]: -4*cos(2*t)/(3*pi) + 2/pi
In [7]: fs1.truncate(3)
Out[7]: -4*cos(2*t)/(3*pi) - 4*cos(4*t)/(15*pi) + 2/pi
In [8]: fs1.truncate(4)
Out[8]: -4*cos(2*t)/(3*pi) - 4*cos(4*t)/(15*pi) - 4*cos(6*t)/(35*pi) + 2/pi
In [9]: fs1.truncate(5)
Out[9]: -4*cos(2*t)/(3*pi) - 4*cos(4*t)/(15*pi) - 4*cos(6*t)/(35*pi) - 4*cos(8*t)/(63*pi) + 2/pi

    Which plot (in numpy) as expected:

    periodic 0-pi FS

    From a table of Fourier Series, I found this formula (in numpy terms) for a rectified sine wave:

    z8 = 1/pi + 1/2*sin(t)-2/pi*np.sum([cos(2*i*t)/(4*i**2-1) for i in range(1,8)],axis=0)

    rectified half sine wave

    This has a similar cos series term, but adds that sin term. That suggests to me that you could approximate this half sin as a sum of a*sin(t)+b(sin(2*t)) (or something like that). I imagine that there are math texts or papers that explore the difficulties in deriving fourier series as sympy does. Have you looked at the Mathworld link?

    Comparing the FS for a rectified half sine with a rectified whole sine

    half sine:

    In [434]: z3 = 1/pi + 1/2*sin(t)-2/pi*np.sum([cos(2*i*t)/(4*i**2-1) for i in range(1,3)],axis=0)

    full sine:

    In [435]: w3 = 1/pi -2/pi*np.sum([cos(2*i*t)/(4*i**2-1) for i in range(1,3)],axis=0)

In [438]: plt.plot(t,sin(t)/2)
In [439]: plt.plot(t,w3)
In [440]: plt.plot(t,z3)
In [441]: plt.plot(t,w3+sin(t)/2)  # full sine + sine/2 == half sine

    rectified half sine vs full sine

    I can imagine transfering insights like this back into sympy, redefining the periodic function in a way that doesn't take so long (or possibly hang).