Symbolic computation of Fourier series 

Wilhelm Werner
Hochschule Heilbronn, Germany
werner@hs-heilbronn.de
2006 

Introduction:  

This worksheet demonstrates the use of Maple for the symbolic computation of Fourier expansions.  

This is an updated version of a package originally published in the Maple Application Center (2000). 

 

Taylor series expansion is provided by Maple via the taylor or series command; there is no corresponding facility for Fourier expansions  

 

however. The module FourierSeries designed for use with Maple 6 (or later) fills this gap. Real as well as complex Fourier coefficients are computed symbolically; furthermore some classroom visualizations are provided. 

We do not give any applications of Fourier series here;  for classroom introductions to this topic we refer to other contributions within the Maple Application Center. Instead we concentrate on the algorithms required to compute those expansions symbolically. 

>

restart;

 

Acknowledgment 

 

The author is indepted to Prof.Dr.Robert Lopez, Maplesoft, for several valuable comments which led to improvements of this package. 

 

Section I: What makes symbolic computation of Fourier series cumbersome? 

There are a lot of calculus texts which deal with the computation of Fourier coefficients; the relevant formulas e.g. on the interval [0,T] are well known: 

 

 

The fact that the computation of these coefficients "only" requires integration seems to indicate that these computations are an easy task for a computer algebra system like Maple. This, however, is not the case as is demonstrated with the following simple example. Let us consider the function 

>	f:=t->t^2+cos(4*t);

 

 

A straightforward implementation of the above formulas on the interval [0,]  

>	T:=2*Pi:k:='k':

 

>	a:=unapply(2/T*int(f(t)*cos(2*Pi/T*k*t),t=0..T),k);

 

 

either results in an error message 

>	a(3);seq(a(k),k=0..6);

 

 

Error, (in a) numeric exception: division by zero 

or in coefficients which only cover the "generic case": 

>	assume(k,integer); simplify(a(k));

 

 

This result is wrong for k=0 and k=4! Assuming k to be integer forces Maple to do some a priori simplifications like 

>	sin(k*Pi); cos(Pi*k);

 

 

 

in the numerator which cause the problem. Probably many users encountered these problems and circumvented them by replacing symbolic integration with numerical integration methods. For graphical purposes this approach is acceptable, it is unsatisfactory in the context of computer algebra though.  

Basically the approach which is taken by FourierSeries mimics hand calculation; the routine computes the generic solution and thereafter looks for critical cases. These exceptional cases are then handled separately.  

The major problem of this approach is that the number of exceptional cases as well as their location is not known in advance. FourierSeries uses Maple's solve command to find exceptional cases and Maple's int command for integration;  the limitations of the package therefore are basically determined by the limitations of these Maple procedures, which, however, are known to be extremely powerful. There is no internal maximum number of exceptional cases within FourierSeries;  any example with an infinite set of exceptional points would result in a failure of FourierSeries, however.  

Section II: Installation issues 

  FourierSeries is provided as a Maple module and therefore requires at least Maple 6 to work (for a version of the software for Maple V Rel.4 or Rel.5 contact the author); the software has been tested within the Windows 9x environment only so no claims whatsoever can be made for other platforms like Linux or Macintosh.  

The easiest way to include FourierSeries into your Maple environment is to generate a directory/folder like 

"C:\programs\Maple 6\UserLibs\FourierLib" 

and copy all files into that directory. Within Maple change the libname according to 

>

restart:

 

>	libname:="C://Programme//Maple 9.6//Users//FourierLib",libname;

 

 

or change your Maple ini file accordingly.  If you encounter an error message like 

>	with(FourierSeries);

 

Error, (in with) invalid input: with expects its 1st argument, pname, to be of type {package, module}, but received FourierSeries 

then most probably you provided an erroneous pathname within libname (check for "/" and "\"!) or you forgot to extend  the libname variable. 

>

restart:

 

>	libname:="C:\\Programme//Maple 9.5//Users//FourierLib2",libname;

 

 

>	with(FourierSeries);

 

 

You get help information in the usual manner: 

>	?FourierCoeff

 

One can call the routines of the module also directly without the with(FourierSeries); statement in the form 

>	FourierSeries:-FourierPlot([t->t,-1..1],6,-2..2,title="Fourier plot");

 

 

as well. 

Section III: Examples 

The following examples demonstrate the various options of  FourierSeries: 

Computation of real Fourier coefficients 

The heart of  FourierSeries is the procedure FourierCoeff for the computation of real and complex Fourier coefficients; FourierCoeff is called by all other routines of the package. 

>	f:=t->t^2+cos(4*t); intv:=-Pi..Pi;

 

 

 

>	FourierSeries:-FourierCoeff([f,intv],"FourierReal"):

 

 

Alternatively you can use the FourierSin resp. FourierCos option: 

>	FourierSeries:-FourierCoeff([f,intv],"FourierCos"):

 

 

In case of a piecewise linear function like 

>	f:=t->piecewise(t<-1,0,t<-1/2,1+t,t<1/2,1/2,t<1,1-t,0);

 

 

>	plot(f,-1..1);

 

 

FourierCoeff supports a simplified input mechanism: The user has to provide a list of points only which describe the polygon: 

>	FourierSeries:-FourierCoeff([[[-1,0],[-1/2,1/2],[1/2,1/2],[1,0]]],"FourierCos"):

 

 

The global variable _FourierSeriesShowFunction controls the display of the underlying function; this option is especially useful, if the user had provided a list of points only. (Note: The global variable _FourierSeriesShowPiecewise used in previous versions of the package is no longer supported.) 

>	_FourierSeriesShowFunction:=true: # default is false

 

>	FourierSeries:-FourierCoeff([[[-1,0],[-1/2,1/2],[1/2,1/2],[1,0]]],"FourierCos"):

 

 

 

Computation of complex Fourier coefficients 

Internally FourierSeries computes the complex Fourier expansion: 

 

The coefficients c[k] are related to a[k] resp. b[k] via the equations 

>	c[k]=(a[k]-I*b[k])/2; # if k>=0, b[0]:=0

 

 

>	c[-k]:=conjugate(c[k]);

 

 

The complex Fourier coefficients can be computed with option FourierExp or FourierComplex: 

>	f:=t->t+t^2+cos(4*t): intv:=-Pi..Pi:

 

>	FourierSeries:-FourierCoeff([f,intv],"FourierExp"): # or "FourierComplex"

 

 

 

>

 

 

>

 

 

>

 

 

>

 

 

Note that FourierCoeff returns a function which may be used to compute specific values of the Fourier coefficients: 

>	a(0), a(4),a(5),b(0), b(4),b(5),c(0), c(4),c(5);

 

 

Even or odd extensions of real functions defined on the positive resp. negative axis can be generated by 

>

 

 

>

 

 

Real coefficients may be derived thereof via the above mentioned relations easily. It should be mentioned that 

FourierCoeff also accepts parameters in the definition of functions or intervals: 

>	FourierSeries:-FourierCoeff([t->sin(alpha*t),-Pi..Pi],"FourierComplex"):

 

 

 

In that case FourierCoeff makes suitable assumptions (e.g. parameters a real) for those parameters; the above formula obviously is not correct for all n if is an integer for example. It is up to the user to check which assumptions must be fulfilled in these cases: 

>	FourierSeries:-FourierCoeff([t->sin(alpha+t),-beta..beta],"FourierComplex"):

 

 

 

The limitations of this feature is basically determined by Maple's int command; e.g.  

>	k:='k':assume(k,integer):int(bernoulli(2*k,x),x=0..1);

 

 

In case Maple can't compute the integral symbolically FourierSeries also fails; of course any concrete instance is handled correctly: 

>	FourierSeries:-FourierCoeff([t->bernoulli(6,t),0..1],"FourierReal"):

 

 

>	FourierSeries:-FourierCoeff([t->bernoulli(7,t),0..1],"FourierReal"):

 

 

Graphic output 

Basically for educational purposes FourierSeries provides three formats for graphical output: 

>	FourierSeries:-FourierPlot([t->t,-Pi..Pi],6,-2*Pi..2*Pi,title="Fourier approximation of order 6",color=[red,blue],thickness=[2,1]);

 

 

FourierPlot as well as FourierAnimate accept any plot options which apply to Maple's plot resp. animate command, e.g. title, legend, color etc. 

>	FourierSeries:-FourierAnimate([t->t,-Pi..Pi],12,-2*Pi..2*Pi,

 

>	title="Fourier approximation up to order 12",                          color=[red,blue],thickness=[2,1],discont=true);

 

 

The spectrum of a function is displayed by FourierSpectrum: 

>	FourierSeries:-FourierSpectrum([t->t,-Pi..Pi],20,axes=normal,title="Spectrum");

 

 

FourierSpectrum uses Maple's histogram command from the stats package; thus the plot options are restricted to those accepted by histogram. 

Note: No graphics is generated if the function definition and/or interval contains additional parameters. 

>	FourierSeries:-FourierSpectrum([t->a*t,-Pi..Pi],20,axes=normal,title="Spectrum");

 

Error, (in FourierSeries:-FourierSpectrum) unspecified parameter(s), f 

The plot procedures FourierPlot, FourierAnimate and FourierSpectrum support the simplified input mechanism for piecewise linear functions, too: 

>

 

 

>	FourierSeries:-FourierPlot([[[-1,0],[-1/2,1/2],[1/2,1/2],[1,0]]],6,-2..2,color=[red,blue],axes=normal);

 

 

 

>

 

 

>	FourierSeries:-FourierAnimate([[[-1,0],[-1/2,1/2],[1/2,1/2],[1,0]]],6,-2..2,color=[red,blue],axes=normal);

 

 

>	FourierSeries:-FourierSpectrum([[[-1,0],[-1/2,1/2],[1/2,1/2],[1,0]]],10,axes=normal);

 

 

Application 

 

It is well known, that global smoothness properties of the periodic continuation of a function determine how fast Fourier coefficients tend to zero; this can be demonstrated easily as follows: 

1. discontinuity ("jump"); coefficients are O(): 

>	FourierSeries:-FourierCoeff([t->t,-1..1],"FourierComplex"):

 

 

2. continous; discontinuous first derivative; the coefficients are O(): 

>	FourierSeries:-FourierCoeff([t->t^2,-1..1],"FourierComplex"):

 

 

3. first, second derivative continous, third derivative has a jump discontinuity; therefore the coefficients are O(): 

>	x0:=1/sqrt(3):

 

>	f:=t->t^2*(t^2-1)^2;

 

>	plot(f,-x0..x0);

 

>	map(expand,Limit(D(f)(x),x=x0,left)=Limit(D(f)(x),x=-x0,right)); map(value,%);

 

>	map(expand,Limit((D@@2)(f)(x),x=x0,left)=Limit((D@@2)(f)(x),x=-x0,right)); map(value,%);

 

>	map(expand,Limit((D@@3)(f)(x),x=x0,left)=Limit((D@@3)(f)(x),x=-x0,right)); map(value,%);

 

 

 

 

 

 

 

 

 

>	c:=FourierSeries:-FourierCoeff([f,-x0..x0],"FourierComplex"):

 

 

>	FourierSeries:-FourierPlot([f,-x0..x0],2,-2*x0..2*x0,title="Fourier approximation",color=[red,blue]);

 

 

>	series(c(n),n=infinity,3); # asymptotic behaviour of coefficients

 

 

References 

 

FourierSeries is one of the routines developped by the author for the (german) book 

W.Werner, Mathematik lernen mit Maple, Band 2, dpunkt.verlag, Heidelberg 1998,  

    ISBN 3-920993-94-2 (633 p., CD included) (out of print) 

 

 

Conclusion:  

FourierSeries extends Maple's capabilities to compute symbolically the coefficients of Fourier series..  

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.