CODEWITHSUNDEEP
pythonpython-3.xfftsympypiecewise
alienflow
alienflow

Reputation: 410

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

Upvotes: 3

Views: 2622

Answers (2)

Helmut Pfeffer
Helmut Pfeffer

Reputation: 1

The piecewise evaluate the conds returning the first that is True. In your definition of p, it evaluates first if 0 < t, and when true, it stops and returns sin(t). So for every positive t, you get sin(t).

Better would be just to put the second condition of each interval, and leave the first one out:

p = Piecewise((sin(t), t < pi  ),
              (0,      t < 2*pi))

Upvotes: 0

hpaulj
hpaulj

Reputation: 231510

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).

Upvotes: 3

Related Questions