## Abstract

This article presents an overview of a transform method for solving linear and integrable nonlinear partial differential equations. This new transform method, proposed by Fokas, yields a generalization and unification of various fundamental mathematical techniques and, in particular, it yields an extension of the Fourier transform method.

## 1. Introduction

The Fourier transform has been, over the last two centuries, one of the most important mathematical tools in every branch of science, and particularly in physics, engineering, computing and image processing.

The one-dimensional Fourier transform formulae are the following. Let *f*(*x*) be a smooth function defined on such that lim_{x→±∞}*f*(*x*)=0 sufficiently fast (I will not give the precise requirement on differentiability and decay). Then the direct and inverse Fourier transforms are defined by(1.1)The function is known as the *Fourier transform* of *f*(*x*).

An important application of the Fourier transform is the construction of the solution of initial value problem for the basic evolutionary partial differential equations (PDEs) of mathematical physics. These equations model processes evolving in time from a given initial state. One fundamental example is the initial value problem for the heat equation, which models how the heat diffuses starting from a certain initial temperature. Mathematically, this is the problem of finding the function *T*(*x*, *t*) such that(1.2)where *x* is a space variable; *t* denotes time; and the solution *T*(*x*, *t*) represents the temperature, at time *t* and position *x*, of an infinite one-dimensional object, which is initially at temperature *T*_{0}(*x*).

The solution of this problem can be found in two steps using the Fourier transform. Since implies , the solution scheme isThe above approach works also when the second derivative *T*_{xx} appearing on the right-hand side of equation (1.2) is replaced by any linear combination of *x*-derivatives of arbitrary order. Indeed, given in general a polynomial such that Re *ω*(*k*)≥0 when *k*∈, the initial value problem(1.3)can be solved by the following algorithmic scheme:(1.4)This method of solution is general, elegant and useful. In particular, we note that the expression for the solution is ‘spectrally decomposed’: the *x* and *t* dependence of the solution is explicit, while the so-called *spectral function* does not depend on *x* or *t*.

The main aim of this paper is to discuss the *Fokas transform* (Fokas 1997, 2001; Fokas & Sung 2005), a generalized Fourier transform that can yield the solution of more complicated problems, retaining the generality and elegance of equation (1.4). Within the general framework we present, the Fourier transform can be generalized to solve certain nonlinear PDEs, and extended to express *effectively* the solution of boundary value problems for linear evolution PDEs such as equation (1.3), posed on bounded domains that are either fixed or time dependent. In this respect, we note that the application of the Fourier transform to boundary value problems has significant limitations. For the simplest example, consider the heat equation (1.2), posed on the half line *x*>0. To specify the problem, in addition to the initial condition, we must prescribe one boundary condition, e.g. *T*(0, *t*)=*f*(*t*), *t*≥0. Taking the Fourier transform of the equation, we now obtain . In this expression, *T*_{x}(0, *t*) is not known, hence simply taking the inverse Fourier transform of this does not yield an effective solution representation. In this particular case, the problem can be circumvented by considering the *sine transform*, rather than the Fourier transform. By using the sine transform, it is possible to find an expression for *T*(0, *t*) in which the unknown boundary condition *T*_{x}(0, *t*) does not appear. Indeed, using this transform, an *effective* representation for the solution of the problem is given by(1.5)However, this formula *cannot be generalized* to more complicated boundaries or boundary conditions, nor to PDEs involving higher order derivatives. Moreover, even as it stands, it presents significant drawbacks. For example, it is *not uniformly convergent* at *x*=0. This implies that it is not possible to verify if *T*(0, *t*)=*f*(*t*) by simply evaluating the integrand at *x*=0. Furthermore, owing to its slow convergence, this integral formula does not provide an efficient algorithm for the numerical evaluation of the solution. In contrast, the formula for the solution of this boundary value problem obtained by the Fokas transform is(1.6)where (figure 1). Although this expression can be reduced to equation (1.5) using straightforward analyticity considerations, it does have several, significant advantages: it is *spectrally decomposed*; it is *uniformly convergent* as *x*→0; and it can be generalized substantially. Indeed, a similar expression can be obtained for boundary value problems posed for linear evolution PDEs *of any order*, in either fixed or time-dependent domains. Moreover, analyticity considerations allow one to deform the contour *Γ* appearing in equation (1.6) in such a way that the integrand on this new contour possesses *fast decay* (figure 1). This fact leads to fast and accurate numerical approximations.

In §2, we illustrate the general formulation of the Fokas method for linear evolutionary problems using the explicit example of a third-order problem. We also comment on the numerical implementation of the formulae obtained.

In §3, we present a class of generalized Fourier transform pairs. In this presentation, an important role is played by the so-called *global relation*, which is a single equation relating all the boundary values of a given boundary value problem. The systematic investigation of this equation reduces the problem of determining the unknown boundary values (the so-called *Dirichlet to Neumann map*) to the problem of inverting an integral transform. For simple boundary value problems, this inversion can be achieved by using the Fourier transform. For more complicated problems, this can be achieved by considering certain generalized Fourier transform pairs. These novel formulae can be derived by performing the spectral analysis of appropriate eigenvalue ODEs. This new approach for deriving integral transforms was originally presented in Fokas & Gelfand (1994) and further developed in Fokas & Pelloni (2000, 2006). The transforms discussed in this section arise in the context of moving boundary value problems.

In §4, we briefly discuss the nonlinear case. Indeed, the methodology leading to the expression (1.6) can be generalized to solve the boundary value problems for an important class of nonlinear evolution PDEs, known as *integrable*. These PDEs are characterized by a remarkable property: they can be written as the compatibility condition of a pair of linear eigenvalue ODEs, one with respect to *x* and the other with respect to *t*, called the associated *Lax pair*. The spectral analysis of the associated ODE in *x* leads to the solution of the initial value problem for these PDEs by the inverse scattering transform. The nonlinear version of the Fokas method is based on the simultaneous analysis of both ODEs forming the Lax pair, and leads to the solution of boundary value problems. Finally, in §5, we indicate current research directions.

## 2. The Fokas transform: a generalization of the Fourier transform

*The spectral analysis of ODEs and a new kind of separability*. The solution expressed by equation (1.4) is based on the separability of the PDE (1.3). This means that, assuming the solution of equation (1.3) is of the form *q*(*x*, *t*)=*Φ*(*x*)*Ψ*(*t*), one can analyse this PDE by using either a transform with respect to *x* or one with respect to *t*. The Fourier transform in *x* is the appropriate transform to solve the initial value problem for this PDE. Indeed, the resulting time evolution has a very simple form and does not depend on unknown data.

The main idea of the method we present here is that one should consider a different kind of separability. The authors of Fokas & Gelfand (1994) observed that it is always possible to associate a linear evolution PDE with a *pair* of first-order eigenvalue ODEs. For example, equation (1.3) is equivalent to the requirement that, for all *k*∈, the following two ODEs are compatible (i.e. *μ*_{xt}=*μ*_{tx}):(2.1)where the coefficients *c*_{j} are determined explicitly by the formulaFor example, equation (1.2), for which *ω*(*k*)=*k*^{2}, this pair of ODEs is(2.2)For the solution of the initial value problem, it is sufficient to use the spectral analysis of the first of these ODEs. This spectral analysis, as we show below, yields the Fourier transform. This is consistent with the fact that, in this case, the Fourier transform in *x* is the appropriate transform for solving the problem. However, in order to solve boundary value problems effectively, *it is necessary to consider both ODEs in this pair simultaneously*. The spectral analysis of the pair yields a transform, which is custom made for solving the given boundary value problem and yields a representation for its solution in the same abstract form as equation (1.4).

We start by presenting a derivation of the Fourier transform based on the spectral analysis of the first of the ODEs (2.1*a*). This derivation is the starting point of the new approach and it requires only the knowledge of some basic complex analysis. Indeed, in what follows, considerations based on analyticity properties and various consequences of Cauchy's theorem will recur many times. We assume some familiarity with these notions.

Consider the ODE (2.1*a*), where *q*(*x*) is an arbitrary function belonging to the Schwartz space () of smooth, rapidly decaying functions defined on (although considerably less regularity is sufficient for this argument). We pose the question of finding a solution *μ*(*x*, *k*) of this ODE that is bounded for all values of the complex parameter *k*, including infinity.

Two particular solutions of equation (2.1) are given byEach of these functions is analytic and bounded in half of the complex *k* plane. Indeed, writing *k*=*k*_{1}+i*k*_{2}, we find . This exponential is bounded for *k*_{2}(*x*−*y*)≥0, hence for *y*≤*x* when *k*∈^{+}, and for *y*≥*x* when *k*∈^{−}. It follows that *μ*^{±} is analytic and bounded in ^{±}, respectively.

Integration by parts shows that, as *k*→∞, *μ*^{±}=*O*(1/*k*). The jump across the contour common to the two regions where *μ*^{±} are bounded, which in this case is , is(2.3)The reconstruction of the function *μ*(*x*, *k*), which is analytic in the entire complex plane except the real axis, and for which the jump condition (2.3) is satisfied, defines a classical problem called a *Riemann–Hilbert problem*, which can be solved in closed form as(2.4)Finally, computing *q*(*x*) using the right-hand side of equation (2.1), we obtaini.e. the transform pair applied to the arbitrary function *q*(*x*). In summary, *the spectral analysis of the ODE* (2.1) *yields the Fourier transform pair applied to the function q*(*x*). This idea can be generalized to other eigenvalue ODEs. The spectral analysis of various cases yields both classical and apparently new transform pairs.

For the purpose of analysing boundary value problems for PDEs, we must perform a *similar* analysis of the *pair* of ODEs defining the Lax pair.

For example, the spectral analysis of the pair of ODEs (2.2) for *x*>0, *t*>0 yields the formula (1.6).

We discuss this construction in §2*a*, but we stress that the basic idea is the derivation of a transform *simultaneously in x and in t* from the spectral analysis of an appropriate ODE pair, generalizing the derivation of the Fourier transform just outlined.

### (a) Boundary value problems for linear evolution partial differential equations on the half line

A powerful abstract result from the early 1970s, known as the Ehrenpreis–Palamodov fundamental principle (Ehrenpreis 1970), guarantees, at least under certain restrictions, the existence of a measure *ρ* and a contour *Γ* such that the solution of a linear boundary value problem for equation (1.3) can be expressed in the form(2.5)This form is similar to the form (1.4) of the solution of the initial value problem, but the abstract theorem does not yield the measure and contour explicitly. Our aim is to give a rigorous procedure to determine explicitly the measure *ρ* and the contour *Γ* in terms of the given data of the problem. To illustrate the difficulties that arise, consider the simplest third-order PDE posed on a half line(2.6)Firstly, it is necessary to determine the number of boundary conditions that must be prescribed in order to define a well-posed problem. It turns out that, in the case of the PDE (2.6), one initial condition at *t*=0 (the condition *q*(*x*, 0)=*q*_{0}(*x*)) and one boundary condition at *x*=0 must be prescribed. For the sake of concreteness, we assume that *q*(0, *t*)=*f*_{0}(*x*) is the prescribed (smooth) boundary condition.

The second problem becomes apparent when trying to solve this problem by using the Fourier transform in *x*. Then the ODE governing the time evolution is This ODE has a simple form, but the two functions *q*_{x}(0, *t*) and *q*_{xx}(0, *t*) are unknown. It turns out that there does not exist any classical transform in *x* that yields an *effective* solution for this problem (where effective means that only the *known* data of the problem are involved). The failure of a transform in *x* such as the Fourier, sine or cosine transforms to yield a solution in the form (2.5) is due to the fact that the associated eigenvalue problem, obtained by the classical separation of variables, is not self-adjoint and does not admit a self-adjoint extension. It is possible to use a transform in *t* such as the Laplace transform, but since the ODE for this transform is of the same order as the PDE, this approach becomes increasingly complicated as the order of the PDE increases, and since *t*∈[0,∞], it imposes unnecessary growth restrictions on the given boundary data. In contrast, the approach proposed by Fokas yields the solution in the form (2.5), with an algorithmic prescription to find *ρ*(*k*) and *Γ*.

The starting point of the new approach is to consider the two ODEs (2.1), rather than the PDE and the classical separation of variables. This formulation involves the two main ideas that make this transform method a genuine generalization of the Fourier transform, namely

Integration is not restricted to the real line, but it is allowed to take place on appropriate

**complex contours**, uniquely determined by the problem.The Lax pair formulation naturally yields a crucial algebraic relation, involving all the initial and boundary values, called the

**global relation**.

The conceptual steps of the Fokas' method can be summarized in the following scheme.

(i) Write the given PDE as the compatibility condition of the two ODEs in the Lax pair (2.1). For equation (2.6), this pair is(2.7)

(ii) The spectral analysis of the pair of ODEs (2.1) yields a representation for the solution of the boundary value problem of the form

(2.8)where ∂*D*^{+} is the boundary of the domain in the upper half plane where Re *ω*(*k*)≤0, and the function is defined in terms of *all* boundary values at *x*=0 (some of which are not known)(2.9)The representation (2.8) is obtained by following exactly the steps used to derive the Fourier transform via the spectral analysis of the ODE (2.1*a*), with the associated Riemann–Hilbert problem determined by using functions *μ*_{j}(*x*, *t*, *k*) that are solutions simultaneously of *both* ODEs in the Lax pair (2.1).

For the PDE (2.6), *ω*(*k*)=−i*k*^{3} and the domain *D*^{+} is shown in figure 2 and it is given by(2.10)while the function is(2.11)

Equations (2.1) also yield the global relation

(2.12)where is the Fourier transform of *q*(*x*, 0) and the function is given by equation (2.11). The relation (2.12) is valid for *k*∈^{−}. For the PDE (2.6), the global relation is the equationThis equation involves also the function , which is not known. However, the terms involving this unknown function can be ignored. Indeed, it can be shown, using the analyticity considerations in *D*^{+}, that they do not contribute to the final solution representation.

(iv) The final step is the use of the transformations that leave the dispersion relation *ω*(*k*) of the PDE invariant. This can be done in general, but we illustrate here this analysis for the PDE (2.6), when the given boundary condition is *q*(0, *t*)=*f*_{0}(*t*).

The global relation is one equation, while there are two unknown functions appearing in equation (2.8) through the definition (2.11) of , namely *q*_{x}(0, *t*) and *q*_{xx}(0, *t*). Moreover, the above relation is valid in ^{−} and not in the domain *D*^{+}, where we need to express the function . The last step uses the two transformations *λ*_{1}(*k*) and *λ*_{2}(*k*) arising as roots of *ω*(*λ*)=*ω*(*k*), *λ*≠*k*. In this example, these transformations are the roots of *λ*^{3}=*k*^{3}, hence they are simply rotations by 2*π*/3 and 4*π*/3 of *k*: *λ*_{1}=*αk*, *λ*_{2}=*α*^{2}*k*, *α*=e^{2πi/3}. In particular, for *k*∈*D*^{+}, *λ*_{1}(*k*), *λ*_{2}(*k*)∈*D*^{−}⊂ (figure 2). These transformation leave the functions invariant (*j*=0,1,2), so that we can evaluate the global relation at these two values and obtain a system of two equations for the two unknown functions and . The solution of this system yields for the spectral function in (2.8) the expression(2.13)Although this expression seemingly involves the unknown value of the solution at the final time *T* (through the Fourier transform ), we already mentioned that the unknown terms can be ignored. Indeed, in this case, the termis *bounded and analytic* in *D*^{+}, hence by Cauchy's theorem its contribution to the representation formula (2.8) vanishes.

*This final step can only be carried out in the complex plane*, and it would not be successful if the spectral functions where only defined for real *k*, as when using Fourier or similar real transforms. Indeed, the crucial fact is that *exactly two of the transformations leaving ω*(*k*) *invariant map the domain D*^{+} *to* ^{−}, *where the global relation is valid*. Hence from the global relation, one can derive a system of two equations for the two unknowns. In addition, the unknown function can only be eliminated using analyticity considerations.

In summary, we have developed the steps necessary to solve the problem in a way analogous to equation (1.4):

Direct map

Global relation + invariance in →explicit determination of , for *k*∈*D*^{+},

Inverse map

The method of solution for the PDE , *σ*=1,*i*, is analogous to the one used for the third-order example we treated above. The domain *D*^{+} is, in general, a union of simply connected regions such as (2.10), with angle *π*/*n*. The invariant transformations *λ*_{1}(*k*), …, *λ*_{n−1}(*k*) are the rotations by 2*π*/*n*, and the number of these transformations mapping each component of *D*^{+} into ^{−} matches exactly the number of unknown boundary values. It turns out that this is asymptotically the case for a generic PDE of order *n* of the form (1.3), hence the elements essential for the proof of the general case are contained in the example we just discussed.

### (b) Two-point boundary value problems for linear evolution partial differential equations

A construction entirely analogous to the one outlined in §2*a* can be used to derive an integral representation of the solution of a two-point boundary value problem for any linear evolution PDE. However, it is well known that in some cases, e.g. for a second-order problem, the solution can be represented as a *Fourier series* over the (infinitely many) eigenvalues of the *x*-differential operator, rather than as an integral. Hence, our construction must allow for this difference. In this section, we illustrate how this difference is expressed in terms of the global relation and of the general integral representation (2.5).

To illustrate the solution method for these types of problems, we consider the PDE (2.6), posed on the domain{0<*x*<*L*}×{0<*t*<*T*}. To define the problem, boundary conditions have to be prescribed at both *x*=0 and *L*. Indeed, it turns out that to obtain a well-posed problem, one needs to prescribe one boundary condition at *x*=0 and two boundary conditions at *x*=*L* (Pelloni 2004). We assume, to fix ideas, that *q*(0, *t*)=*f*_{0}(*t*), *q*(*L*, *t*)=*g*_{0}(*t*) and *q*_{x}(*L*, *t*)=*g*_{1}(*t*) are prescribed, smooth functions which are compatible with *q*_{0}(*x*) at *x*=0 and *L*.

Performing the step (ii) of the scheme outlined in §2*a*, we obtain an integral representation of the solution in the form(2.14)where *D*^{±}=*D*∩^{±}, and *D* is the region shaded in figure 2 and given by(2.15)The function is defined by equation (2.11), and has a similar definition(2.16)Therefore, in the representation (2.14), three terms are unknown, namely the terms involving *q*_{x}(0, *t*), *q*_{xx}(0, *t*), *q*_{xx}(0, *t*) and *q*_{xx}(*L*, *t*). However, the spectral functions and also satisfy the global relation, which is now valid in the entire complex plane and, in this case, is given by(2.17)Recall that *k*^{3}=(*α*^{j}*k*)^{3}, where *α*=e^{2πi/3} and *j*=1, 2. Evaluating the global relation at *k*, *αk* and *α*^{2}*k*, and using the invariance with respect to these rotations, we obtain the system of three equations(2.18)where and . This system can be solved for the three unknown functions (as in the case of the problem posed on a half line, analyticity considerations imply that the contribution of the unknown terms involving can be neglected, and we have dropped them from equation (2.18)). The solutions of this system are of the formwhere is the determinant of the system (2.18), and the functions *ρ*^{+}(*k*) and *ρ*^{−}(*k*) are explicitly expressed in terms of *q*_{0}(*x*) *f*_{0}(*t*), *g*_{0}(*t*) and *g*_{1}(*t*). For example,whereA similar expression holds for *ρ*^{−}(*k*).

The difference between this system and the one obtained in the case of boundary value problem posed on the half line is that the determinant *Δ*(*k*) of this system *has infinitely many complex zeros k*_{n}. In correspondence with these complex zeros, the solution of this system has poles. Hence, if the zeros of this determinant are inside the closed domains *D*^{+} or *D*^{−}, then the contour integrals must include the corresponding residue contributions. In such cases, e.g. for second-order problems, the integrals can be computed entirely in terms of residues and the integral representation yields the series solution. However, for the third-order example we are considering, the zeros of the determinant *Δ*(*k*) can be shown to lie precisely on the lines **L**_{1}, **L**_{2} and **L**_{3} shown in figure 2 (this is true for any choice of *uncoupled* boundary conditions; Pelloni 2005). These lines are outside the domains *D*^{±}, and it follows that the solution formula does not involve a series contribution. In contrast, if the given boundary conditions are coupled, for example, if the given conditions are(2.19)then the zeros of the determinant of equation (2.18) are precisely on the six half lines comprising ∂*D*^{±}, and the integral representation reduces to an infinite series.

To summarize, for third-order problems, when the zeros of the determinant function *Δ*(*k*) are inside the domain *D*, the integral representation is equivalent to a series one, but if these zeros are outside *D*, then in general the solution does not involve a series. (We remark for second-order problems, it is *always possible* to write the integral representation purely as an infinite series.)

We stress that the integral representation can *always be* derived, regardless of the existence of a series representation for the solution, and it presents important advantages. For example, in contrast with the series one, it is *uniformly convergent* both at *x*=0 and at *x*=*L*. It is also more convenient for the study of the long-time asymptotic behaviour of the solution. Furthermore, it is remarkably well suited to numerical approximations, as we discuss later. It should also be noted that the integral representation of the solution does not require a detailed analysis of the determinant of the system (2.18). This is to be contrasted with the fact that, since the zeros of this determinant define the discrete spectrum of the associated *x*-differential operator, their determination is crucial for the construction of the classical eigenfunction series representation.

### (c) Numerical computations

One of the main advantages offered by the complex integral representation of the solution is that one can use the analyticity properties of the integrand to deform the contour of integration in a way suited to speed up numerical convergence and achieve greater accuracy. For example, for problems posed on the half line, part of the contour of integration involving the boundary condition (namely the part of ∂*D*^{±} which is different from the real axis) can always be deformed in such a way that the integrand possesses *fast decay*. This is to be contrasted with the evaluation of oscillatory exponential integrals, such as Fourier type integrals on the real line, which poses practical difficulties.

The above remarks apply also to the case of two-point boundary value problems. In this case, the contour of integration can always be deformed in such a way that the exponential kernel to be computed is *decreasing* and moreover avoids the spectrum of the *x* operator, at which the integrand is singular. This is to be contrasted with the computation of the eigenfunction series, which is known to be unworkably slow in some cases.

It should be emphasized that for the solution of simple boundary value problems like those discussed in §§2*a*,*b*, it is possible to avoid the simultaneous spectral analysis and to derive directly from the global relation the relevant integral representations, like those in equations (2.8) and (2.14). However, we have chosen to present the simultaneous spectral analysis, since this is the only approach that generalizes to integrable nonlinear PDEs.

## 3. Analytic integral inversion formulae and boundary value problems for linear evolution partial differential equations on time-dependent domains

In this section, we illustrate how the approach of Fokas can be used to derive an inversion formula for more general integral transforms. The results are based on the spectral analysis of appropriate ODEs. The first significant application of this method appeared in the work of R. Novikov, who was able to invert the attenuated Radon transform by using a simple extension of a novel derivation of the Radon transform (Novikov 2002). This analytic inversion has important applications in medical imaging, namely in the technique known as single particle emission computerized tomography.

Here, we consider an application of this method to a problem that arises naturally when considering a boundary value problem for the general linear evolution PDE (1.3) posed on a time-dependent domain of the formwhere we assume that *l*(*t*) is a given smooth convex function.

Consider, for example, the PDE i*q*_{t}+*q*_{xx}=0 posed on * B*, with

*q*(

*x*,

*t*) and

*q*(

*l*(

*t*),

*t*) given initial and boundary conditions. This PDE is the linear limit of the equation (4.2), an integrable equation which we discuss in §4. To solve this boundary value problem, one needs to characterize the unknown boundary value

*q*

_{x}(

*l*(

*t*),

*t*), and this requires the inversion of the integral(3.1)with

*f*(

*t*)=

*q*

_{x}(

*l*(

*t*),

*t*). The expression (3.1) involves a simple variation of the usual Fourier transform (with

*k*replaced by

*k*

^{2}). However, it appears that such transforms have not been considered previously in the mathematics literature.

The inversion of the transform (3.1), i.e. the reconstruction of the function *f*(*t*) from the knowledge of the transform *F*(*k*), is achieved by performing the spectral analysis of the ODE(3.2)What distinguishes the spectral analysis of equation (3.2) from the analysis of an equation such as equation (2.1*a*) is that, in the present case, it is not possible to find a function *μ*(*x*, *t*, *k*) which is sectionally analytic in . As a consequence, the inversion formula arises from the solution of a generalization of the Reimann–Hilbert problem, known as the *d*-*bar problem*. This implies that the inversion cannot be given by a closed formula, but it is obtained through the solution of a linear Volterra integral equation; namely, we obtain the following result (Fokas & Pelloni 2006).

*Let F(k) be defined in terms of f(t) by equation* *(3.1)*, *where f(t) is a sufficiently smooth function. Then f(t) satisfies the Volterra integral equation*(3.3)*where J(s, t) is the well-defined kernel**with* , *and* *where*(3.4)*is shown in* *figure 3*.

The derivation of the kernel *J*(*s*, *t*) involved in the Volterra integral equation (3.3) is given in Fokas & Pelloni (2006). Here, we only note that this kernel is defined as an integral with respect to the complex variable *k*, and that with respect to this variable, its integrand possesses *fast decay*. As we have already stressed, this fast decay property is important from the point of view of efficient numerical computations of the transform.

The transform (3.1) plays the role of the Fourier transform for the boundary value problem(3.5)Indeed, the solution of this boundary value problem requires the determination of the boundary value *q*_{x}(*l*(*t*), *t*). The global relation for this boundary value problem isThis equation, together with the equation (3.3), yields *q*_{x}(*l*(*t*), *t*) as the solution of a Volterra linear integral equation of the form(3.6)wherewhere *A* is a numerical constant and the kernel *J*(*s*, *t*) defined in the statement of theorem 3.1.

## 4. Integrable evolution partial differential equations

We now turn to the generalization of the Fourier transform to a nonlinear setting. This generalization first appeared in the famous paper (Gardner *et al*. 1967), where the initial value problem for the KdV equation(4.1)was solved analytically by what is now known as the inverse scattering transform. This transform is in essence a nonlinearization of the Fourier transform in *x* that yields the solution of the corresponding linearized problem (2.6).

Soon after this result was obtained, it was realized that the only essential property of the PDE needed in order to apply this methodology is that the PDE possess a certain kind of separability, namely, that it can be written as the compatibility condition of a pair of linear ODEs (Lax 1968). This pair is now known as the *associated Lax pair*. The existence of a Lax pair can be taken as the defining property of an integrable PDE. Several important nonlinear integrable PDEs of mathematical physics, such as the nonlinear Schrödinger equation and the sine-Gordon equation, were solved by this approach. To solve the initial value problem, this nonlinear transform is considered as a transform in *x*, just as the Fourier transform in *x* is applied to a linear evolution PDE after separation of variables. The ODE in *x* of the Lax pair is used to derive the transform, and then the ODE in *t* is used to compute the time evolution of the ‘scattering data’. However, it turns out that in order to solve boundary value problems, one needs to consider *both ODEs in the Lax pair simultaneously*, in analogy with the approach to solve boundary value problems for linear PDEs discussed in §2. The generalization to boundary value problems of the inverse scattering methodology had to wait over 30 years until the work of Fokas clarified this crucial difference (Fokas 2002).

We describe below how this approach is applied to the solution of both the initial value and the boundary value problem for the particular example of a very important integrable PDE, the defocusing nonlinear Schrödinger (NLS) equation(4.2)For brevity, we do not include the focusing case, which supports soliton solutions, though their inclusion would not pose any particular difficulty.

This equation is the compatibility condition of the Lax pair(4.3)(4.4)where *Ψ*(*x*, *t*, *k*) is a 2×2 matrix, and

### (a) Initial value problem

We first pose equation (4.2) on the real line and prescribe a decaying initial condition *q*(*x*, 0)=*q*_{0}(*x*)∈().

**Direct map**

*q*_{0}(*x*)→{*a*(*k*), *b*(*k*)}, with(4.5)where (*Φ*_{1}, *Φ*_{2}) is defined in terms of *q*_{0}(*x*) through the solution of the Volterra integral equation(4.6)where *x*∈, Im(*k*)≥0. The vector (*Φ*_{1}, *Φ*_{2}) is a particular solution of the ODE (4.3) with *q*(*x*, *t*) replaced by *q*_{0}(*x*).

**Inverse map**(4.7)where and *M*(*x*, *t*, *k*) is defined in terms of *γ*(*k*) through the solution of the following matrix problem:

*M*is analytic in \*M*=*I*+*O*(1/*k*) as*k*→∞On ,

*M*satisfies the jump condition

The equations (2.2) and (2.4) provide an analogue of the direct and inverse Fourier transform.

### (b) Boundary value problems

We now consider the equation (4.2) posed on the half line *x*>0. We assume that one decaying initial condition *q*_{0}(*x*) and one boundary condition *q*(0, *t*)=*f*_{0}(*t*) are prescribed, and that the functions *f*_{0}(*t*) and *q*_{0}(*x*) are smooth and compatible (Fokas *et al*. 2005).

The solution scheme, obtained by combining the generalization and ideas illustrated in §§2 and 3, is the following.

**Direct map**

{*q*_{0}(*x*), *q*(0, *t*), *q*_{x}(0, *t*)}→{*a*(*k*), *b*(*k*); *A*(*k*), *B*(*k*)}, with *a*(*k*), *b*(*k*) defined in terms of the solution (*Φ*_{1}, *Φ*_{2}) of the linear Volterra integral equations (4.6) (but with −∞ replaced by 0), and *A*(*k*), *B*(*k*) defined in terms of *f*_{0}(*t*) and *q*_{x}(0, *t*) through the solution (*Ψ*_{1}, *Ψ*_{2}) of the linear Volterra integral equations(4.8)The vector (*Ψ*_{1}, *Ψ*_{2}) is a particular solution of the ODE (4.4) with *q*(*x*, *t*) and *q*_{x}(*x*, *t*) replaced by *f*_{0}(*t*) and *q*_{x}(0, *t*).

**Inverse map**(4.9)where , *Γ*(*k*) is explicitly defined in terms of *a*(*k*), *b*(*k*), *A*(*k*), *B*(*k*), and *M*(*x*, *t*, *k*) is explicitly defined in terms of *γ*(*k*) and *Γ*(*k*) through the solution of a matrix problem similar to the one associated with the initial value problem, but with jumps on both the real and the imaginary axes.

As in the linear case, in order to express the solution effectively it is necessary to analyse the *global relation* to determine the unknown boundary conditions, in this case *q*_{x}(0, *t*)=0. The global relation isThis procedure now yields the characterization of this unknown boundary value through the solution of a pair of nonlinear ODEs. For example, if *q*_{0}(*x*)=0, the global relation is explicitlyIt is possible to solve this nonlinear integral equation *analytically* for the function *q*_{x}(0, *t*) in terms of *f*_{0}(*t*), and the vector (*Ψ*_{1}, *Ψ*_{2}). Substituting the resulting expression in equation (4.8), one obtains two nonlinear ODEs for *Ψ*_{1} and *Ψ*_{2}.

We stress that in spite of the considerable technical difficulties involved, in particular, in the last step above, the method is conceptually the analogue of the method described in §2.

## 5. Conclusions

We have discussed how the Fokas transform method for solving boundary value problems for linear and integrable nonlinear PDEs can be viewed as an extension of the Fourier transform method, and indeed, how in simple cases it reduces to the Fourier transform. The unifying character of the steps involved in the Fokas method makes it attractive from the theoretical and formal point of view. For nonlinear integrable problems, this approach is, at present, the only existing method yielding results in a general context.

From the point of view of numerical computations, the most important feature of this method is the decay of the integrands appearing in the integral representations in the complex *k* plane.

The Fokas method has a much broader domain of applicability than it is possible to present here. For example, elliptic problems can also be treated by this general approach. In this case, the analysis of the global relation, which is the crucial step in the methodology, may involve the solution of additional Reimann–Hilbert problems. As regards numerical approximations, preliminary results indicate that, using this approach, the Dirichlet to Neumann map for linear elliptic problems can also be evaluated, at least for some important examples, with exponential accuracy.

We mention also that the global relation appears to play a fundamental role in other contexts. For example, it appears in several problems in applied mathematics including water waves theory (Ablowitz *et al*. 2006) and Hele-Shaw type problems (Crowdy 2002).

Current research directions include implementation of the new method for the solution of the basic linear elliptic PDEs not only in polygonal domains, but also in circular and elliptic ones, as well as the solution of the basic elliptic PDEs in three dimensions in polyhedral, spherical and ellipsoidal domains.

Work in progress also concerns the extensions of this transform method to nonlinear integrable evolution equations in two space dimensions, such as the so-called KP and DS equations.

## Acknowledgements

I gratefully acknowledge the support of the Leverhulme Trust through the award of a Research Fellowship.

## Footnotes

One contribution of 23 to a Triennial Issue ‘Mathematics and physics’.

- © 2006 The Royal Society