Calcolare la serie di Fourier con l'approccio trigonometria

Question

provo ad implementare la funzione serie di Fourier secondo le seguenti formule:

... dove ...

... e ...

Ecco il mio approccio al problema:

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

... e qui è quello che ottengo dopo aver tracciato il risultato:

ho letto il How to calculate a Fourier series in Numpy? e ho già implementato questo approccio. Funziona alla grande, ma usa il metodo espotenziale, dove voglio concentrarmi sulle funzioni di trigonometria e il metodo rettangolare nel caso di calcolare le integrazioni per i coefficienti a_{n} e b_{n}.

Grazie in anticipo.

UPDATE (RISOLTO)

Infine, ecco un esempio operativo del codice. Tuttavia, passerò più tempo su di esso, quindi se c'è qualcosa che può essere migliorato, sarà fatto.

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

considerazione lo sviluppo il codice in un modo diverso, blocco per blocco. Dovresti essere sorpreso se un codice come questo funzionerebbe al primo tentativo. Il debug è un'opzione, come ha affermato @ tom10. L'altra opzione è la prototipazione rapida del codice passo dopo passo nell'interprete, ancora meglio con ipython.

Sopra, ci si aspetta che b_1000 sia diverso da zero, poiché l'ingresso f(x) è una sinusoide con un 1000 in esso. Ti stai anche aspettando che tutti gli altri coefficienti siano pari a zero?

Quindi è necessario concentrarsi solo sulla funzione b(n, L, accuracy = 1000). Guardandolo, 3 cose stanno andando male. Ecco alcuni suggerimenti.

• la moltiplicazione di dx è all'interno del ciclo. Sicuro su questo?
• nel loop, i dovrebbe essere un numero intero giusto? È davvero un numero intero? con la prototipazione o il debug scoprirai che questo
• stai attento quando scrivi (1/L) o un'espressione simile. Se stai usando python2.7, probabilmente stai sbagliando. In caso contrario, utilizzare almeno uno from __future__ import division nella parte superiore dell'origine. Leggi this PEP se non sai di cosa sto parlando.

Se si indirizza questi 3 punti, b() funzionerà. Quindi pensa a a in un modo simile.