Spectral Solver for (Subscript[u, t])+L[u]+N[u]=0 on a 1D periodic domain

In the above L is a linear operator which can involve spatial derivatives etc. while N is a nonlinear operator which can also involve spatial derivatives.  The code is surprisingly easy to implement.  Lets look at a simple example:
	(Subscript[u, t])+a(x) (Subscript[u, x])=f(x,t)
we need to take trade in the functions u for a vector 𝕦∈(ℝ^n)of values sampled from the function u on a regular grid
	(Subscript[x, i])=a+i Δx
for i=0 to n-1 which  we immediately trade in for the Discrete fourier Transform (DFT) (OverTilde[𝕦]) of 𝕦 because it makes derivatives much simpler.   The transformed ea

Discrete Fourier Transforms

The Dscrete Fouier Transform turns a sampled signal from a function f[x] into a vector containing all the frequency components in the sample.  It is most naturally thought of as complex exponentials rather than real sines and cosines.   Different software packages and fields (math, physics, signal processing, etc.) number the components and scale differently.  Mathematica implelemts the most common varient by default in the Fourier function.  We are going to work out what the various parts mean.

If I sample a real function f[x] (imagined to be periodic) on an interval a<x<b with n equally spaced points then I want to sample at
	(Subscript[x, i])=a+i Δx
where Δx=(b-a)/n and i=0,1,…,(n-1)  or equivalently i=1,2,…,n-1,n and leave the periodic (or warp around) point off.  If this is a simple sine or cosine I can count the number of wiggles (which could be an integer) in the curve.

Manipulate[
{a,b}={0.0, 2π};
n = 123;
Δx=(b-a)/n;
f[x_]:= (Subscript[a, 0])+(Subscript[a, 1])Cos[(Subscript[m, 1]) x]+(Subscript[b, 1])Sin[(Subscript[m, 1]) x]+(Subscript[a, 2])Cos[(Subscript[m, 2]) x]+(Subscript[b, 2])Sin[(Subscript[m, 2]) x]+(Subscript[a, 3])Cos[(Subscript[m, 3]) x]+(Subscript[b, 3])Sin[(Subscript[m, 3]) x];
fS = Table[ f[ a + i Δx], {i, 0, n-1}];
GraphicsGrid[{{ListPlot[ fS,
PlotRange->{-10,10}], ListLogPlot[ Abs[Fourier[fS]],
PlotRange->{{0, 20},All}]}},
PlotLabel->"(Subscript[a, 0])+(∑_(i = 1)^3)(Subscript[a, i])Cos[(Subscript[m, i])x]+(Subscript[b, i])Sin[(Subscript[m, i])x]"],
{{(Subscript[a, 0]),1}, -3,3},
{{(Subscript[a, 1]),1},-3,3},
{{(Subscript[a, 2]),1}, -3,3},
{{(Subscript[a, 3]),1},-3,3},
{{(Subscript[b, 1]),0}, -3,3},
{{(Subscript[b, 2]),0}, -3,3},
{{(Subscript[b, 3]),0}, -3,3},
{{(Subscript[m, 1]),1}, -10,10},
{{(Subscript[m, 2]),2}, -10,10},
{{(Subscript[m, 3]),3}, -10,10}
]

Convolutions and Products

Suppose we multiply two functions together to get a new function
	f1[x]*f2[x]
we can sample it by just multiplying samples from f1 and f2 together.  The new sample inherits wiggles from both f1 and f2.

In fact, if we think about muliplying two sinusoids together we can use the product identity
	Sin[θ] Sin[ϕ] =(1/2)(Cos[θ-ϕ]-Cos[θ+ ϕ]) 
to identify exactly how the frequencies get mixed.  There are similar identities for other products of sines and cosines.   All the formula collapse when written as complex exponentials which is what the DFT does. Written in terms of the "wiggle counts" m1 and m2 that we are using on our reference interval this says 
	Sin[m1 x] Sin[m2 x] =(1/2)(Cos[(m1-m2)x]-Cos[(m1+m2)x]) 
in other words "wiggle counts" of m1 and m2 combine to give  m1-m2 and m1+ m2.   This is clear in a graph.

{m1,m2}={120,110}
Plot[ Sin[m1 x] Sin[m2 x], {x, 0, 2π}]

The question is if how do we combine DFTs of two signals to get the DFT of the product signal.  The answer is called a convolution that sums up all the possible ways to get the

{a,b}={0., 2π};
n= 124;
Δx=((b-a)/n);
x[i_]:= a+i *Δx
{f1[x_],f2[x_]}= {1.2 Sin[3x],3.2 Cos[x]-2.3Sin[2x]};
𝕗1=Table[ f1[x[i]],{i,0, n-1}];
𝕗2=Table[ f2[x[i]],{i,0, n-1}];
{ListPlot[
{𝕗1,𝕗2,𝕗1*𝕗2},
Joined->True,
Mesh->None],
ListLogPlot[
 Map[Abs[Fourier[#]]&,{𝕗1,𝕗2,𝕗1*𝕗2}] ,
PlotRange->{{1,10},All}]}

The question is if how do we combine DFTs of two signals to get the DFT of the product signal.  The answer is called a convolution which sums up all the possible ways to get a specific output frequency.  This is a very common operation and has its own function "ListConvolve" the only thing you need to know is that you need to tell it to wrap around and that the default version of Fourier needs a scaling of (Sqrt[n]).  You should not worry about the stuff down in the weeds of (10^(-15))-> (10^(-17)) scale on this plot those are floating point zeros.

a=Fourier[ 𝕗1 𝕗2];
b=ListConvolve[ Fourier[𝕗1] ,Fourier[𝕗2],{1,1}]/(Sqrt[n]);
ListLogPlot[ {Abs[a], Abs[b]},
PlotRange->{{0,20},All}]

In words: "The DFT of the product is the convolution of the DFTs.  Well, you need to scale by the square root of the dimension!

This also works the other way "The DFT of the convolution is the product of the DFTs.  Well, you still need to scale by the square root of the dimension!

a=Fourier[ListConvolve[ 𝕗1 ,𝕗2,{1,1}]]/(Sqrt[n]);
b= Fourier[𝕗1]*Fourier[𝕗2];
ListLogPlot[ {Abs[a],Abs[b]},
PlotRange->{{0,20},All}]

Now this signal (10^(-15)) lokks pretty small but if I look at the plot of the list convolve I can see that it is correct and now it is the (10^(-33)) that is the floating point zero stuff!

ListPlot[ ListConvolve[ 𝕗1 ,𝕗2,{1,1}]]