Automated solving of constant-coefficients second-order linear PDEs using Fourier analysis

04/02/2020 ∙ by Emmanuel Roque, et al. ∙ Cinvestav 0

After discussing the limitations of current commercial software, we provide the details of an implementation of Fourier techniques for solving second-order linear partial differential equations (with constant coefficients) using a computer algebra system. The general Sturm-Liouville problem for the heat, wave and Laplace operators on the most common bounded domains is covered, as well as the general second-order linear parabolic equation with constant coefficients, which includes cases such as the convection-diffusion equation by reduction to the heat equation.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The need for explicit solutions of differential equations (be them ordinary or partial) arises in the everyday practice of physicists, engineers and scientists in general. The class of ordinary differential equations that can be solved by elementary methods, mainly clever changes of variables and linear algebra, is quite ample and, accordingly, there exists a vast amount of software devoted to finding both their explicit and numerical solutions. This is in contrast with the case of partial differential equations (PDEs). Here, solution techniques are more advanced, involving Fourier analysis for bounded domains and integral transforms in the unbounded case

111In this paper we restrict our attention to classical solutions, thus leaving aside considerations regarding weak solutions and the use of such tools as distribution theory. and, although there is a good deal of software packages related to numerical computations (most of them using the finite element method), there is a shortage of options for finding explicit analytic solutions. In fact, to the best of our knowledge, only Maplesoft’s Maple™ and Wolfram’s Mathematica™ –both closed and proprietary software– are capable of this task. Both are excellent general purpose computer algebra systems (CAS) and do their job fast and efficiently, so one could wonder what is the motivation for considering the problem of automated solution of PDEs or thinking it is not already solved. To this question we can offer the following answers:

  1. Although they are very good at finding solutions, still there are examples where Maple™ or Mathematica™ can not find it (see examples in subsection 4.2). In those cases, it is very difficult to determine the origin of the failure, as both are closed source packages. The user can have a hard time trying to determine whether there is a problem in the formulation of the problem or in the method of solution.

  2. Proprietary software usually has a cost, and this case is no exception. While both packages are worth their price, this can be a serious hurdle for people working in underdeveloped countries, students in general or even freelancers for whom an open source alternative would be desirable.

  3. From a pedagogical point of view, the use of closed source software is like using a black box for finding answers. The teacher can expose an algorithm in the classroom and then say it is exemplified by the output, but really there is no clue about what is being actually computed, if there is any mathematical trick used which is not contemplated in the algorithm, or if use is made of pre-computed tables of particular cases, and so on.

  4. Lastly, the problem is interesting in itself from a mathematical and computational perspective, due to the blending of purely symbolic and numeric algorithms it requires, as will be seen in detail below.

These reasons led us to write a library, called pdefourier, for solving the general Sturm-Liouville problem for the heat, wave and Laplace operators with sources, as well as some particular problems reducible to these, such as the convection-diffusion equation. Our particular implementation uses the open source Maxima CAS [7], which has a programming language very closely tied to LISP and should be translatable to any other CAS or programming language in general, such as Python or Symbolic C++, without difficulty.

In Section 2 we comment on the challenges posed by the automated computation of solutions using Fourier techniques in general. By including examples of computations using Mathematica™ , we present a panorama of the current state of the art222We must insist in this aspect, since nothing could be further from our minds than to do a comparison test or a benchmark against a commercial software, supported by a whole team of scientists..


In what follows, all Mathematica™ commands have been executed using version on a Linux platform. The Maxima commands have been executed on the same Linux machine, with version 5.42.2 compiled against the LISP implementation SBCL 1.4.4. It is assumed that the package pdefourier has already been loaded. For installation instructions see Appendix A.

Section 3 describes how we handle the issues mentioned in Section 2, detailing some of the core features of the library, and we finish in Section 4 by showing some examples and applications.

2 Sturm-Liouville problems and eigenfunctions expansion

Notice that, from a mathematical point of view, the basic models given by the heat, wave and Laplace equations are all second-order linear partial differential equations with constant coefficients. They are supplemented by some initial and/or boundary conditions (the general situation contemplates both, and it is known as an IBVP or initial-boundary value problem), and in this work we will assume that the problem is set on a bounded spatial domain , so these conditions will be given at the boundary points (the only exception being Laplace’s equation, which will be considered on a plane domain). The general problem for the heat and wave equations can be written as:


where denotes either the heat or the wave operator. In the case of the Laplace’s equation we consider a wider class of plane domains (rectangles, disks, wedges, annuli) and boundary conditions. In a rectangular domain we consider the problem with mixed boundary conditions of the form:


On the remaining domains, the Laplace equation is solved in polar coordinates with Dirichlet and Neumann boundary conditions; the only exception being the case of an annular domain, where only Dirichlet boundary conditions are considered. Here is a summary of the cases:


A common procedure to solve these equations is to apply separation of variables: Writing

and denoting generically or by , leads to a regular Sturm-Liouville problem whose general form is


where are positive real functions, and are continuous on . The abstract theory of Sturm-Liouville operators is very well developed (see [1, 3]), and its main result can be summarized as follows:

Theorem 2.1.

Let be a Sturm-Liouville operator as in (4), subject to the boundary conditions derived from those in (1) in the homogeneous case. Assume that, for ,


Then, the following hold:

  1. All the eigenvalues

    are real and form a denumerable set.

  2. The set of all eigenfunctions,

    , is an orthogonal set.

  3. The set is complete in the space of square-integrable functions .

  4. The sequence diverges to and is bounded from below.

  5. If are all non-negative, then the eigenvalues are all non-negative.

Notice that the preceding theorem applies to homogeneous initial and boundary conditions and . To deal with the non-homogeneous case , we use the method of eigenfunction expansion or finite cosine or sine transform (see [6, 9]). Once we know how to deal with the general Sturm-Liouville problem for the heat, wave and Laplace operators, we can solve some other problems that can be reduced to one of these. The cases of the convection-difussion equation and heat loss through lateral boundaries are then particular cases. Thus, consider the general expression of a second-order linear parabolic equation with constant coefficients with sources:


Let us apply the change of variable


Substituting into (6) we get

hence, the original equation in terms of the new dependent variable reduces to the heat equation with sources

Of course, the initial and boundary conditions for must be transformed accordingly into conditions for .

The Sturm-Liouville theory is the basic tool used in our Maxima package pdefourier. The algorithm can be summarized as follows: we first solve the eigenvalue problem associated to (4), normalize the corresponding eigenfunctions and then proceed to develop the solution to (1) in term of eigenfunctions, taking non-homogeneous conditions into account with the use of the methods commented above. These steps involve some delicate issues that we discuss next.

Singular values of the Fourier coefficients

. The computation of Fourier coefficients in a CAS can be done using integration formulas. Particular care must be taken to separate singular values that lead to divergences in the general expressions of the Fourier coefficients. There are two types of special input functions whose Fourier coefficients must be calculated separately, namely, the functions belonging to one of the cosine, sine, trigonometric or complex orthogonal systems or a product of a polynomial with them. In these cases, the Fourier coefficients have a general expression

with a polynomial, which is only valid for the non-singular values of . The singular values must be treated separately. This process is, of course, transparent to the user, but lies at the bottom of the difficulties experienced by such mature software as Mathematica™ . Due to the widespread of this CAS, we will use it to get a benchmark for our library performance.

To get a sense of the kind of problems we are talking about, recall the definition of the Fourier coefficients of a periodic function (we consider the case for simplicity):

where is a natural number. In practice, the explicit computations of these coefficients depends heavily on the orthogonality relations ()

However, most computer algebra systems fail to recognize the possibility that . Thus, Mathematica™ gives

In[1]:= Integrate[Cos[m x] Cos[n x], {x, -Pi, Pi}]
Out[1]= (2 m Cos[n \[Pi]] Sin[m \[Pi]] - 2 n Cos[m \[Pi]] Sin[n \[Pi]])/(m^2 - n^2)

and the same is true of Maxima:


This fact has many consequences due to the presence of singular values of , as pointed above. Take as an example the computation of Fourier coefficients for


Using the trigonometric identity

and integrating by parts, it is straightforward to obtain

Mathematica™ fails to compute correctly the Fourier coefficients directly from the integral formulas because it does not recognize the singularity corresponding to :

In[2]:= Assuming[Element[n,Integers],1/Pi Integrate[Cos[n x] 3 x^2 Cos[7 x],{x,-Pi,Pi}]]
Out[2]= -((12 (-1)^n (49 + n^2))/(-49 + n^2)^2)

and so does Maxima:


Equivalence of trigonometric expressions

. Another source of concern when dealing with trigonometric expressions is the existence of multiple, equivalent ways of writing them. This is a well-known issue in different CASs. The set of heuristics rules Mathematica™  appears to use inside its Fourier coefficients built-in functions try to support common special cases, but are highly sensitive to the way the input is written and do not take advantage of the linearity of the integral nor distribute over subintervals of piecewise-defined functions. For instance, using equation (

7) in the case we get:


For illustration purposes, here are the results given by MathematicaTM for the Fourier cosine coefficients of each side of the equation.

In[3]:=FourierCosCoefficient[Cos[x]^2, x, n]
In[4]:=FourierCosCoefficient[(1+Cos[2x])/2, x, n]
Out[4]:=1/2 (DiscreteDelta[-2 + n] + 2 DiscreteDelta[n])

Moreover, its routines do not take advantage of the linearity, as a slight change in the input function prevents them from working:

In[3]:=FourierCosCoefficient[(2+Cos[2 x])/2,x,n]

The strategy followed in our implementation of the package is to internally transform any trigonometric function into its canonical form, so it becomes easy to decide whether or not the input contains an expression whose Fourier coefficients have singular values using the pattern matching capabilities of the Maxima CAS, as described in the following section.

Piecewise-defined functions. Let us point that, in applications, many of the functions we deal with are piecewise defined, and this poses its own challenges. Maxima currently does not have built-in support for piecewise-defined expressions which are of particular importance in engineering and physics applications. We developed another package called piecewise automatically loaded by pdefourier to deal with this kind of functions. In particular, we are able to detect the parity of a piecewise-defined function defined on an interval

by the means of comparison of subintervals allowing us to simplify the computation of Fourier coefficients of odd or even functions. Also, we take advantage of the linearity and interval addition properties of the integral; because of this, our pattern matching rules work even with piecewise-defined functions and equivalent trigonometric expressions as we will see in Section


3 Computation of the Fourier coefficients and series

The aim of this section is to explain the strategies followed in the implementation of our package to tackle the different challenges we have showed in the previous section.

3.1 Piecewise-defined functions

Inherited from Lisp, Maxima’s main data structures are lists, thus the natural way to proceed was to translate a piecewise-defined function to a list and viceversa. Having a way to convert a piecewise expression to a list, it is straightforward to write some other functions to make operations between them or to compute the derivative and the integral of a piecewise expression written as a list. Although the package piecewise is able to work in more general settings, we restrict ourselves here to show only how the package works in a bounded interval. The remaining details are available in the package documentation. For illustration purposes, consider an if-else expression of the form:

if x>=a and x<=a then expr elseif ... elseif x>a and x<=a then expr

We impose the restrictions to avoid potential issues with the detection of valid intervals in all the equivalent ways of writing them, and to do it efficiently. The output returned by pw2list has the form:

The approach of internally working with lists instead of if-else expressions is particularly handy since the Maxima language supports the functional programming paradigm so once we have proper detection of special cases with pattern-matching, we can easily map the routines to each subinterval.

3.2 Computation of Fourier coefficients

As mentioned in the preceding Section, the product of a polynomial by a trigonometric function gives rise to singular values when determining Fourier coefficients, having as particular cases trigonometric functions whose wave number is a multiple of . Because of the linearity of the integral, using trigonometric canonical forms it is enough to detect the following patterns:

To do so, we used the Maxima built-in commands defmatch and matchdeclare. Then, two similar strategies were followed depending on the input expression being piecewise-defined or not. Algorithm 1 describes the steps followed in our implementation to calculate the Fourier coefficients. Of course, when the input function is piecewise-defined the procedure described in the algorithm must be applied to each subinterval of the domain.

Input: expression expr; variable var, semi-length of interval L
Output: A list of Fourier coefficients of the form
l.s.v=list of singular values
Apply simplification functions to convert trigonometric expressions appearing in expr into their canonical form. Expand expr fully and convert the expanded expression into a list . Apply a heuristic routine with pattern matching capabilities to each element of the list to compute and apply an auxiliary function which searches the set of indices having singular values, if any, of and store it in a list . Sum over both to obtain the final answer for and compute separately the Fourier coefficients of the indices appearing in storing the result in the list of singular values. The list of singular values, if any, is written as follows , otherwise, an empty list [ ] is returned.
Algorithm 1 Computation of Fourier coefficients

The same ideas also apply for the case of the complex, sine or cosine coefficients; the only difference is in the way the output is returned to the user (see Table 1).

Type Command Answer format List of s.v format
Trigonometric fouriercoeff(expr, var, L)
Complex cfouriercoeff(expr, var, L)
Cosine fouriercoscoeff(expr, var, L)
Sine fouriersincoeff(expr, var, L)
Table 1: Output formats for the Fourier coefficients (s.v= singular values).

Also notice that, by using expressions instead of declared functions, we are capable to treat arbitrary, non-evaluated functions such as a symbolic . The capabilities of Maxima, however, allows us to evaluate these when required.

3.3 Fourier series

Following our policy of efficiency, Fourier series are obtained using an expansion routine of the Fourier coefficients. We can suspect that Mathematica™  does not perform an expansion of the coefficients to obtain the Fourier series, but computes them one by one inside the FourierSeries function, for two reasons: First, because Mathematica™  only works with truncated series and so a general expression for the coefficients is not needed and second, because this approach avoids the problems associated to simplification of symbolic integrals using assumptions. Nevertheless this strategy is inefficient because it does not take advantage of previously computed Fourier coefficients. As we have seen before, Mathematica™  can not obtain correctly the cosine coefficients of , but if we compute the first five terms of the cosine series of this function we get the right answer:


In our case the upper limit of summation can be a positive integer or infinite. In the first case, a truncated series is returned; in the second, a symbolic series is displayed. In Table 2 we summarize the syntax for the expansion routines along with the different Fourier series commands, see also Table 1 for comparison.

Expansion routine Series
fouriercoeff_expand(list of coeff,var,L,N) fourier_series(expr,var,L,N)
cfouriercoeff_expand(list of coeff,var,L,N) cfourier_series(expr,var,L,N)
fouriersincoeff_expand(list of coeff,var,L,N) fouriersin_series(expr,var,L,N)
fouriercoscoeff_expand(list of coeff,var,L,N) fouriercos_series(expr,var,L,N)
Table 2: Expansion routines and Fourier series syntax in pdefourier.

3.4 Numerical aspects

The determination of Fourier coefficients in the case of the wave equation requires computing the zeros of Bessel’s functions. Unfortunately, Maxima ’s numerical capabilities are rather limited, in particular, it does not have any built-in function for this task (although it has implemented the Bessel functions of first and second kind). Thus, it has been necessary to write a dedicated function for computing these zeros, and in the process of doing so several surprises regarding the implementation of Bessel functions in commercial CASs arose.

Our solution for determining the zeros of and is based on the papers [2] and [4], which allows for a simple and efficient implementation of their algorithms, the main step being the computation of eigenvalues of certain matrices (we do this by importing the function dgeev from Netlib’s lapack package). As we are interested exclusively in real zeros, we restrict the order in to (Lommel’s theorem), but otherwise we impose no restrictions on . Thus, for example, the third zero of is returned by the function BesselJzeros (mimicking the names in Maple™  and Mathematica™ ):

(%i1) load(pdefourier)$

and the third zero of is


While Maple™  gives the result (even with better precision) effortlessly, to our surprise Mathematica™  could not calculate it, neither did the web-based engine Wolfram Alpha. Probably, this is due to the fact that they are using some variant of MacMahon formula, which fails for higher orders (see the comments in

[10]). Other alternatives, such as Halley’s method, also put a bound on the order, as it occurs in Matlab [8].

As mentioned, the algorithms in [2, 4] can be used to generate the zeros of . This is done by the function BesselJdiffZeros. The following table gives the first 5 zeros of the derivatives of for , and is to be compared with the one in [11]:


Up to our knowledge, Maple™  does not have a similar function for computing the zeros of derivatives of Bessel functions.

3.5 Solution of PDEs

In developing our library, we focused on solving the three main second-order linear partial differential equations with constant coefficients, namely, the heat, wave and Laplace equations. The mathematical details can be found in any standard textbook on PDEs or Fourier Analysis (for instance, see [6, 9]). Here, we just briefly discuss some technical aspects of the implementation, a list of examples will be given in the next Section.

It has been mentioned that Fourier series are obtained by performing an expansion of the Fourier coefficients. This approach is multipurpose. It is also useful to code the solution of the PDEs in an easier way. Although the package can solve the heat equation with a heat source , we want to illustrate how the expansion sub-routines facilitate the process of finding solutions with a simpler instance of the general equation. Consider the following IBVP:

the solution is then given by

Notice that for our purposes, it will be sufficient if we create a list of the form where and then use the expansion routine corresponding to a Fourier sine series on the interval and the space variable . For instance, if we consider , a solution to the IBVP can be obtained in a few lines of code:

(%i2)Ψg(x):=if (0<=x and x<=1) then x^2*(1-x)$

Here the list of singular values was empty (obviously), and no further work was required. The library command mixed_heat automates all the process in the most general setting.

Similar methods to obtain solutions for the three equations with different types of boundary conditions and domains have been implemented, dealing with general expressions (piecewise-defined or not) and taking care of the possible singular values in the coefficients.

4 Some examples

4.1 Fourier coefficients and series

Let be as in (7). Then we have the following:


Indeed, notice that .

The case where Mathematica™  fails to compute correctly the cosine coefficients of equivalent trigonometric expressions (8), is readily solved by our package:


Next, we offer an example of a piecewise-defined function having singular values. We will use the command cfouriercoeff to get the complex Fourier coefficients, so that we can compare our answer to the one returned by Mathematica™.

(%i5)Ψf(x):=if x>= -%pi and x<0 then 0 elseif x>=0 and x<%pi then sin(3*x)$

However, Mathematica TM is not able to detect the singular value of the coefficient when :

In[8]:=f[x]:= Piecewise[{{0,-Pi<=x<0}, {Sin[3x],0<=x<Pi}}]
In[9]:=FourierCoefficient[f[x], x, n]

Now, we consider an example of how Fourier series are displayed symbolically. If on we get the answer in the form used in textbooks:


Lastly, we show an example of how to use the expansion routines to obtain the Fourier series, truncated or not, and how they handle singular values in the Fourier coefficients when displaying an infinite series.

(%o10)Ψ"The sum is over |N-{3}"

It is important to notice that displaying infinite series correctly has been a source of troubles in different CAS. This is mainly due to the fact that they are not evaluated, only displayed symbolically, so a simplification evident for a human might not be performed by a CAS.

For illustration purposes, suppose that we want to compute the Fourier series of on the interval . It is obvious that its Fourier series is exactly equal to , because its Fourier coefficients are and we have:

However, infinite sums of expressions containing Kronecker delta functions are really hard to simplify to a single term, because it must be verified that the only index that does not vanish is indeed contained in the set of indices over which we are considering the sum. Integration routines in different CAS sometimes return the result for in terms of some sort of Kronecker delta (Mathematica TM returns it in terms of DiscreteDelta), which complicates the task of obtaining the Fourier series by an expansion of the coefficients. Following our approach, we avoid the issue of the evaluation of Kronecker delta functions inside an infinite sum:

(%o11) sin(15*x)

Notice that, as we just saw in the examples above, when the list of singular values is not empty, the library prints a message warning about the set of indices excluded in the infinite sum.

4.2 Partial differential equations

We will show some examples involving each type of equation (parabolic, hyperbolic and elliptic). Although pdefourier is able to solve these equations with many more boundary conditions, the ones given here are meant to give a panorama about how to use the package to solve PDEs.

4.2.1 The linear second-order parabolic equation

Consider the general problem:

The syntax for solving it is

where ord can be inf, for the complete series solution, or a natural number, for the truncated series.

As an example, we will solve a problem with homogeneous boundary conditions


In comparing results with Maple™  2020 and Mathematica™  we imposed a time limit of seconds (real time) to get an answer from both. If no answer was returned in such time, we considered it a failure. Maple™  reached that time limit without giving an answer; however, Mathematica™  solved it giving the same answer as below.

With pdefourier we would do


4.2.2 Laplace equation

Here we consider the Laplace equation on a wedge with Neumann conditions:

Mathematica™  is unable of solving it (tested on versions 12.0 and 12.1). The syntax for Neumann problems in pdefourier is

In this example we will consider333The use of the percent sign is useful in some Maxima’s graphical interfaces to get the symbol displayed as a Greek character. , and an arbitrary function :


4.2.3 Wave equation

For this example, we will consider the wave equation with homogeneous boundary conditions and a driving term, a problem whose general expression is:

The syntax for solving it using pdefourier is

If and , then:


Since no assumptions were made about , the solution corresponds to the case without resonance (). This is Exercise 8.5.2(b) in [6]. For a different example, consider a vibrating clamped circular membrane, which is modeled in polar coordinates as

The general case is handled by the function


where is the maximum order of the Bessel functions desired in the solution and is the number of positive zeros taken into account. In particular, for a circular membrane of unit radius, , initial shape described by , and no initial velocities, we would issue the following commands:


There is a difference between this function and the remaining ones in the library. While we used expressions as arguments before (as sin(5*%pi*x)), here we need to enter declared functions (like f,g). Ultimately, this has to do with the fact that we are using numerical computations, but this is not an essential aspect, and as such will surely be modified in future versions to allow for expressions in the arguments.

We have suppressed the output of the last command, since it contains a lot of terms with the aspect

Obviously, these are numerical artifacts coming from rounding errors (as revealed by exponents like ). In these cases, we can employ the function chop (provided in the package) analogous to the Chop[] encountered in Mathematica™ , to avoid rounding errors in the plot function. By default, chop discards terms of absolute value less than . An additional argument can be supplied, as in chop(expr,exponent), to chop those terms whose absolute value is less that :


The output has a format suitable for graphical representation, or even to create an animation of the motion444Here we use the command wxanimate_draw, provided by the graphical frontend wxMaxima.:

Ψ   surface_hide=true,zrange=[-1.5,1.5],ytics=0.4,xtics=0.4,color=orange,
Ψ   parametric_surface(r*cos(theta),r*sin(theta),subst(t=s,expr),r,0,1,theta,0,2*%pi)

Appendix A Installation of pdefourier

The package pdefourier is available at Once the repository has been cloned or downloaded, the package can be installed by putting a copy of the files in side 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)$ once inside Maxima.


  • [1] M. A. Al-Gwaiz: Sturm-Liouville Theory and its Applications. Springer Verlag, London, 2008.
  • [2] J. Grad and E. Zakrajsek: Method for evaluation of zeros of Bessel functions. J. Inst. Math. Appl. 11 (1973) 57-72.
  • [3] G. Hellwig: Differential Operators of Mathematical Physics. Addison-Wesley, Reading (Ma.), 1967.
  • [4] Y. Ikebe,Y. Kikuchi, I. Fujishiro: Computing zeros and orders of Bessel functions J. of Comp. and Appl. Math. 38 Issues 1-3 (1991) 169-184.
  • [5] F. John: Partial Differential Equations. Springer Verlag, New York, 1982.
  • [6] R. Haberman: Applied Partial Differential Equations. Pearson Education, Upper Saddle River (NJ), 2013.
  • [7], Maxima, a Computer Algebra System. Version 5.42.2 (2019).
  • [8] J. Nicholson (2020). Bessel Zero Solver (, MATLAB Central File Exchange. Retrieved March 31, 2020.
  • [9] J. S. Walker: Fourier Analysis. Oxford University Press, New York, 1998.
  • [10] G. N. Watson: The zeros of Bessel functions. Proc. of the Royal Soc. of London. Series A94 (1918) 190-207.
  • [11] E. W. Weisstein: Bessel Function Zeros. From MathWorld–A Wolfram Web Resource. Retrieved March 31, 2020.