Vedeni tepla na kruznici 

Rovnice vední tepla:  (∂_t) u[x] = Δ u[x]

1. DFT pořebuje jistou magii ohledně záporných frekvencí

ns = 32; (* počet vzorků *)
ff=2;  (* frekvence *)
Table[ {x,Sin[ ff x]},{x,0,2Pi,2Pi/ns}];
Show[
Plot[ Sin[ff x],{x,0, 2Pi},PlotTheme->"Detailed"],
ListPlot[%,PlotStyle->Red]
]

(* Fourierovská báze odpovídající DFT *)
fb[n_,x_] := Module[{n2=Quotient[n-1,2],odd=OddQ[n]},
basis = Table[Cos[k x]+I Sin[k x],{k,-n2,n2}];
If[ !odd,AppendTo[basis,0]];
RotateLeft[ basis,n2]
]

fb[7,x]

f0 = Exp[ - 4 Sin[s/2]^2]-4/5Exp[ - 12 Sin[s/2]^2];
Plot[f0,{s,-Pi,Pi}]

(*kolik bodů*)
nF=20;  (* MENTE ME *)
nFhalf=Quotient[nF-1,2]; (*max. frekvence, co se vejde, viz sampling thm.*)

(*nějaká analytická funkce*)
tstF=Function[s,Evaluate@f0];

(* periodicky samplovaný interval, tedy zahodím 2π, když už tam je 0, proto [[1;;-2]] *)
Slist = Range[0,2Pi,2.Pi/nF]    [[1;;-2]]

(*hodnoty F v požadovaných bodech*)
Flist=tstF /@ Slist  

(*spočtu DFT, pro větší nF (>~30) nutno vyhodit Chop*)
fourF=Fourier[Flist]//Chop

(*Inverzni Four. transf. má opačné znaménko, proto -s *)
fourSum = 1/Sqrt[nF] fb[nF,-s].fourF//Expand//Chop

Plot[ {tstF[s],fourSum},{s,-Pi,Pi},PlotLabel->"Porovnání obou fcí"]
Plot[ {tstF[s]-fourSum},{s,-Pi,Pi},PlotLabel->"Chyba"]
ListPlot[ Log10@Abs@({Re[#],Im[#]}&)@fourF[[1;;nFhalf]],PlotTheme->"Detailed",PlotLabel->"Spektr. konvergence"]

(* samozřejmě, takto spočtené koeficienty nejsou přesné, ale ... *)
2Pi/Sqrt[ nF ] fourF[[1]] - Integrate[ tstF[s] , {s,0,2Pi}]

fbEvol[n_,x_,t_] := Module[{n2=Quotient[n-1,2],odd=OddQ[n]},
basis = Table[Exp[-k^2 t]( Cos[k x]+I Sin[k x]),{k,-n2,n2}];
If[ !odd,AppendTo[basis,0]];
RotateLeft[ basis,n2]
]

fourSumEvol[t_Real] = 1/Sqrt[nF] fbEvol[nF,-s,t].fourF

Animate[
Rasterize@Plot[ fourSumEvol[t],{s,-Pi,Pi},PlotTheme->"Detailed",PlotLegends->None,PlotRange->{0,1/2},ImageSize->555],
{t,0,1}]

fsum[ak_] := Compile[ {{x,_Real}},
Sum[ak[[k+1]]Cos[k x],{k,0,Length[ak]-1}]
]

fbEvol2:=Compile[ {{n,_Integer},{x,_Real},{t,_Real}}, Module[{basis,n2=Quotient[n-1,2],odd=OddQ[n],k},
basis = Table[Exp[-k^2 t] (Cos[k x]+I Sin[k x]),{k,-n2,n2}];
If[ !odd,AppendTo[basis,0]];
RotateLeft[ basis,n2]
]

]

(*CompilePrint[fbEvol2] *)
 (* nebytne: ukazat priklad - kdyz zapomenu basis jako local!! *)

fbEvol3:=Compile[ {{n,_Integer},{x,_Real},{t,_Real}}, Module[{et,w,basis,n2=Quotient[n-1,2],odd=OddQ[n],k},
basis = Table[ 1.0 I, {k,-n2,n2}];
basis[[n2+1]]=1.0+0.0I;
w=Cos[x]+I Sin[x];
For [k=1,k<=n2,k++,
basis[[n2+k+1]]=basis[[n2+k]]*w;
basis[[n2-k+1]]=Conjugate[basis[[n2+k+1]]];
];
For [k=1,k<=n2,k++,
et=Exp[-k^2*t];
basis[[n2+k+1]]*=et;
basis[[n2-k+1]]*=et;
];
If[ !odd,AppendTo[basis,0]];
RotateLeft[ basis,n2]
]
]

fbEvol3[3,0.5,0.01]

Needs["CompiledFunctionTools`"]
CompilePrint[fbEvol3]

Timing[fbEvol[300,0.5,0.01];]
Timing[fbEvol2[300,0.5,0.01];]
Timing[fbEvol3[300,0.5,0.01];]

fourSumEvol3[t_Real] := 1/Sqrt[nF] fbEvol3[nF,-s,t].fourF;
Animate[
Rasterize@Plot[ fourSumEvol3[t],{s,-Pi,Pi},PlotTheme->"Detailed",PlotLegends->None,PlotRange->{0,1/2},ImageSize->555,PlotLabel->"t="<>ToString[t]],
{t,0,1,0.01}]

fbEvol2[3,0.5,0.1]
fbEvol3[3,0.5,0.1]