pdefourier: A package for doing Fourier analysis and solving partial differential equations in Maxima CAS

Fourier analysis provides a set of techniques for solving partial differential equations (PDEs) in both bounded and unbounded domains, and various types of initial conditions. In the bounded domain case, the method of separation of variables leads to a well-defined algorithm for developing the solution in a Fourier series, making this problem tractable with a CAS.

Table of Contents

Overview
Installation
Fourier coefficients and series
Frequency spectrum
The heat equation
The wave equation
The Laplace equation
Bessel functions

<a name="overview">Overview<a/>

This Maxima package computes symbolically the Fourier of piecewise-smooth functions. Using the method of separation of variables it is also able to symbolically solve the one-dimensional heat and wave equations on a domain [0,L], with regular Sturm-Liouville conditions, that is, general boundary conditions of the form:

<p align="left"> &alpha;₁u(0,t) + &beta;₁u_x(0,t) = h₁(t) <br> &alpha;₂u(L,t) + &beta;₂u_x(L,t) = h₂(t) </p>

Moreover, the package can solve the two-dimensional Laplace equation for a variety of domains (rectangles, disks, annuli, wedges) either with Dirichlet or Neumann boundary conditions. In the case of a rectangular domain [0,a] x [0,b], the package can solve the Laplace equation with mixed boundary conditions of the form

<p align="left"> (1-&alpha;)u(x,0) + &alpha;u_y(x,0) = f₀(x), 0 &le; x &le; a<br> (1-&beta;)u(x,b) + &beta;u_y(x,b) = f_b(x), 0 &le; x &le; a<br> (1-&gamma;)u(0,y) + &gamma;u_x(0,y) = g₀(y), 0 &le; y &le; b<br> (1-&delta;)u(a,y) + &delta;u_x(a,y) = g_a(y), 0 &le; y &le; b<br> </p>

where α,β,γ,δ ∈ {0,1}.

Of course, in all cases it is possible to truncate a series to make numerical calculations.

Remark: pdefourier automatically loads other packages, including piecewise, draw, simplify_sum, and syntactic_factor. The Maxima package syntactic_factor was written by Stavros Macrackis. Currently, pdefourier loads the files piecewise and syntactic_factor contained in this repository.

The Documentation folder folder contains a pdf file explaining some technical details of the implementation and a description of many of the functions contained in the package, as well as their syntax. Also, there is a Maxima session (in wxm format) with lots of examples, graphics, animations and tips for use. Here, we only give a quick introduction to the main commands used for solving typical problems.

<a name="installation">Installation<a/>

The package can be installed by putting a copy of the files pdefourier.mac, syntactic_factor.mac, and piecewise.mac inside a folder contained in the environment variable file_search_maxima. In a Linux box, such a system-wide location could be something like /usr/share/maxima/5.42.2/share/contrib/, while in a Windows environment typically it will be c:\Program files (x86)\Maxima-sbcl-5.42.0/share/maxima/5.42.0/share/contrib (you may need administrator rights in order to do that in either case). The package can then be loaded with the command load(pdefourier) inside a Maxima session.

<a name="foucoeff">Fourier coefficients and series<a/>

The package can deal with piecewise functions, defined in natural notation:

<p align="left"> <code>(%i1) load("pdefourier.mac")$ </code><br> <code>(%i2) v(x):=if (-%pi<=x and x<0) then x^2 elseif (0<=x and x<=%pi) then sin(3*x)$</code> </p>

It is possible to detect the parity of such a functions, with paritycheck; possible outcomes are even, odd or none:

<p align="left"> <code>(%i3) paritycheck(v(x),x);</code><br> <code>(%i2) none</code> </p>

Let us draw the curve to chek the answer:

<p align="left"> <code>(%i4) plot2d(v(x),[x,-%pi,%pi],[ylabel,"v(x)"]);</code><br> <code>(%t4)</code> </p> <img src="img/Example-01.png">

The Fourier coefficients are computed with fouriercoeff, whose syntax is

<p align="center"> <code>fouriercoeff(expr,var,p)</code> </p>

Here p=(b-a)/2 if the whole interval of definition for expr is [a,b]. The output has the form

<p align="center"> <code>[[a₀,a_n,b_n],svlist]</code> </p>

where svlist is a sublist containing the singular values of the coefficients, again with the format

<p align="center"> <code>[m,a_m,b_m]</code> </p>

In the example we are considering, notice that the function is defined on [-ππ]:

<p align="left"> <code>(%i5) fouriercoeff(v(x),x,%pi);</code><br> <img src="img/fcoeff_vx.png">

This example illustrates the presence of singular values of the coefficients (for n=3). We can approximate the function by its Fourier series truncated to order 15:

<p align="left"> <code>(%i6) vseries15:fourier_series(v(x),x,%pi,15)$</code><br> <code>(%i7) wxplot2d([v(x),vseries15],[x,-%pi,%pi],[legend,false]);</code><br> <code>(%t7) </code> </p> <img src="img/Example-02.png">

Here is a well-known example of an unbounded function:

<p align="left"> <code>(%i8) absolute(x):=if (x<=0) then -x elseif (x>0) then x$</code><br> <code>(%i9) paritycheck(absolute(x),x);</code><br> <code>(%o9) even</code> </p>

and its bounded version, for which we compute the Fourier series:

<p align="left"> <code>(%i10) absolute0(x):=if ( x>=-1 and x<=0) then -x elseif (x>0 and x<=1) then x$</code><br> <code>(%i11) paritycheck(absolute0(x),x);</code><br> <code>(%o11) even</code><br> <code>(%i12) fourier_series(absolute0(x),x,1,inf);</code><br> </p> <img src="img/abs_series.png">

<a name="frequency">Frequency spectrum<a/>

Frequency analysis is very useful in Engineering applications (but also in Physics). This technique requires that the Fourier series be first re-expressed using the identity

<p align="left"> a cos(w)+b sin(w)=r cos(w-u) </p>

where the modulus and the phase shift are given, respectively, by r=sqrt(a²+b²) and u=atan(b/a). Then, we can rewrite the terms of the series in the so-called harmonic form:

<p align="left"> c_n cos(nwx-u_n) </p>

In the theory of sound, the first harmonic (corresponding to n=1) is called the fundamental harmonic. The remaining ones are called overtones. The coefficients c_n are the harmonic amplitudes, and the absolute value |c_n| is a measure of the relative importance of the nth harmonic in a given signal (sound). The frequency analysis is done in pdefourier with the aid of the function fourier_freq. The function fourier_harm returns a list with the first n harmonics of a given function, with the syntax

<p align="center"> <code>fourier_harm(function(variable),variable,p,n)</code> </p>

where p=(b-a)/2 and [a,b] is the interval of definition of the function, while fourier_freq_list (with the same syntax) gives the list of the amplitudes.

Here is the example of the square pulse with compact support on [-2,2]:

<p align="left"> <code>(%i1) load(pdefourier)$</code><br> <code>(%i2) square0(x):=if (-2<=x and x<-1) then 0 elseif (-1<=x and x<=1) then 1 elseif (1 < x and x< =2) then 0$</code><br> <code>(%i3) plot2d(square0(x),[x,-2,2],[y,-0.1,1.1],[ylabel,"square pulse"]);</code><br> <code>(%t3) </code> </p> <img src="img/Example-harm00.png"> <p align="left"> <code>(%i4) fourier_harm(square0(x),x,2,5);</code><br> </p> <img src="img/fourier_harm_square00.png"> <p> <code>(%i5) fourier_freq_list(square0(x),x,2,5);</code><br> <code>(%o5) [[1.0,0.6366197723675814],[2.0,0.0],[3.0,0.2122065907891938],[4.0,0.0],[5.0,0.1273239544735163]]</code><br> </p>

For smooth or mildly non-smooth functions, the Fourier series converges very fast, and so the first few terms suffice to get a good approximation. This is reflected in the relative weight of each harmonic (fastly decreasing), and can be visualized (using the wxMaxima frontend) with the command wxfourier_freq:

<p align="left"> <code>(%i6) wxfourier_freq(square0(x),x,2,10);</code><br> <code>(%t6) </code> </p> <img src="img/Example-harm01.png">

For very discontinous functions, the amplitude of harmonics does not decrease that fast (or does not decrease at all). Consider the following function as an example:

<p align="left"> <code>(%i7) f0(x):=if (x>=-4 and x<-3 ) then 0 elseif (-3<=x and x<=-2) then 1 elseif (-2<x and x<-1) then 0 elseif (-1<=x and x<=1) then -x elseif (1<=x and x<=2) then 0 elseif (2<x and x<3) then -1 elseif (x>=3 and x<=4) then 0$</code> </p>

Here is its graphical representation, along with its Fourier approximation up to order 15:

<p align="left"> <code>(%i8) seriesf0:fourier_series(f0(x),x,4,15)$</code><br> <code>(%i9) plot2d([f0(x),seriesf0],[x,-4,4],[legend,false]);</code><br> <code>(%t9) </code> </p> <img src="img/Example-harm02.png">

And its frequency spectrum (containing the first 10 harmonics):

<p align="left"> <code>(%i10) wxfourier_freq(f0(x),x,4,10);</code><br> <code>(%t10) </code> </p> <img src="img/Example-harm03.png">

Finally, notice that adding terms to the Fourier series reflects in the harmonic content of the signal, as the following animation shows (this too requires the wxMaxima frontend):

<p align="left"> <code>(%i11) ramp(x):=if (-5<=x and x<=-1) then (x+3)/2 else 0$</code><br> <code>(%i12) define(ramp_series(x,n),fourier_series(ramp(x),x,2,n))$</code> <code><pre>(%i13) with_slider_draw(k,makelist(j,j,1,10), dimensions=[900,450], xrange=[-5.25,10.25], yrange_secondary=[-1.45,1.45], axis_top