Fourier Series

Define a function that is periodic  with wavelength λ using the modulus function:

Boomerang[f_ ,  x_ ,  λ_ ]  :=  
      f[Mod[x, λ]]

AFunction[x_ ]  :=  ((3 - x)^3)/27

Plot[Boomerang[AFunction, x, 6], {x, -12, 12}, 
      PlotRange -> All]

Fourier Series

Demonstration of the orthogonality of the trigonometric functions

coscos  =  
      Integrate[
        Cos[((2 π Minteger x)/λ)] Cos[((2 π Ninteger x)/λ)] , {x, xo, xo  +  λ},  
        Assumptions -> {Minteger ϵ Integers,  
            Ninteger  ∈  Integers, 
            xo  ∈  Reals,  λ  ∈  
              Reals}]

Simplify[
      coscos /. {Minteger -> 4 ,  Ninteger -> 34}]

Assuming[
      Minteger ϵ Integers  &&  
        Ninteger  ∈  Integers  && 
        xo  ∈  Reals  &&  λ  ∈  Reals,  
      Limit[coscos, Minteger -> Ninteger]]

cossin  =  
      Integrate[
        Cos[((2 π Minteger x)/λ)] Sin[((2 π Ninteger x)/λ)] , {x, xo, xo  +  λ},  
        Assumptions -> {Minteger ϵ Integers,  
            Ninteger  ∈  Integers, 
            xo  ∈  Reals,  λ  ∈  
              Reals}]

Simplify[
      cossin /. {Minteger -> -7 ,  
          Ninteger -> 35}]

Assuming[
      Minteger ϵ Integers  &&  
        Ninteger  ∈  Integers  && 
        xo  ∈  Reals  &&  λ  ∈  Reals,  
      Limit[cossin, Minteger -> Ninteger]]

sinsin  =  
      Integrate[
        Sin[((2 π Minteger x)/λ)] Sin[((2 π Ninteger x)/λ)] , {x, xo, xo  +  λ},  
        Assumptions -> {Minteger ϵ Integers,  
            Ninteger  ∈  Integers, 
            xo  ∈  Reals,  λ  ∈  
              Reals}]

sinsin /. {Minteger -> 10 ,  Ninteger -> 9}

Assuming[
      Minteger ϵ Integers  &&  
        Ninteger  ∈  Integers  && 
        xo  ∈  Reals  &&  λ  ∈  Reals,  
      Limit[sinsin, Minteger -> Ninteger]]

Example of determination of fourier coefficients

Define functions based on the formulas derived for the fourier amplitudes

The constant term:

EvenTerms[0, function_ ,  
        wavelength_]  :=  (1/wavelength) (∫_0^wavelength) 
              function[dummy] \[DifferentialD]dummy

A function that defines each even amplitude individually (this is not very efficient, it would be better to evaluate the integral once and use that result)

EvenTerms[SP_Integer, function_ ,  
        wavelength_]  :=  EvenTerms[SP, function ,  
          wavelength] = (2/wavelength) (∫_0^wavelength) 
                function[
                dummy] Cos[((2 SP π dummy)/wavelength)] \[DifferentialD]dummy

Define the zeroeth odd term as zero for symmetry with the even terms:

OddTerms[0, function_ ,  wavelength_]  :=  0

OddTerms[SP_Integer, function_ ,  
        wavelength_]  :=  OddTerms[SP, function ,  
          wavelength] = (2/wavelength) (∫_0^wavelength) 
                function[
                dummy] Sin[((2 SP π dummy)/wavelength)] \[DifferentialD]dummy

A function to create a vector of amplitudes for the odd terms and one for the even terms

OddAmplitudeVector[NTerms_Integer,  function_,  
        wavelength_]  := 
Table[
        OddTerms[i,  function,  wavelength],  {i, 0, NTerms}]

EvenAmplitudeVector[NTerms_Integer,  function_,  
        wavelength_]  := 
Table[
        EvenTerms[i,  function,  wavelength],  {i, 0, NTerms}]

Define a function to use as an example:

myfunction[x_ ]  := (x*(2 - x)*(1 - x)^2)

OriginalPlot =  
      Plot[myfunction[x], {x, 0, 2}, 
        PlotStyle -> {Hue[1], Thickness[0.015]}]

EvenAmplitudeVector[6,  myfunction, 2]

The odd terms should all vanish because the function is even around the center x=(λ/2)

OddAmplitudeVector[6,  myfunction, 2]

OddBasisVector[NTerms_Integer, var_,   wavelength_]  :=  
      Table[Sin[((2 π  i var)/wavelength)],  {i, 0, 
          NTerms}]

OddBasisVector[6, x, 2]

EvenBasisVector[NTerms_Integer, var_,   wavelength_]  :=  
      Table[Cos[((2 π  i var)/wavelength)],  {i, 0, 
          NTerms}]

EvenBasisVector[6, x, 2]

FourierTruncSeries[n_,  function_, var_ ,  wavelength_] :=  
      EvenAmplitudeVector[n, function, wavelength] . 
          EvenBasisVector[n, var, 
            wavelength]  +  
OddAmplitudeVector[n, 
            function, wavelength] . 
          OddBasisVector[n, var, wavelength]

FourierTruncSeries[6, myfunction, x, 2]

FourierPlot  = 
      Plot[FourierTruncSeries[3, myfunction, x, 2], {x, -2, 
          4}]

Show[OriginalPlot, FourierPlot]

FourierPlot  = 
      Plot[FourierTruncSeries[6, myfunction, x, 2], {x, -2, 
          4}]

Show[OriginalPlot, FourierPlot]

Fourier Series Coefficients using Mathematica's package:

<< Calculus`FourierTransform`

AFunction[x_]  :=  (x - 3)^(3/27)

Plot[AFunction[x], {x, 0, 6}]

Mathematica's Fourier Series functions are defined for function that are periodic in the domain x ∈ (-1/2,1/2).  So we need to map the perioidic functions to this domain

ReduceHalfHalf[f_ ,  x_ ,  λ_ ]  :=  
      f[(x + 1/2)*λ ]

ReducedFunction  =  
      ReduceHalfHalf[AFunction, x, 6] // Simplify

Plot[ReducedFunction, {x, -1/2, 1/2}, 
      PlotRange -> All]

FourierCosCoefficient[ReducedFunction, x, n]

FourierSinCoefficient[ReducedFunction, x, n]

FourierTrigSeries[ReducedFunction, x, 5]

AnimateTruncatedFourierSeriesVersionOne[
          function_,  {truncationstart_,  truncationend_,  
            truncjump_}]  :=  
        Table[Plot[
            Evaluate[FourierTrigSeries[function, x, itrunc]], {x, -1, 1}, 
            PlotRange ->  {-2, 2}], {itrunc, truncationstart, 
            truncationend,  truncjump}];

The function above will work, but it is horribly inefficient!  Because it asks FourierTrigSeries to calculate one more term each time, it is doing some redundant work.  We can fix this up by having
 it calculate one new term and adding to the sum calculated previously. Here it is:

AnimateTruncatedFourierSeriesVersionTwo[
        function_,  {truncationstart_,  truncationend_}]  :=  
      Module[{coscof, sincof,  currentappx , n,  
          TwoPi =  2 π},  
currentappx  =  
          FourierCosCoefficient[function, x, 0]; 
coscof[
            n_]  =  
          Simplify[
            FourierCosCoefficient[function, x, n]]; 
sincof[
            n_]  =  
          Simplify[
            FourierSinCoefficient[function, x, 
              n]]; 

Table[
          Plot[Evaluate[
              currentappx  +=  
                coscof[itrunc]*Cos[TwoPi itrunc x]  +   sincof[itrunc]*
                    Sin[TwoPi itrunc  x]], {x, -1, 1}, 
            PlotRange ->  {-2, 2}], {itrunc, truncationstart, 
            truncationend}];]

The following will demonstrate how convergence is difficult where the function changes rapidly---this is known as Gibbs Phenomenon

AnimateTruncatedFourierSeriesVersionTwo[
      ReducedFunction, {1, 60}]

Recall the definition used above for myfunction

myfunction[x]

Reduce myfunction to a domain so that we can use Mathematica's Fourier Package on it.

ReducedMyFunction  =  
      ReduceHalfHalf[myfunction, x, 2] // Simplify

General form of the even amplitudes

myfunccos  = 
      FourierCosCoefficient[4 (x^2) - 16 (x^4), x, n]

FourierSinCoefficient[4 (x^2) - 16 (x^4), x, n]

Illustrate how amplitude of each term behaves. This plot indicates that the short wavelength contribtuions are smaller and smaller.

ListPlot[Table[myfunccos, {n, 1, 50}]]

ListPlot[Table[myfunccos, {n, 1, 10}], PlotJoined -> True, 
      PlotRange -> All]