Ряд Фурье в Numpy. Вопрос о предыдущем ответе

Question

Я пытаюсь повторить ответ, данный в предыдущей теме: Как вычислить ряд Фурье в Numpy?

import numpy as np
import matplotlib.pyplot as plt
import itertools

def func(x):
    if x >= 1.0 or x <= -1.0:
        return 0
    else:
        return (abs(x) - 1.0)
a = 1.0
b = -1.0
N = 128.
time = np.linspace( a, b, N )
y = (np.fromiter(itertools.imap(func, time), 
                  dtype=time.dtype, count=time.shape[0]))
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(time,y) 
period = 2.
def cn(n):
    c = y*np.exp(-1j*2*n*np.pi*time/period)
    return c.sum()/c.size
def f(x, Nh):
    f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/period) for i in range(1,Nh+1)])
    return f.sum()
y2 = np.array([f(t,10).real for t in time])
ax.plot(time, y2)
plt.show()

Я получаю решение, близкое кправильный ответ, но сместился.Я не был уверен, что я делаю неправильно.

Answer 1

Ошибка, по-видимому, связана с вашей суммой Римана метод (справа / посередине / слева) - обозначается регулярный фрай . Использование среднего метода дает:

Код:

import numpy as np
import matplotlib.pyplot as plt
import itertools

def func(x):
    if x >= 1.0 or x <= -1.0:
        return 0
    else:
        return (abs(x) - 1.0)

a = 1.0
b = -1.0
N = 128.
time = np.linspace( a, b, N )
y = (np.fromiter(itertools.imap(func, time), 
                  dtype=time.dtype, count=time.shape[0]))

period = 2.
def cn(n):
    c = y*np.exp(-1j*2*n*np.pi*time/period)
    return c.sum()/c.size
def f(x, Nh):
    rng = np.arange(.5, Nh+.5)
    f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/period) for i in rng])
    return f.sum()

y2 = np.array([f(t,10).real for t in time])

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(time, y)
ax.plot(time, y2)
plt.show()

Как отметил Свен в другом вопросе использование вами списочных представлений (и imap) вместо массивов и ufuncs довольно неэффективно (если вы столкнетесь с проблемами производительности)

Answer 2

Похоже, ваш термин в округе Колумбия где-то потерян. Я не могу проверить себя сейчас, но вы уверены, что 1 в диапазоне (1, Nh + 1) в f () правильно?

Answer 3

Векторизованная версия вашего кода:

import numpy as np
import matplotlib.pyplot as plt
from optparse import OptionParser

def func(x):
    return np.where(np.abs(x) >= 1, 0., np.abs(x) - 1.0)

def cn(x, y, n, period):
    c = y * np.exp(-1j * 2. * np.pi * n * x / period)
    return c.sum()/c.size

def f(x, y, Nh, period):
    rng = np.arange(.5, Nh+.5)
    coeffs = np.array([cn(x,y,i,period) for i in rng])
    f = np.array([2. * coeffs[i] * np.exp(1j*2*i*np.pi*x/period) for i in rng])
    return f.sum(axis=0)

if __name__=='__main__':

    Version = '0.1'
    usage = "usage: %prog [options]"

    parser = OptionParser(usage = usage,version="%prog "+Version)
    parser.add_option("-a", dest='a', type='float', default=1., help="initial time")
    parser.add_option("-b", dest='b', type='float', default=-1., help="end time")
    parser.add_option("-N", "--Nt", dest='N', type='int', default=128, help="number of time steps")
    parser.add_option("-p", "--period", dest='period', type='float', default=2., help="period [time span]")
    parser.add_option("--Nh", dest='Nh', type='int', default=10, help="number of fourier series terms")

    (options, args) = parser.parse_args()

    for key,value in options.__dict__.iteritems():
        exec key + ' = ' + repr(value)

    time = np.linspace( a, b, N )
    y = func(time)
    period = np.abs(a-b)

    y2 = f(time,y,Nh,period).real

    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    ax.plot(time, y)
    ax.plot(time, y2)
    plt.show()

Учитывая, что код был сохранен с именем "fourier_series.py", вы можете попробовать:

python fourier_series.py -N 512 --Nh 128

в обычном терминале или:

% run fourier_series.py -N 512 --Nh 128

в консоли ipython