## Abstract

A new method, combining complex analysis with numerics, is introduced for solving a large class of linear partial differential equations (PDEs). This includes any linear constant coefficient PDE, as well as a *limited* class of PDEs with variable coefficients (such as the Laplace and the Helmholtz equations in cylindrical coordinates). The method yields novel integral representations, even for the solution of classical problems that would normally be solved via the Fourier or Laplace transforms. Examples include the heat equation and the first and second versions of the Stokes equation for arbitrary initial and boundary data on the half-line. The new method has advantages in comparison with classical methods, such as avoiding the solution of ordinary differential equations that result from the classical transforms, as well as constructing integral solutions in the complex plane which converge exponentially fast and which are uniformly convergent at the boundaries. As a result, these solutions are well suited for numerics, allowing the solution to be computed at any point in space and time without the need to time step. Simple deformation of the contours of integration followed by mapping the contours from the complex plane to the real line allow for fast and efficient numerical evaluation of the integrals.

## 1. Introduction

A new method for analysing boundary-value problems for linear and integrable nonlinear partial differential equations (PDEs) was introduced by one of the authors in Fokas (1997) and implemented for a variety of problems in Fokas & Sung (2005). For linear evolution PDEs with spatial derivatives of arbitrary order formulated either on the half-line, 0<*x*<∞, or in the finite interval, 0<*x*<*L*, this method expresses the solution as an integral in the complex *k*-plane (the Fourier plane). This integral involves the *x*-Fourier transform of the initial condition and certain *t*-transforms of the given boundary conditions (Fokas 2002; Fokas & Pelloni 2005; Pelloni 2005). It should be noted that this method constructs novel representations for even simple problems that are traditionally solved in terms of classical transforms. For example, for the heat equation on a finite interval with Dirichlet boundary conditions, the new method yields a novel *integral* representation instead of the classical Fourier sine *series* representation. The new representation, in contrast to the classical one, is uniformly convergent at the boundaries.

It will be shown here that these novel integral representations are suitable for the numerical evaluation of the solution. This is a consequence of the fact that it is possible, using simple contour deformations in the complex *k*-plane, to obtain integrals involving integrands with a strong decay for large *k*.

In order to present this novel approach in its simplest form, we will only consider initial and boundary conditions for which the associated Fourier and *t*-transforms can be computed analytically. In this case, the numerical implementation consists only of computing an integral in the complex *k*-plane involving a decaying integrand for large *k*.

This new technique will be illustrated by integrating numerically certain initial-boundary-value problems on the half-line for the heat equation(1.1)for the first version of the Stokes equation(1.2)and for the second version of the Stokes equation(1.3)Equations (1.2) and (1.3) are the linear limits of the celebrated Korteweg–de Vries equation and equation (1.3) corresponds to the case of dominant surface tension.

The paper is organized as follows: §2 gives a brief overview of the analytical method of Fokas (2002) and derives the integral representations of the solutions for the above PDEs; §3 presents the technique for the numerical integration of the integrals presented in §2; §5 summarizes the methodology and §6 the advantages and limitations of the method. Appendix A gives the Matlab and Mathematica code for implementing the novel numerical method for an example of the heat equation. Appendix B compares the new method with the classical Laplace transform method.

## 2. Novel integral representations on the half-line

Here, we give only a brief overview of the analytical method that is used to derive the novel integral representations of the solution. A detailed discussion is given in Fokas (2002). The one-parameter family of solutions(2.1)where *w*(*k*) is a polynomial of degree *n* and Re *w*(*k*)≥0 for *k* real (the latter restriction ensures that the initial-value problem for the associated PDE is well posed) is admitted by the following linear evolution PDE:(2.2)For example, for the heat equation, *w*(*k*)=*k*^{2} and . We next seek a function *X*(*x*, *t*, *k*) such that we can rewrite the PDE as the following family of divergence forms:(2.3)Simplifying this equation and replacing *q*_{t} by , we find(2.4)The r.h.s. of this equation involves derivatives of order up to *n*−1, *n* being the order of the polynomial *w*(*k*). Hence, we can write *X*(*x*, *t*, *k*) as(2.5)Equating (2.4) and (2.5) for *X*(*x*, *t*, *k*) and letting −i*∂*_{x}=*l* (an arbitrary complex parameter), we can compute the coefficients *c*_{j}(*k*) in terms of *w*(*k*)(2.6)where *k* and *l* are arbitrary complex parameters. By employing Gauss's theorem, equation (2.3) yields in Fourier space the time evolution of the solution in terms of its behaviour on the time–space boundary, as given by the Fourier transform of the initial condition and certain *t*-transforms of the boundary values (see (2.9)). Then, taking the inverse Fourier transform, we find(2.7)where and are defined as follows (Fokas 2002): is the Fourier transform of the initial condition, i.e.(2.8)and represents certain *t*-transforms of *X*(0, *t*, *k*), namely the transforms of the boundary values at *x*=0,(2.9)(2.10)with *c*_{j}(*k*) defined in (2.6).

The representation (2.7) *cannot* be used directly for the solution of a given initial-boundary-value problem because some of the boundary values are unknown. Indeed, the number of boundary conditions *N* required for a well-posed problem is given in Fokas (2002) by(2.11)where *c* is a positive constant and *α*_{n} denotes the coefficient of *k*^{n} in *w*(*k*).

Suppose that *q*(0, *t*) and the first *N*−1 derivatives of *q*(0, *t*) are given as boundary conditions. Then, the *t*-transforms of the remaining unknown derivatives can be obtained by solving a system of *n*−*N* algebraic equations. These equations are obtained by replacing *k* with *ν*(*k*) in the equation(2.12)where *ν*(*k*) denotes any of the roots of the equation *w*(*ν*)=*w*(*k*), which map *D*^{−} into *D*^{+}, with *D*^{−} denoting the part of the domainin the lower half complex *k*-plane and *D*^{+} the part of *D* in the upper half complex *k*-plane. We *emphasize* that the boundary values are never needed, only their *t*-transforms are required, which can be obtained from the transforms of the given boundary conditions and from (2.12).

The solution can then be represented in the form(2.13)

In what follows, we will apply the above construction to equations (1.1)–(1.3).

### (a) The heat equation

Substituting the expression (2.1) into equation (1.1), we find *w*=*k*^{2}. Hence equation (2.6) with *n*=2 implies(2.14)giving for equation (2.10)(2.15)Using , it follows that(2.16)

For mapping *D*^{−} into *D*^{+}, the equation *ν*^{2}=*k*^{2} implies *ν*=−*k*. Replacing *k* by −*k* in equation (2.12), we find(2.17)

We now consider the Dirichlet problem(2.18)where the functions *q*_{0}(*x*) and *g*_{0}(*t*) have appropriate smoothness. The first of these functions has appropriate decay, and the functions *q*_{0} and *g*_{0} are compatible at *x*=*t*=0, i.e. *q*_{0}(0)=*g*_{0}(0).

In this case we solve equation (2.17) for (the *t*-transform of the unknown boundary derivative) and substitute the result in the expression for given in (2.15). This yields(2.19)

Hence, if *q*(*x*, *t*) satisfies the heat equation (1.1) and the initial and boundary conditions given by equations (2.18), then the solution for {0<*x*<∞, *t*>0} is given by(2.20)where *∂D*^{+} is the boundary of the domain *D*^{+} defined in (2.16) and depicted in figure 1, and and are given by(2.21)

(2.22)

### (b) The first version of the Stokes equation

For equation (1.2), *w*(*k*)=−i*k*^{3}. Hence by (2.6),(2.23)In this case(2.24)Now, the equation *ν*^{3}=*k*^{3} implies *ν*=*αk* and *ν*=*α*^{2}*k*, . Furthermore, if *k*∈*D*^{+} then and (figure 2). Hence, replacing *k* by *αk* and *α*^{2}*k* in equation (2.12), we find the following two equations coupling the three boundary values:(2.25)(2.26)Thus, for a well-posed problem, we need one boundary condition, i.e. *N*=1.

We now consider the problem defined by equations (2.18). In this case we solve equations (2.25) and (2.26) for and and substitute the results in the expression for defined by the second of equations (2.23).

Hence, if *q*(*x*, *t*) satisfies equation (1.2) and the initial and boundary conditions given by equations (2.18), then the solution for {0<*x*<∞, *t*>0} is given by(2.27)where ; *∂D*^{+} is the boundary of the domain *D*^{+} defined by the first of equations (2.24) and depicted in figure 2; and and are defined by equations (2.21) and (2.22), respectively.

### (c) The second version of the Stokes equation

For equation (1.3), *w*(*k*)=i*k*^{3}, hence(2.28)In this case(2.29)

If then *αk*∈*D*^{−}, and if then *α*^{2}*k*∈*D*^{−} (figure 3). Thus, for the determination of for , we replace *k* by *αk* in equation (2.12),(2.30)similarly, for the determination of for , we replace *k* by in equation (2.12),(2.31)Each of equations (2.30) and (2.31) is one equation coupling the three boundary values. Thus for a well-posed problem, we need two boundary conditions, i.e. *N*=2.

We now consider the initial-boundary-value problem defined by equation (2.18), supplemented with the equation(2.32)where the function *g*_{1}(*t*) has sufficient smoothness and .

In this case we solve equations (2.30) and (2.31) for and substitute the resulting expressions in equation (2.28).

Hence, if *q*(*x*, *t*) satisfies equation (1.3) and the initial and boundary conditions are given by equations (2.18) and (2.32), then the solution for {0<*x*<∞, *t*>0} is given by(2.33)where the contours are the boundaries of the domains and defined by the first expression in (2.29) and depicted in figure 3; and are defined by equations (2.21) and (2.22), respectively; and is defined by(2.34)

In equation (2.33), we have simplified the expressions for by employing the identity 1+*α*+*α*^{2}=0.

## 3. Numerical evaluations on the half-line

In this section, we present a novel technique for numerically evaluating the integral representations of the solutions to the initial-boundary-value problems discussed in §2. The main ideas common to all cases are as follows.

To perform simple contour deformations in the complex

*k*-plane such that the deformed integration paths are in regions where the integrands decay exponentially for large*k*. This yields rapid convergence of the numerical scheme.For algorithmic convenience and simplicity, to make a change of variables, which maps the contours from the complex plane to the real line.

Below, we will discuss in detail how we choose the integration path in each case.

### (a) Solutions of the heat equation

Recall the solution to the heat equation on the half-line for {0<*x*<∞, *t*>0},(3.1)

Let us first consider the integral whose contour runs along . Using (2.22), the first term of the integral can be written as(3.2)Since exp[i*kx*] is bounded and analytic for Im *k*>0 and with *t*≥*s* is bounded and analytic for Re (*k*^{2})>0, it follows that, for this term, the contour can be deformed to any contour *L* in the unshaded domain of the upper half complex *k*-plane of figure 1. The next term involves the factors e^{ikx} and which are bounded and analytic for Im *k*>0, and the factor which is bounded and analytic for Re *k*^{2}>0. Hence, the contour *∂D*^{+} for this term can also be deformed to *L*.

Splitting the term defined by (see (2.8)) into two terms and then substituting these terms into the integral along the real axis in (3.1) results in the following integrand:(3.3)

In the first term, *x*−*ξ*≥0, thus for this term the contour along the real axis can be deformed to *L*. For the second term in (3.3), *x*−*ξ*≤0, thus for this term the contour along the real axis can be deformed to a contour *L*^{−} in the unshaded domain of the lower half complex *k*-plane.

The term is in general analytic only for Im *k*<0; however, depending on the properties of *q*_{0}(*x*), it is sometimes possible to extend the domain of analyticity to the upper half of the complex *k*-plane. For example, this is possible if *q*_{0}(*x*) is such that *q*_{0}(*x*)e^{αx}, *α*>0, is square integrable on [0,∞). If *q*_{0}(*x*) belongs to this restricted class, then the contour along the real axis can also be deformed to the contour *L*.

(3.4)where *a* and *b* are real numbers. Then,(3.5)Hence, equation (2.20) yields(3.6)where the contour *L* that is actually used in the computation is depicted in figure 4*b*.

Technically, any valley-to-valley path in the upper half plane of figure 4*a*, where and , will give exactly the same result. However, in order to have rapid convergence for large *k*, it is natural to choose a path which for large *k* aligns with the rays of steepest descent, where the integrand decays most rapidly. In the current problem, these lie at . Hyperbolas are simple curves that possess the property of having asymptotic directions and are therefore a natural choice for an integration path.

There are many ways to parametrize hyperbolas. One way that is particularly well suited for our work in the complex plane is to use the analytic function(3.7)which maps the points *θ* on the real line to hyperbolas *k*(*θ*) in the complex plane. The factor i multiplying the sin in (3.7) rotates the contour *π*/2 degrees in the complex plane, while *γ* (set to 1 in this case) is a scaling factor and *α*=*π*/8 so that the asymptotes of the hyperbola *k*(*θ*) lie along the path of steepest descent. An equispaced grid in the *θ*-direction corresponds to points along the hyperbola that are concentrated near the origin and become exponentially sparse for large *k* (figure 4*b*). This type of mapping is particularly well suited for the numerical integration of functions that feature their greatest variations near the origin.

Using the mapping (3.7), the integral (3.6) becomes(3.8)Now that the integral has been mapped to a straight line (the real line in this case), the simple trapezoidal rule can be applied to (3.8), giving exponential accuracy, since the integrand is analytic and decays rapidly on an unbounded domain. This type of numerical quadrature has been used to evaluate a variety of contour integrals that appear in applied mathematics, from the inversion of Laplace integrals to fast methods for computing special functions and matrix exponentials (Sidje 1998; Gil *et al*. 2003; Lopez-Fernandez & Palencia 2004; Weideman 2006; Weideman & Trefethen 2007).

Before numerically computing (3.8), we first discuss a truncation estimate for the required integration interval, which is an estimate for *N* in truncating to , where *f*(*θ*) is the integrand in (3.8). It is easy to show that the imaginary part of the integrand is antisymmetric and thus will not contribute to the integral. If we expand *f*(*θ*) into real and imaginary parts, simplify and then take the leading order term, we find that the leading behaviour of |*f*(*θ*)| is given by *g*(*θ*), where(3.9)Plotting |*f*(*θ*)| against |*g*(*θ*)| in figure 5, we see that this estimate is very accurate. In addition, |*g*(*θ*)| has decayed halfway from its maximum to zero near the ends of the interval when , i.e.(3.10)For *x*=10^{−13}, this yields . For *x*=10^{−43}, the transition occurs around *θ*=±100, giving an integration interval of *θ*∈[−101,101]. Thus, for any *x*>0, we can easily and accurately compute the integral since the integration interval increases at a slow logarithmic rate. However, for *x*=0 the integration interval is infinite and any truncation will result in the integral not converging exactly to the boundary condition.

The above discussion shows that the integration interval depends on *x* (the dependence on *t* is innocuous as noted in (3.9)). Furthermore, figure 6, where the real part of the integrand for decreasing *x* is plotted, indicates that the contributions to the integral come from two ‘humps’ of constant height near the origin which are not a function of *x* and two ‘humps’ of constant height which move out from the origin as *x*→0.

The fact that the real part of the integrand contains a lot of ‘fine structure’ near the outer part of its support renders the use of the trapezoidal rule cumbersome: as *x*→0 the number of evaluation nodes increases dramatically since they are equispaced. Hence, an adaptive quadrature method is more suitable. We next discuss such an approach.

#### (i) An alternative approach

Equation (3.8) defines an ordinary integral with an exponentially decaying integrand as *θ*→±∞. Any language with a built-in numerical integrator, such as Mathematica, Maple or Matlab, can provide a much simpler approach to evaluating the integral. This is particularly true for the integrand of (3.8) which, as discussed earlier, exhibits rapid variations for small *x*. Since these numerical integrators use adaptive quadrature schemes that ‘sense’ these variations, they are most suitable for the integrals we are dealing with. In the table below, the trapezoid rule, programmed in Matlab, is compared with the Mathematica NIntegrate command, displaying the time in seconds it took to run the code for evaluating the integral at (*x*, *t*)=(0.5, 0.5), (*x*, *t*)=(1×10^{−3}, 1×10^{−3}) to six digits of accuracy and to produce figure 7.

The evaluation of a single point takes the same amount of time in either program. To produce figure 7 takes about 14 s longer in Mathematica than Matlab. However, in Mathematica one does not have to worry about the number of evaluation nodes needed; furthermore, it is simpler than the Matlab code (see appendix A). Thus, we will use Mathematica in the examples analysed in this paper (table 1).

### (b) Solutions to q_{t}+q_{xxx}=0 on the half-line

Recall the general solution to (1.3) on the half-line given by the representation (2.27),(3.11)The term involving can be treated as the corresponding term of the heat equation in the previous subsection. However, since the real axis is *not* surrounded by the domain satisfying , it is *not* possible, by splitting , to deform the real axis to a contour in the lower half of the complex *k*-plane. On the other hand, if *q*_{0}(*x*) belongs to a restricted class, then the real axis can be deformed to the same contour that *∂D*^{+} will be deformed to in the upper half plane; similar considerations apply to the terms involving and .

Let *q*_{0}(*x*) and *g*_{0}(*t*) be defined by equation (3.4). Then equation (2.27) becomes(3.12)

Regarding the determination of the integration path *L*, this case is almost identical to that of the heat equation except that the rays defining the boundary of the domain *D*^{+} are at arg *k* equals *π*/3 and 2*π*/3. As a result, we will use the same change of variables defined by equation (3.7), except that *α* will now be *π*/6. The integration path and the modulus of the integrand are depicted in figure 8.

After using the substitution , the built-in numerical integrator in Mathematica is used to evaluate the r.h.s. of (3.12). The solution is plotted in figure 9.

### (c) Solutions to *q*_{t}−*q*_{xxx}=0 on the half-line

The rays and 2*π*/3 in the representation of the solution (2.33) can be deformed to the unshaded domain of the upper half complex *k*-plane of figure 3, but if exp[i*kx*] and exp[−i*k*^{3}*t*] are treated *separately*, the real axis *cannot* be deformed to the unshaded domain of the lower half of the complex *k*-plane, since exp[i*kx*] is unbounded for Im *k*<0. However, this problem can be bypassed if the term exp[i*kx−*i*k*^{3}*t*] is treated as a *single* term. In a follow-up paper for the finite interval case, we will also address the issue of inhomogeneous boundary conditions for this PDE on the half-line as certain novel treatments are needed in order to numerically handle the deformation of the contours.

Let(3.13)where *a* is a positive real number. Then,(3.14)Hence,(3.15)The regions of the *k*-complex plane where determines the three sectors (unshaded areas in figure 3) where and hence each of the integrands in (3.15) is bounded. The contours, and , and the real axis can therefore be deformed to run from the valley of one sector into the valley of another where the integrands are exponentially small for large *k*. The only precaution that must be taken is that the valleys widen or narrow as a function of *x* and *t*. To illustrate this point, let us plot in figure 10 the absolute value of for two different values of *x* and *t*, (*x*=15, *t*=15) and (*x*=10, *t*=0.1), keeping the scale the same on both plots for comparison.

In order to overcome this difficulty, we require that the integration path *k*=*f*(*θ*) runs straight down the centre of the valleys and also passes through the saddlepoints of the surface defined by the integrand. The saddlepoints for all the integrands occur when , i.e. . We then define a shift to our contour by *k*=*f*(*θ*)+*α* and ask for what value of *θ* and *α* will *k* be equal to .

#### (i) Deformation of and

The shift to the contour in this case is trivial since we only have to pass through one saddlepoint that occurs when *θ*=0 in each case. As a result, for (3.16)and for (3.17)

The factor of rotates the contour *π*/6 degrees to orient it along the correct sectors (for , the addition of the minus signs to this factor simply reflects it across the imaginary axis). The contour for the deformation is plotted in figure 11*a* and the contour for the deformation is plotted in figure 11*b*. The poles do not occur in the sectors where the contours run.

#### (ii) Deformation of the real axis

This case is more complicated since the contour is deformed into the lower half plane and must pass through both saddlepoints . Our original deformed contour is . Splitting this expression into real and imaginary parts, adding the shift *α* and setting the expressions equal to , we have two equations in two unknowns, *α* and *θ*. Solving the system, we arrive at the contour(3.18)

In this case there is an additional difficulty due to the pole at *k*=i*a*^{2}: if we fix *x* and let *t* become small, the contour gets shifted upward, interacting with the pole and adding an unwanted contribution. To illustrate this difficulty, we plot in figure 12*a* and figure 12*b* the absolute value of , as well as the contour for *x*=1, *t*=1 and then for *x*=1, and *t*=0.01.

This difficulty can be overcome by subtracting out the pole to obtain a singularity-free integrand. This is done by subtracting from the function, , a function that decays in the lower half plane (such as e^{−ik}) but has the exact same pole character. The function is derived by enforcing that the coefficients (*a*_{−1}, *a*_{−2}, *a*_{−3}) for the three components of the pole in the Laurent expansion of the integrand (i.e. ) vanish when it is subtracted. Such a function is given by(3.19)

Figure 12*c* and figure 12*d* show the surfaces with the pole removed and the corresponding integration path. The integrand to be evaluated is now , where the mapping (3.18) has been used.

Now that all three contours have been mapped to the real axis, the integrands decay exponentially for large *k* and do not involve any poles; we can numerically evaluate the integrals using Mathematica. The solution is given in figure 13.

#### (iii) A note on time–space corner singularities

One of the striking features of the solution is the wave pattern that emanates from the corner of the domain (*x*=0, *t*=0). This phenomenon results from the fact that unless the initial and boundary data are compatible for all orders (i.e. unless they satisfy an infinite set of compatibility conditions), the solution will feature irregularities that emanate from the corner and behave according to the nature of the given PDE (i.e. diffuse for the heat equation; propagate for the current problem; Boyd & Flyer 1999). In our case, the initial condition (IC) and boundary conditions (BCs) do match at the corner, satisfying the lowest order compatibility condition. However, there exists an incompatibility in the next order, *q*(BC)_{t}≠*q*(IC)_{3x}, since *q*(BC)_{t}=0 and *q*(IC)_{3x}=−27/2. All higher order conditions are derived by simply taking derivatives with respect to the time of the PDE and substituting in the IC and BCs. An extensive study on the nature of such singularities and their numerical implications for dissipative, dispersive and convective PDEs is given in Flyer & Swarztrauber (2002) and Flyer & Fornberg (2003*a*,*b*).

## 4. An example of a finite interval case

We will present details of the discussion for the finite interval case in a follow-up paper. Hence, we discuss only briefly a simple example.

For equation (1.1) the numerical evaluation of *q*(*x*, *t*) involves steps similar to those used in §3*a*. For example, let the initial and boundary conditions be(4.1)(4.2)(4.3)Then, for the boundary condition at *x*=0,(4.4)and for that at *x*=1(4.5)The solution on 0<*x*<1 and *t*>0 is then given by(4.6)

The details of the derivation will be presented in a follow-up paper that will focus on the finite interval case as well as on inhomogeneous boundary conditions for the second version of Stokes equation. As with the heat equation on the half-line, the mapping is used to deform the *∂D*^{+} contour to the real line and for the *∂D*^{−} contour. For the domain *x*∈[0,1], *t*∈[0,2*π*], only the short integration interval *θ*=±10 was needed before the integrand became negligible. The solution is given in figure 14.

## 5. Summary of methodology

The steps for the analytical derivation of the solution can be summarized as follows.

Compute the Fourier transform of the initial condition for 0≤

*x*≤∞.Determine the number of necessary boundary conditions , where

*n*is the degree of the polynomial*w*(*k*) appearing in (2.1).Compute the

*t*-transforms of the given boundary data.Find the roots of

*w*(*ν*)=*w*(*k*) which map the part of the domain in the lower half plane to the upper half complex*k*-plane. For example, for the heat equation the transformation*ν*(*k*)=−*k*maps the domain*D*^{−}to the domain*D*^{+}as depicted in figure 1. Similarly, for the first version of the Stokes equation, the transformations and map the domains and depicted in figure 2 to the domain*D*^{+}.Substitute these roots in the equation . This will give a system of

*n*−*N*linear algebraic equations for the*t*-transforms of the unknown boundary values, where*c*_{j}(*k*) is defined in (2.6). We emphasize that we do*not*need the unknown boundary values , only their*t*-transforms are required, which can be expressed in terms of the transformed boundary and initial data.The general solution is then given by

(5.1)

Once we obtain the solution, the methodology for numerically integrating the contour integrals for the heat and first version of the Stokes equation is as follows.

Deform the integration paths, i.e. the contours

*∂D*^{+}and the real axis to a hyperbola in the upper half*k*-complex plane by the mapping , choosing*α*and*γ*so that the path aligns with the rays of steepest descent. Real values of*θ*are then chosen as the numerical evaluation nodes.Insert the mapping

*k*(*θ*) into the analytical solution and then evaluate numerically the resulting expression using the simple trapezoidal rule, which will give spectral accuracy for integrands that are analytic and decay exponentially fast on an unbounded domain. Alternatively, numerical integrators built into such languages as Mathematica, Maple or Matlab can be used.

For the second version of the Stokes equation, the same steps as above apply except the following.

The mapping

*k*(*θ*) will need to be shifted by a constant*α*(*x*,*t*) that is determined by finding the saddlepoints of the relevant integrand.For the deformation of the integral path along the real axis, any poles in the upper half

*k*-complex plane will need to be subtracted out as seen in figure 3.This is accomplished by subtracting from the integrand a function such that the resulting expression is free from singularities.

## 6. Conclusion

We have introduced a novel hybrid analytical–numerical method for solving certain linear PDEs. The starting point of this method is the derivation of a novel analytical representation for the solution (Fokas 2002). The advantages as well as the limitations of the method are discussed below.

For evolution PDEs, the main advantage of the new method is that it can be used to compute solutions at arbitrary points in the (

*x*,*t*) plane. Neither time stepping nor spatial differencing is required. In this respect, the new method has some similarities to the Laplace transform technique; however (i) for an evolution PDE with an*n*th order spatial derivative, the Laplace transform involves exp[−*st*+*λ*(*s*)*x*], where*λ*(*s*) solves an*nth order algebraic equation*, whereas our formulation involves , where*w*(*k*) is an*explicit*polynomial of degree*n*. For example for the equation*λ*(*s*) solves the cubic equation*λ*^{3}+*λ*−*s*=0, whereas . (ii) The explicit and analytic dependence on*k*allows us to deform contours which in turn yields exponentially decaying integrals and thus efficient numerical computations. (iii) The Laplace transform requires*t*going to ∞, which is*not*natural for evolution equations (one usually attempts to overcome this difficulty by appealing to causality arguments). Even for problems where*λ*(*s*) can be found explicitly, such as the case of the heat equation, the new method has advantages over the Laplace transform method as shown in appendix B.The novel representation, in contrast to the classical Fourier integral (half-plane) or the Fourier series (finite interval), has the advantage that it converges uniformly at the boundaries.

The semi-analytical nature of the new method has a pedagogical advantage: the usual numerical techniques for spatial discretizations Fornberg (1996), Trefethen (2000) and Boyd (2001) are constructed independent of the analytical treatment of the given PDE. This often raises questions about the reason for teaching students analytical techniques. This should be contrasted with the new method, where the numerical integration is the last step of an approach that is based on the analytical treatment of the given PDE.

Although this method is obviously very efficient, the applications presented in the paper have the following limitations:

involve constant coefficient PDEs,

involve PDEs with derivatives of first order in time, and

require initial and boundary conditions of sufficient simplicity, so that their Fourier transforms can be computed explicitly.

However, among these limitations, it appears that only (i) is a significant limitation. Indeed, the analytical method introduced in Fokas (2002) can be applied to any PDE that admits a Lax pair formulation. This includes linear PDEs with constant coefficients, such as the wave, the Laplace and the Helmholtz equations, as well as a *limited* class of PDEs with variable coefficients (such as the Laplace and the Helmholtz equations in cylindrical coordinates). For these equations it is possible to obtain an integral representation in the complex *k*-plane, which involves the Fourier transforms of the given boundary conditions and of the *unknown* boundary values. For simple problems, the Fourier transforms of the unknown boundary values *can* be eliminated and then the numerical technique introduced in this paper can be immediately applied. For example this approach is implemented in Kalimeris (in preparation) and Spence (in preparation) for the particular problems of the Laplace and modified Helmholtz equations formulated in the interior of either an orthogonal isosceles triangle or an equilateral triangle. For more complicated problems, such as the Laplace and Helmholtz equations in the interior of an arbitrary convex polygon, one must first use the novel numerical technique of Sifilakis *et al*. (in press) to determine the unknown boundary values and then one can apply the numerical technique introduced here (see Smitheman *et al*. in preparation).

Regarding the limitation (iii), there exist approximations currently under investigation which yield an *explicit* integrand for the integral formulated in the complex *k*-plane; furthermore, this integrand has similar analytic properties to those discussed in this paper. There exist powerful methods for the numerical evaluation of Fourier-type integrals when they cannot be computed analytically (see for example Iserles 2004).

We emphasize that the results presented in this paper are of exploratory nature and one can expect further improvements to be made, particularly if these techniques are pursued by other researchers.

## Acknowledgments

The work was supported by NSF grant no. ATM-0620100.

## Footnotes

- Received January 28, 2008.
- Accepted February 28, 2008.

- © 2008 The Royal Society