Alcuni sviluppi in serie di Fourier

Riferimento : Matematica per l' Ingegneria dell' Informazione, par. 3.2

Jean Baptiste Joseph Fourier (1768-1830)

Lo scopo di questo notebook è principalmente quello di mostrare  i polinomi di Fourier per alcune classiche funzioni periodiche, e di verificare che la successione dei coefficienti di Fourier tende a 0  tanto più rapidamente quanto più la funzione è regolare.

Clear["Global`"]

Costruzione di un diagramma a barre:

diagBarre[lista_List,n_:1]:=Module[{nm1=n-1},Show[Table[Graphics[{Thickness[0.009],RGBColor[1,0,0],Line[{{i+nm1,0},{i+nm1,lista[[i]]}}]}],{i,Length[lista]}],PlotRange->{{0,Length[lista]+1},{Min[lista]-1,Max[lista]+1}},Ticks->{Automatic,Range[Floor[Min[lista]],Ceiling[Max[lista]]]},Axes->True,AxesOrigin->{0,0},ImageSize->{400,200},BaseStyle->{FontSize->12}]];

Esempio 1: la funzione identità sull'intervallo [–π, π], prolungata con periodo 2 π

(∫_0^π)x Sin[k x]\[DifferentialD]x

Simplify[%, k ∈ Integers]

Polinomi di Fourier:

s1[n_,x_]:=2  Sum[(-1)^(k+1) Sin[k x]/k, {k,1,n}]

Funzione che viene sviluppata:

f1[x_]:=If[Abs[x]<π,x,f1[Mod[x+π,2 π]-π]]

Alcuni grafici (monometrici):

Manipulate[Plot[{f1[x],s1[n,x]}, {x,-5,5},
		AspectRatio->Automatic,
	PlotRange->{{-5,5},{-4,4}},
	PlotStyle->{{RGBColor[0,0,1],Thickness[0.012]},  {RGBColor[1,0,0],Thickness[0.01]}}],{n,1,10,1}]

Spettro di ampiezza:

amp1 = Table[2/k, {k,1,20}];

diagBarre[amp1]

Esempio 2: la funzione segno sull'intervallo [-π, π], prolungata con periodo 2 π

a=.;b=.;

(∫_0^π)Sin[k x]\[DifferentialD]x

Simplify[%,k∈Integers]

Polinomi di Fourier:

s2[n_,x_]:=(4*Sum[Sin[(2*k+1)*x]/(2*k+1),{k,0,n}])/π

Funzione che viene sviluppata:

f2[x_]:=If[Abs[x]<π,Sign[x],f2[Mod[x+π,2 π]-π]]

Alcuni grafici (monometrici):

Manipulate[Plot[{f2[x],s2[n,x]},{x,-5,5},AspectRatio->Automatic,PlotRange->{{-4,4},{-1.5,1.5}},PlotStyle->{{RGBColor[0,0,1],Thickness[0.012]},{RGBColor[1,0,0],Thickness[0.01]}}],{n,0,6,1}]

coeff2=Table[((2 (1+(((-1))^k)))/(π k)),{k,1,20}]//N;

Spettro di ampiezza:

diagBarre[coeff2]

Esempio 3: la funzione valore assoluto di x sull'intervallo [–π, π], prolungata con periodo 2 π

a=.; b=.;

Commento: si tratta di una funzione continua su R ; i coefficienti di Fourier tendono a 0 come 1/(k^2).

(∫_0^π)x Cos[k x]\[DifferentialD]x

Simplify[%, k∈Integers]

Polinomi di Fourier (si tratta dei polinomi di indice 2n +1):

s3[n_, x_] := π/2 - (4*Sum[Cos[(2*k + 1)*x]/(2*k + 1)^2, {k, 0, n}])/π

Funzione che viene sviluppata:

f3[x_]:=If[Abs[x]<π,Abs[x],f3[Mod[x+π,2 π]-π]]

Alcuni grafici (monometrici):

Manipulate[Plot[{f3[x],s3[n,x]}, {x,-5,5},
		AspectRatio->Automatic,
	PlotRange->{{-5,5},{-4,4}},
	PlotStyle->{{RGBColor[0,0,1],Thickness[0.01]},  {RGBColor[1,0,0],Thickness[0.01]}}],{n,0,3,1}]

Tracciamo il grafico dell'errore (funzione data meno polinomio di Fourier del terzo ordine):

Plot[{f3[x]-s3[3,x]},{x,-4,4}]

Attenzione alle unità di misura!
Spettro di ampiezza:

amp3=N[Table[((2 Abs[((-1+(((-1))^k))/(k^2))])/π),{k,1,20}]];

amp3=Join[{N[(π/2)]},amp3];

diagBarre[amp3,0]

Esempio 4: la funzione x -> max{0, x} sull'intervallo [–π, π], 
prolungata con periodo 2 π

a=.; b=.;

Calcoliamo direttamente a[0] per evitare formula in cui k compare a denominatore:

(∫_0^π)t \[DifferentialD]t

a[0] = %/π

Simplify[(∫_0^π)t Cos[k t]\[DifferentialD] t, k∈Integers]

a[k_]:=((-1+(((-1))^k))/(π (k^2)))

Simplify[(∫_0^π)t Sin[k t]\[DifferentialD] t, k∈Integers]

b[k_]:= -(((((-1))^k) π)/(π k))

Polinomi di Fourier:

s4[n_,x_]:= a[0]/2 + Sum[a[k]  Cos[ k x] + b[k]  Sin[k x], {k,1,n}]

Funzione che viene sviluppata:

f4[x_]:=If[Abs[x]<π,Max[0,x],f4[Mod[x+π,2 π]-π]]

Alcuni grafici (monometrici):

Manipulate[Plot[{f4[x], s4[n,x]}, {x,-5,5},
AspectRatio->Automatic,
PlotStyle->{{RGBColor[0,0,1],Thickness[0.01]},  {RGBColor[1,0,0],Thickness[0.01]}}],
{n,0,5,1}]

Studiamo la serie di Fourier per x = π/2, punto in cui essa converge con somma π/2, e per x = π in cui la somma è ancora π/2. Per x = π/2 i termini della serie si scrivono (usando n al posto di k):

a[n]*Cos[n π/2]

Per n pari si annulla il coefficiente, per n dispari si annulla il coseno:

Simplify[a[2 k],k∈Integers]

Simplify[Cos[(2k+1) π/2],k∈Integers]

Gli altri termini sono

b[n]*Sin[n π/2]

che sono ovviamente nulli per n pari:

Simplify[%/.n->2k,k∈Integers]

mentre per n dispari diventano

Simplify[%%/.n->2k+1,k∈Integers]

In definitiva si ottiene un modo indiretto per sommare la serie di termine gennerale ((((-1))^k)/(1+2 k)), cosa che peraltro Mathematica sa fare direttamente:

(∑_(k=0)^∞)(((-1)^k)/(2k+1))

Per x = π i termini contenenti i seni si annullano, gli altri si scrivono

a[n] Cos[n π]

Simplify[%,n∈Integers]

Tali termini si annullano per n pari:

Simplify[%/.n->2k,k∈Integers]

mentre per n dispari si ha

Simplify[%%/.n->2k+1,k∈Integers]

Poiché per  x  = π la serie di Fourier converge al valore π/2, si ha un modo indiretto per sommare la serie di termine generale 1/(2k+1)^2, operazione che Mathematica peraltro sa fare:

(∑_(k=0)^∞)(1/(((2 k+1))^2))

Esempio 5: la funzione x^3 - x sull'intervallo [–1,1], 
prolungata con periodo 2

a=.;b=.

Commento: si tratta di una funzione derivabile con derivata prima continua su R; i coefficienti di Fourier tendono a 0  come (1/(k^3)).

Plot[(x^3)-x,{x,-1,1},AspectRatio->Automatic]

(∫_0^1)((t^3)-t) Sin[k π t]\[DifferentialD]t

Simplify[%, k ∈Integers]

b[k_]:=((12 (((-1))^k))/(((k π))^3))

Polinomi di Fourier:

s5[n_,x_]:=Sum[b[k] Sin[k Pi x],{k,1,n}]

Funzione che viene sviluppata:

f5[x_]:= If[Abs[x]<1,x^3-x, f5[Mod[x+1,2]-1]]

Alcuni grafici (monometrici):

Manipulate[
Plot[{f5[x], s5[n,x]},{x,-2,2},
AspectRatio->Automatic,
PlotRange->{{-2,2},{-1,1}},
PlotStyle->{{Thickness[0.012],RGBColor[0,0,1]},{RGBColor[1,0,0],Thickness[0.01]}}],
	{n,1,4,1}]

Tracciamo il grafico dell'errore (funzione data meno polinomio di Fourier del quarto ordine):

Plot[{(x^3)-x-s5[4,x]},{x,-1,1}]

Commento: attenzione alla scala sull'asse delle ordinate!
Spettro di ampiezza:

amp5 = Table[Abs[b[k]],{k,1,20}]//N;

diagBarre[amp5]

Commento: la convergenza a 0 dei coefficienti di Fourier è velocissima!

Esempio 6: impulso di ampiezza unitaria e durata τ su un periodo T ( τ < T ).

a=.;b=.

La funzione è pari: calcoliamo i coefficienti  (Subscript[a, k]):

a[k_]:=((4 (∫_0^(τ/2))Cos[((2 k π t)/T)]\[DifferentialD]t)/T)

a[0]/2

a[k]

Polinomi di Fourier:

s6[n_,t_,τ_,T_]:=τ/T + Sum[((2 Sin[k π τ/T])/(k π))  Cos[2 k π t/T],{k,1,n}]

Funzione sviluppata:

x0[t_,τ_] := If[Abs[t]≤τ/2, 1, 0]

x[t_,τ_,T_]:=If[Abs[t]<T/2, x0[t,τ], x[Mod[t+T/2,T]-T/2,τ,T]]

Alcuni grafici monometrici (con τ = 1, T = 2):

Manipulate[
Plot[{x[t,1,2], s6[n,t,1,2]},{t,-3,3},
AspectRatio->Automatic,
PlotRange->{-1,2},PlotStyle->{{Thickness[0.01],RGBColor[0,0,1]},{RGBColor[1,0,0],Thickness[0.01]}}],
{n,4,10,2}]

Altri grafici monometrici (con τ = 0.5, T = 2):

Manipulate[
Plot[{x[t,0.5,2], s6[n,t,0.5,2]},{t,-3,3},
AspectRatio->Automatic,
PlotRange->{-1,2},PlotStyle->{{Thickness[0.01],RGBColor[0,0,1]},{RGBColor[1,0,0],Thickness[0.01]}}],
{n,4,10,2}]

Se T è multiplo di τ, come negli esempi appena mostrati, cioè T = m τ con m intero, allora i coefficienti  a[k] sono nulli per k multiplo di m. Verifichiamolo sullo spettro di ampiezza:

amp6 =Table[ Abs[((2 Sin[k π τ/T])/(k π))], {k, 1, 12}];

amp6 = Join[{(τ/T)}, amp6]/.{T->2, τ->1};

diagBarre[amp6,0]

Esempio 7:  funzione seno raddrizzata

a=.; b=.;

a[k_]:=(((∫_0^π)Sin[t] Cos[k t]\[DifferentialD]t)/π)

b[k_]:=(((∫_0^π)Sin[t] Sin[k t]\[DifferentialD]t)/π)

a[k]

Simplify[%, k ∈Integers ∧ k ≠ 1]

a[0]/2

a[1]

b[1]

b[k]

Simplify[%, k ∈Integers ∧ k > 1]

s7[n_,t_]:=(1/π)+ Sin[t]/2 + Sum[-((1+(((-1))^k))/((-1+(k^2)) π))  Cos[k t],{k,2,n}]

Funzione sviluppata:

x[t_] := If[Sin[t] >0, Sin[t], 0]

Alcuni grafici (monometrici):

Manipulate[
Plot[{x[t], s7[n,t]},{t,-5,5},
AspectRatio->Automatic,
PlotRange->{-1,2},PlotStyle->{{Thickness[0.01],RGBColor[0,0,1]},{RGBColor[1,0,0],Thickness[0.01]}}],
{n,1,4,1}]

Spettro di ampiezza:

Table[Sqrt[a[k]^2+b[k]^2], {k,1,12}];

amp7 = Join[{a[0]/2}, %];

diagBarre[amp7,0]