Я получаю, с секундой или двумя задержками:
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))))
Это просто настройка.
_.truncate(8)
занимает (слишком) много времени.Это должно делать оценку.
Работает ли другое усечение лучше?Других элементов управления я не вижу.
.truncate(1)
возвращает sin(t)
..truncate(2)
виснет.Смешивание этого простого sin(t)
с плоским сегментом должно создать сложный случай, который является аналитически сложным.Но я немного заржавел в этой области математики.
В поисках ряда Фурье с numpy я нашел:
Как вычислить ряд Фурье в Numpy?
Для FS, определенной на (0, пи) 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
Какой участок (в виде цифры), как и ожидалось:
Из таблицы ряд Фурье я нашел эту формулу (в терминах numpy
) для выпрямленной синусоидальной волны:
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)
Этот термин имеет аналогичный термин cos
, но добавляет термин sin
.Это наводит меня на мысль, что вы можете приблизить эту половину греха к сумме a*sin(t)+b(sin(2*t))
(или что-то в этом роде).Я предполагаю, что есть математические тексты или статьи, в которых исследуются трудности получения рядов Фурье, как это делает sympy
.Вы смотрели на ссылку Mathworld?
Сравнение FS для выпрямленного полусинуса с выпрямленным полным синусом
полусинус:
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)
полный синус:
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
Я могу представить себе передачу таких идей обратно в sympy
, переопределение периодической функции таким образом, чтобы это не занимало много времени (или возможно повесить).