Serie di Fourier Matlab

[Edo_Rm]

Salve a tutti!
Prima di tutto mi presento dato che è il mio primo post: mi chiamo Edoardo e sono al 4° anno di università, studio ingegneria chimica.

I corsi di Matlab che ho seguito sono sempre stati specifici sugli argomenti inerenti l'ingegneria chimica, e per conto mio sto cercando di imparare ad usarlo al meglio anche per altri scopi, purtroppo da autodidatta avanzo lentamente...

L'ultimo problema che mi ero posto era di scrivere un programma che mi calcolasse l'espansione in serie di Fourier di una generica funzione f(x).
Incollo qui il programma che ho scritto (mi scuso per quanto possa essere "rozzo"):

clc
close all
clear all

syms x

L=2; %voglio valutare la funzione x^2 nell'intervallo [-1,1], quindi periodo L=2
y=x^2;

%Qui calcolo i coefficienti Ak, metto dentro il ciclo anche A0 tanto non dipende da k e non cambia nulla, essendo la funzione f(x)=x^2 pari, i coefficienti Bk sono nulli.

for k=1:100
a0=(2/L)*int(y,-1,1);
a(k)=(2/L)*int(y*cos(2*pi*k*x/L),-1,1);
end

%Qui calcolo i coseni che poi dovrò implementare nella serie di Fourier

for k=1:100
c(k)=cos(2*pi*k*x/L);
end

%Qui calcolo i singoli elementi della serie di Fourier

for k=1:100
fourier(k)=a(k)*c(k);
end

%$A_0$/2+$\sum_{k=1}^N $A_k$*cos(2*k*pi*x/L)$ Qui sommo tutti i singoli elementi della serie di Fourier

Fourier=0;
for k=1:100
Fourier=(Fourier+fourier(k))+(a0/2);
end


Risultato=Fourier

%Qui verifico la validità della mia serie di Fourier vedendo che risultato mi da per x=0.3, dovrebbe venire qualcosa di simile a 0.09 ma invece mi viene un 33.09...

x=0.3;

Risultatobis=Fourier



Ora mi chiedo, dov'è l'intoppo del programma che scrivo? Inoltre, come mai il programma è abbastanza lento a girare?

Grazie a tutti per l'aiuto che mi darete!

[apatriarca]

metto dentro il ciclo anche A0 tanto non dipende da k e non cambia nulla

Ma certo che cambia!!!! Stai calcolando 100 volte quello che potrebbe essere calcolato una sola volta e non sono certo che Matlab sia abbastanza intelligente da estrarlo dal ciclo. Sono anzi abbastanza certo che non sia in grado di farlo con un'operazione complessa come quella che stai facendo. In ogni caso il calcolare un integrale in modo simbolico è sempre lento. Considerando che alla fine vuoi calcolare dei numeri, perché non calcoli direttamente gli integrali in modo numerico? Un'altra ragione per la lentezza del tuo codice è che matlab è più veloce quando si scrive tutto in modo vettoriale invece di usare cicli. Per cui il ciclo che calcola il vettore fourier si poteva riscrivere con la seguente formula:

fourier = a .* c;

o ancora meglio si potevano unire gli ultimi due cicli osservando che stavi semplicemente calcolando il prodotto scalare tra i vettori a e c e aggiungendo il termine a0/2. Per cui sarebbe stato meglio scrivere

Fourier = a0/2 + a' * c;

Un'altra nota in fatto di efficienza, riguarda l'inizializzazione dei vettori che vengono usati. È più efficiente infatti informare matlab della dimensione di una matrice o vettore prima di usarlo in modo che allochi direttamente la memoria necessaria.

Per quanto riguarda la correttezza, l'errore è invece probabilmente nella seguente riga:

Fourier=(Fourier+fourier(k))+(a0/2);

Stai sommando (a0/2) ad ogni iterazione del ciclo, cioè 100 volte!!! Non credo che il tuo obiettivo fosse quello...

[Edo_Rm]

Grazie mille per la risposta tempestiva!
Il problema era esattamente dove dicevi te! Ho modificato il codice in questo modo:

Fourier=0;
for k=1:10
Fourier=(Fourier+fourier(k));
end

Fourier=Fourier+(a0/2);

Così che mi conti una sola volta il termine a0/2.

Per quanto riguarda la lentezza ho provato a modificare i cicli for, ma togliendo i ";" mi sono accorto che la parte più lenta e pesante del programma è il calcolo dell'integrale per il coefficiente Ak, quindi anche fondendo i due cicli for insieme non cambia nulla!

L'integrale indefinito lo voglio lasciare dato che voglio che il programma sia in grado di calcolare lo sviluppo in serie di Fourier di una qualsiasi funzione pari, modificando semplicemente la y all'inizio e l'ampiezza del periodo L.

Grazie ancora

[apatriarca]

Non vedo come tuo possa riutilizzare il risultato per altre funzioni. I coefficienti che ti calcoli nel ciclo dipendono dalla funzione, non sono uguali per altre funzioni pari. Quindi non credo che cambi qualcosa utilizzando il tuo metodo al posto di una qualche versione numerica.

[Edo_Rm]

"apatriarca":Non vedo come tuo possa riutilizzare il risultato per altre funzioni. I coefficienti che ti calcoli nel ciclo dipendono dalla funzione, non sono uguali per altre funzioni pari. Quindi non credo che cambi qualcosa utilizzando il tuo metodo al posto di una qualche versione numerica.

Beh ma cambiando la funzione y= mi cambiano i coefficienti Ak e di conseguenza cambia tutto lo sviluppo.
Tenendo lo stesso programma ho sostituito x^2 con abs(x) e gira lo stesso!

[apatriarca]

Ma questo sarebbe vero anche nel caso in cui risolvessi il problema in modo numerico. Sarebbe sufficiente scrivere

y = @(x) x^2;

o

y = @(x) abs(x);

al posto della tua definizione della funzione e poi usare quad o quad8 al posto di int per risolverlo in modo numerico.

[Guerz1]

ho trovato questo post molto interessante perche approfondisce un esame che sto preparando proprio in questo periodo. Avrei però una domanda molto stupida da fare: come faccio a vedere i valori che restituisce l'approssimazione della serie di Fourier? (ad esempio per poter fare un grafico dell'errore plottando f(x)-fourier(x) )

grazie 1000!