Oblicz szereg Fouriera z podejściem trygonometryczne

Question

staram się realizować funkcję szeregu Fouriera według następujących formuł:

... Gdzie ...

... i ...

Oto moje podejście do problemu:

import numpy as np 
import pylab as py 

# Define "x" range. 
x = np.linspace(0, 10, 1000) 

# Define "T", i.e functions' period. 
T = 2 
L = T/2 

# "f(x)" function definition. 
def f(x): 
    return np.sin(np.pi * 1000 * x) 

# "a" coefficient calculation. 
def a(n, L, accuracy = 1000): 
    a, b = -L, L 
    dx = (b - a)/accuracy 
    integration = 0 
    for i in np.linspace(a, b, accuracy): 
     x = a + i * dx 
     integration += f(x) * np.cos((n * np.pi * x)/L) 
    integration *= dx 
    return (1/L) * integration 

# "b" coefficient calculation. 
def b(n, L, accuracy = 1000): 
    a, b = -L, L 
    dx = (b - a)/accuracy 
    integration = 0 
    for i in np.linspace(a, b, accuracy): 
     x = a + i * dx 
     integration += f(x) * np.sin((n * np.pi * x)/L) 
    integration *= dx 
    return (1/L) * integration 

# Fourier series. 
def Sf(x, L, n = 10): 
    a0 = a(0, L) 
    sum = 0 
    for i in np.arange(1, n + 1): 
     sum += ((a(i, L) * np.cos(n * np.pi * x)) + (b(i, L) * np.sin(n * np.pi * x))) 
    return (a0/2) + sum  

# x axis. 
py.plot(x, np.zeros(np.size(x)), color = 'black') 

# y axis. 
py.plot(np.zeros(np.size(x)), x, color = 'black') 

# Original signal. 
py.plot(x, f(x), linewidth = 1.5, label = 'Signal') 

# Approximation signal (Fourier series coefficients). 
py.plot(x, Sf(x, L), color = 'red', linewidth = 1.5, label = 'Fourier series') 

# Specify x and y axes limits. 
py.xlim([0, 10]) 
py.ylim([-2, 2]) 

py.legend(loc = 'upper right', fontsize = '10') 

py.show() 

... i oto co mam po kreślenia Rezultat:

Czytałam How to calculate a Fourier series in Numpy? i już wdrożyłem to podejście. Działa wspaniale, ale wykorzystuje metodę ekspotencjalną, w której chcę skupić się na funkcjach trygonometrycznych i metodzie prostokątnej w przypadku obliczania współczynników całkowania dla współczynników a_{n} i b_{n}.

Z góry dziękuję.

UPDATE (SOLVED)

Wreszcie, tutaj jest przykład pracy z kodem. Jednak poświęcę na to więcej czasu, więc jeśli jest coś, co można poprawić, będzie to zrobione.

from __future__ import division 
import numpy as np 
import pylab as py 

# Define "x" range. 
x = np.linspace(0, 10, 1000) 

# Define "T", i.e functions' period. 
T = 2 
L = T/2 

# "f(x)" function definition. 
def f(x): 
    return np.sin((np.pi) * x) + np.sin((2 * np.pi) * x) + np.sin((5 * np.pi) * x) 

# "a" coefficient calculation. 
def a(n, L, accuracy = 1000): 
    a, b = -L, L 
    dx = (b - a)/accuracy 
    integration = 0 
    for x in np.linspace(a, b, accuracy): 
     integration += f(x) * np.cos((n * np.pi * x)/L) 
    integration *= dx 
    return (1/L) * integration 

# "b" coefficient calculation. 
def b(n, L, accuracy = 1000): 
    a, b = -L, L 
    dx = (b - a)/accuracy 
    integration = 0 
    for x in np.linspace(a, b, accuracy): 
     integration += f(x) * np.sin((n * np.pi * x)/L) 
    integration *= dx 
    return (1/L) * integration 

# Fourier series. 
def Sf(x, L, n = 10): 
    a0 = a(0, L) 
    sum = np.zeros(np.size(x)) 
    for i in np.arange(1, n + 1): 
     sum += ((a(i, L) * np.cos((i * np.pi * x)/L)) + (b(i, L) * np.sin((i * np.pi * x)/L))) 
    return (a0/2) + sum 

# x axis. 
py.plot(x, np.zeros(np.size(x)), color = 'black') 

# y axis. 
py.plot(np.zeros(np.size(x)), x, color = 'black') 

# Original signal. 
py.plot(x, f(x), linewidth = 1.5, label = 'Signal') 

# Approximation signal (Fourier series coefficients). 
py.plot(x, Sf(x, L), '.', color = 'red', linewidth = 1.5, label = 'Fourier series') 

# Specify x and y axes limits. 
py.xlim([0, 5]) 
py.ylim([-2.2, 2.2]) 

py.legend(loc = 'upper right', fontsize = '10') 

py.show() 

Answer 1

rozważyć opracowanie kodu w inny sposób, blok po bloku. Powinieneś być zaskoczony, jeśli taki kod zadziała przy pierwszej próbie. Debugowanie jest jedną z opcji, jak powiedział @ tom10. Inną opcją jest szybkie prototypowanie kodu krok po kroku w tłumaczu, jeszcze lepiej z ipythonem.

Powyżej, oczekujesz, że b_1000 ma wartość niezerową, ponieważ wejście f(x) jest sinusoidą z 1000 w nim. Oczekujesz również, że wszystkie pozostałe współczynniki wynoszą zero, prawda?

Następnie należy skoncentrować się tylko na funkcji b(n, L, accuracy = 1000). Patrząc na to, 3 rzeczy idą źle. Oto kilka wskazówek.

• mnożenie dx jest w pętli. Pewny tego?
• w pętli, i ma być liczbą całkowitą w prawo? Czy to naprawdę jest liczba całkowita? poprzez prototypowanie lub debugowanie można by się przekonać, że przy każdym pisaniu (1/L) lub podobnym wyrażeniu należy zachować ostrożność. Jeśli używasz python2.7, robisz prawdopodobnie źle. Jeśli nie, użyj przynajmniej numeru from __future__ import division u góry źródła. Przeczytaj this PEP, jeśli nie wiesz o czym mówię.

Jeśli rozwiążesz te 3 punkty, zadziała b(). Następnie pomyśl o a w podobny sposób.