Fokas’ method (Fokas, 1997  and 2002  ) is a recently introduced approach that allows solving a large class of PDEs. Also known as the unified transform, this method extends classical approaches, such as separation of variables or the scattering transform, contains them as special cases, and gives solutions in situations where the classical methods cannot.
Classical methods work only poorly for spatial domains with boundary, when dealing with inhomogeneous boundary conditions. In this context, the key contribution of Fokas is the so-called “global relation”, which combines specified and unknown values of the solution and its derivatives on the boundary.
Fokas’ method gives the solutions as integrals on an unbounded contour in the complex plane. For example, for the heat and Stokes (first kind) equations on the half-line
respectively, this method constructs the solutions as integrals (2.1) and (2.2) in the complex plane involving an x-transform of the initial condition and a t-transform of the boundary condition .
To evaluate these integrals numerically, Flyer-Fokas (2008,  ), Fokas et al. (2009,  ), and Papatheodorou-Kandili (2009,  ) used the steepest descent method, deforming the contour in the complex plane, then used the simple trapezoid rule, after parametrizing the curve, to compute the integral.
Following Flyer-Fokas  and for the sake of enabling a comparison, we make the choice
But, any reasonable (e.g. tempered distribution) Dirichlet boundary condition can be decomposed into a superposition of sines, using the sine Fourier transform. Thus, we can solve the above equations for general boundary data. Initial data are simpler, to the point that one does not need Fokas’ method to handle them and will be treated elsewhere.
In this paper, we evaluate the integrals (2.1) and (2.2) by an alternative approach. First, we obtain a closed formula for the solution of the heat Equation (1.1) in terms of elementary functions and one special function, namely the imaginary error function erfi:
where and are defined later. Second, we obtain an exact expression for the solution of Stokes’ Equation (1.2) in terms of elementary functions and another special function, the incomplete Airy functions:
where are functions of the incomplete Airy functions and will be determined later. The quantities and .
The above formulas lend themselves well to applications, since there exist fast and highly accurate methods for computing the imaginary error function and incomplete Airy functions. For example, the imaginary error function is a standard function in MATLAB, where it is estimated with 10−20 accuracy by means of Padé approximants (Cody, 1969  ). Hence, our numerical solution is more accurate than Fokas’ solution and even than that of Papatheodorou-Kandili (2009,  ), who obtained a solution with 10−15 accuracy for the heat equation on a finite spatial domain. We also obtained speed improvements of at least an order of magnitude, see the Appendix A.
A related problem is that, as stated in Flyer-Fokas: “for , the integration interval is infinite and any truncation will result in the integral not converging exactly to the boundary condition”. By contrast, our exact expressions for the solutions extend all the way to the lateral boundary without this issue of convergence (they are still ill-behaved as , though).
In this paper, we also derive the asymptotic expansion for the heat equation solution with precise bounds for the error term, which allows one to compute the solution with arbitrarily high precision. This only requires evaluating elementary functions and the imaginary error function erfi.
The paper is organized as follows. In the next section, we briefly describe Fokas’ method as applied to Equations (1.1) and (1.2). In Section 3 we derive the exact formulas for the solutions. In Section 4 we obtain an asymptotic expansion for the heat equation solution, together with an estimate for the error size. Finally, in the Appendix B, we describe the numerical implementation of our scheme and compare the results to those of Flyer-Fokas.
2. Fokas’ Integral Solutions
In this section, we present the solutions of Equations (1.1) and (1.2) given by Fokas and that take the form of integrals. For the heat equation on the half-line with initial condition ( ) and a Dirichlet boundary condition ( ), the solution is:
For the Stokes’ equation of the first kind with the same initial and boundary conditions as above the solution is:
To evaluate numerically these integrals, Flyer-Fokas (2008,  ) and Fokas et al., (2009,  ) deformed the contour of integration to a path in the region where the integrand decays exponentially fast for large k. Specifically, in order to get rapid convergence of the numerical scheme, they deformed the contour of integration to a hyperbola above the real axis and asymptotically below . After that, they mapped the hyperbola from the complex plane to the real line using the following parametrisation:
which maps points on the real line to points on a hyperbola in the complex plane. The integral in Equation (2.1) becomes (also see Formula (3.8) in Flyer-Fokas):
The parameters and were set to 1. After mapping the integral to the real line, Flyer-Fokas truncated the infinite integration interval to a finite one and used the simple trapezoidal rule. The same deformation-mapping procedure was used with the stokes Equation (2.2) except that was replaced with .
3. Exact Formulas for the Solutions
3.1. The Heat Equation
Using several variables and contour changes, as well as Cauchy’s residue theorem, we obtain a more manageable expression for the solution of the heat Equation (1.1). Our starting point is identity (2.1). Define the following four quantities , , which will appear in the computation:
More precisely, let
They have the property that and . Consequently, and , We also note for future reference that
Lemma 3.1. The solution of Equation (1.1) is given by:
Proof. First, we convert the integral along the infinite contour to an integral on the real line. From (2.1), we have
We consider the contour , where is the part of C on the boundary of the domain and and are circular arcs of radius R. Let denote the integrand. On the contour C there are two removable singularities at and at . There are no poles either inside or on the contour. Using the analyticity of the integrand, we have:
The contributions of the integrals along and vanish as , by Jordan’s lemma. Equation (3.1) can be written as:
Along , , for and . The first integral in (3.3) becomes
We will show that the integral of the modulus of the integrand converges to 0 as and therefore the integral of the integrand converges to 0 as well.
Therefore, when , the integral converges to 0 since for .
For the second integral in (3.3), we have
which also converges to 0 as .
For the third integral, we have
which converges to 0 as well as since and for .
Therefore, the three integrals vanish also along since and for . Also, as , the contour becomes . Therefore, as ,
Now, we evaluate each of these three integrals. The first and second one are computed directly using the residue theorem. The third one is computed using a different strategy.
For the first integral, we consider a contour that is the boundary of an upper semicircular region of a circle of a radius large enough to include the pole . Using the residue theorem and Jordan’s lemma, we get:
In this integral, we need not consider the residue at the other pole since it is not inside the contour. Same holds for the second integral. We consider a contour that is the boundary of an upper semicircular region of a circle of radius large enough to include the pole . Using the residue theorem, we get:
In this integral as well, we need not consider the residue at the other pole since it is not inside the contour. Summing (3.5) and (3.6) gives:
For the third integral in (3.4), we use the partial fractions decomposition of as:
We next compute these four integrals. After changing the variable k by and completing the square in , we have:
• The first integral
Another change of variable: , gives:
The integrand in the right hand side of (3.10) has as a pole. Deforming the integral back to the real line will result in a residue at exactly when
is non-positive (Cauchy’s theorem). If , then there is a residue. If , then the pole is not inside the contour. Therefore, the first integral is:
where the integral is interpreted in the principal value sense when and
• The second integral
where . Note that in this case there is no residue
associated with the pole since it is outside the contour when deforming the integral to the real line.
• The third integral
where . Note also that deforming the integral on the
real line in this case will not result in a residue since the pole is never inside the contour.
• The fourth integral
where . Deforming the integral on the real line will
result in a residue depending on the location of the pole . The residue at is:
The difference of the residues at and is
which cancels in the boundary case . Therefore, the sum of (3.7) and (3.9) produces (3.1), which is what we wanted to prove.
We note that
is odd and . Since and , it follows that and . Then
Therefore Equation (3.1) becomes:
3.2. The Stokes Equation of the First Kind
Lemma 3.2. The solution of Equation (1.2) is given by:
where , ; , ; , .
Proof. From Equation (2.2), we have that
We convert the integral along the infinite contour to an integral on the real line. Consider the contour , where is the part of C on the boundary of the domain and and are circular arcs of radius R. Let denote the integrand. On the contour C, there are
four poles: , , and . Using the analyticity of the integrand, we obtain:
because the residues at , , and are all equal to 0. (It is easy to show that. The proof is available upon request).
The contributions of the integral vanish along and . In fact, can be written as:
Along , , for and . Therefore:
• The first integral becomes
when , the integral converges to 0 since for .
• For the second integral, we have:
which converges to 0 as .
• The third integral
which converges to 0 as since and for .
Along , the three integrals vanish also since and for .
Also, as , the contour becomes .
Therefore, as , .
Therefore, the contour is deformed to the real axis:
The first and second integrals are computed directly using the residue theorem. For the first one, the roots of the denominator of the integrand are:
, and . We therefore consider
a contour that is the boundary of an upper semicircular region of a circle that has a radius large enough to include the poles and . We do not need to consider the residue at the other pole since it is not in the considered contour. Therefore:
where (after some computations)
Substitution of and in the integral gives
For the second integral, the roots of the denominator of the integrand are: , and . We therefore consider a
contour that is the boundary of an upper semicircular region of a circle that has a radius large enough to include the poles and . We do not need to consider the residue at the other pole since it is not in the considered contour. Therefore:
where (after some computations),
Substitution of and in gives
For the third integral, we use the partial fractions decomposition of as follows:
In light of identities (3.1) and (3.15), to obtain exact formulas for the solutions of Equations (1.1) and (1.2) we are left with evaluating integrals of the form:
These Cauchy integrals can be computed by means of the Hilbert transform. When restricted to the positive half-line , the first integral is known as the Goodwin-Staton integral, see Abramowitz and Stegun ; also see Dawson’s integral. We compute both integrals in the following section in terms of special functions.
3.3. The Hilbert Transform
Lemma 3.3. Consider two entire functions f and g, such that is bounded on and is bounded on , then . In particular, consider two entire functions f and g such that f and g are real-valued on the real line and is bounded on . Then .
Proof. Since is bounded on , its Fourier support is contained in , (that is its Fourier transform vanishes over ). Therefore, for , . Likewise, the boundedness of implies that for . Hence, for , , and therefore . For , under our boundedness assumptions, one cannot get anything worse than a constant, and the Hilbert transform is only unique up to a constant.
For the second conclusion, by the Schwartz reflection principle and , so . Hence the boundedness of in the upper half-plane implies that of in the lower half-plane and we can fall back on the previous argument.
Regarding the error function, since and are real-valued on the real axis and is bounded in the upper half-plane, it follows that the Hilbert transform of is
which is also an analytic function, and real-valued on the real line. Here the imaginary error function is defined as and is the usual error function:
For , is real-valued too, as seen from the fact that
This Hilbert transform is closely related to Dawson’s function and to the Mittag-Leffler functions : , see .
Writing the Cauchy kernel as a combination of the Poisson and Hilbert kernels, we get:
where denotes the Cauchy principal value, the Hilbert transform, the imaginary error function, the upper half complex plane, and the lower half complex plane.
The expression (3.19) is analytic, for and for separately, and is bounded on the complex plane, but is not continuous, due to the jump
discontinuity across the real axis. It decays at a rate of as . At this stage, we obtain an exact formula for the solution of the heat equation.
Proposition 3.4. The exact solution of Equation (1.1) is given by
The above formula is valid for all with , (since are not defined when ). This agrees with the fact that, although the expression (3.1) defines a discontinuous function, the solution is continuous up to the boundary and smooth for .
Proof. Of the four values , and are always in the upper half-plane, while and can be in , , or . Thus, we have to distinguish the following three cases:
• If , then and are in and and are in . Therefore,
and from Equation (3.1) above
Here we used the fact that and , so .
• If , then and are in ; and and are in . Therefore,
• If , then all the roots are in . Therefore,
The Formulas (3.21) and (3.23), which correspond respectively to the regions and , are the same and reduce to the Formula (3.22) when x approaches from both regions. In all three cases, we have proved (3.20).
For arbitrary boundary data , this leads to the following solution of the heat equation:
where is the sine transform of the boundary condition :
Next, we compute the integral and obtain an exact expression for the solution of the Stokes equation of the first kind. The change of variable in the expression of gives:
where and . This is the Cauchy integral of , which is an analytic function bounded on the real line—so, for fixed x and , .
Proposition 3.5. The exact solution of Equation (1.2) is given by:
where are computed later.
Proof. By analogy with (3.18) and following Constandinides-Marhefka , let the incomplete Airy functions be given by:
For the integral to converge, the upper integration limit in (3.26) can be taken to be , for . Over each range, the integral is constant.
Due to the integrand’s rapid decay, the integral (3.26) converges for any , defining an analytic function of both. is known as the incomplete Airy integral, introduced by Levey-Felsen (1969,  ). Also see Michaeli (1996,  ) for a so-called “Airy-Fresnel integral”.
When , , , and solve the inhomogeneous Airy equation . More generally, for any they solve it with a source term of
The three incomplete Airy functions are then related, modulo the usual Airy functions and :
Since (by a change of variable )
for and . Hence for and .
By the Schwartz reflection principle, it follows that , so . In particular, for and .
For and , is purely imaginary:
When , is a positive real number. is an absolute constant that can be computed in terms of the function:
When , can also be expressed in terms of the incomplete function:
However, for and , the incomplete Airy functions cannot be reduced to a simpler expression.
Lemma 3.6. The function
is uniformly bounded for .
• If , then the integral converges and gives a continuous function of k and therefore it is uniformly bounded.
• If , then we have and the first part of the above function becomes:
which is bounded given that . The absolute value of the integral of the above function is also bounded. In fact:
Now, let and then implies and the above function becomes:
which is bounded for .
Now, let , where , then the above function will be expressed in terms of the two other incomplete Airy functions. In fact:
where , is bounded for in some half-plane. On the other half plane, the function becomes:
where , is also bounded.
Now Let’s apply the same change of variables to the derivative of the above function with respect to x. We have:
Then the change of variables and gives:
which is bounded. In terms of the two other incomplete Airy functions, we let , where . Therefore:
where: , and
where: , which are both bounded.
Now, we have the general formula of the Hilbert transform for the function :
Taking the Fourier transform and using the above results, we get:
Using Equation (3.28), the above system becomes:
Therefore, the solution for , , and :
Using the formulas for , , , , , and the Wronskian , we get after simplification:
Now, we have:
4. Asymptotic Expansion of the Heat Equation’s Solution
4.1. Deriving the Asymptotic Expansion
Our exact expression (3.20) for the solution q is not very transparent. In order to better understand its properties, we next derive an asymptotic expansion for . It is well known, see Abramowitz and Stegun , that the imaginary error function has the following asymptotic expansion:
Here is the semifactorial function, with the convention that .
For completeness, we rederive this formula below, together with some precise bounds for the error in this expansion. We start from the definition of the error function:
By repeated integration by parts, we obtain
where is the remainder after N terms and is given by
Replacing (4.1) in (3.20), we formally get the asymptotic expansion of the solution
However, this computation is slightly misleading, because it does not give a good estimate for the size of the error. If we use Formula (4.4), the error will be large in some cases of interest. In fact, Formula (3.1) is better in a certain sense than Formula (3.20), because the sine term in the former, which is absent in the latter, can sometimes be a leading term, see below.
4.2. Estimating the Error in the Asymptotic Expansion
A more straightforward approach also leads to a better error estimate. For this purpose, we discard the exact Formula (3.20), obtained by estimating the Cauchy integral using the Hilbert transform, and go back to (3.1).
Lemma 4.1. The following asymptotic expansion is valid as for , :
where the error is bounded by
which converges to 0 for fixed N as .
Proof. The fraction can be written as:
The values of the integrals in (4.7) are known: for odd they are zero because the integrands are odd. For even values of we obtain, by repeated integrations by parts,
An upper bound for the error term , based on the exact formula above, is
Substituting (4.7) in Formula (3.1) and taking into account (3.14), we get the following asymptotic expansion of the solution :
Proposition 4.2. The solution to Equation (1.1) admits the following asymptotic expansion:
where the error is bounded by
Proof. From (3.1) and (4.4) we get the following rough bound for the error :
The first line can also be used to obtain other bounds, e.g.
Thus, the asymptotic expansion (4.9) may be inaccurate near the diagonal , but covers both the case when x is fixed and t is large and the case when t is fixed and x is large.
Interestingly, the solution behaves differently in these two regions. In the large x regime, the leading order approximation to the solution is given by (after simplification)
The leading term and the error are both of size .
The next term in this expansion is (after simplification)
The second term is of size and the error is of size .
On the other hand, in the large t regime, the term
dominates (at least on average, since the sine can be zero), so as the first two terms are
Again, has the same size, , as the remainder. Considering one more term, we have
Here has a size of at most .
Finally, let us mention that along the diagonal a similar analysis implies that
M.B. is partially supported by a grant from the Simons Foundation (No. 429698, Marius Beceanu) and by startup funds provided by the University at Albany, SUNY.
Appendix A. A Comparison with Fokas’ Method
For comparison with Fokas’ method, we repeatedly ran both our code and Flyer-Fokas’ code on the same publicly available MATLAB/Octave online implementation http://www.tutorialspoint.com/execute_matlab_online.php, and averaged the times we obtained. The running times and averages are listed in TableA1.
We found that it takes much longer to produce the figure below for the heat equation using Fokas’ method and that the difference in running times becomes more pronounced for a bigger grid.
Table A1. Running time difference in seconds between Fokas' method and our method in solving the heat equation on a half-line.
Appendix B. Code
For reference, this is the MATLAB/Octave code we used to compare the running times of Fokas’ method and our method and to generate FigureB1.
Figure B1. The solutions of the heat equation displayed on and (left) and of Stokes equation displayed on and (right).
B.1. Our Matlab code for the heat equation
B.2. Flyer-Fokas’ code for the heat equation (Adapted from  )
 Flyer, N. and Fokas, A.S. (2008) A Hybrid Analytical-Numerical Method for Solving Evolution Partial Differential Equations in the Half-Line. Proceedings of the Royal Society A, 464, 1823-1849.
 Fokas, A.S., et al. (2009) A Semi-Analytical Numerical Method for Solving Evolution and Elliptic Partial Differential Equations. Journal of Computational and Applied Mathematics, 227, 59-74.
 Papatheodorou, T.S. and Kandili, A.N. (2009) Novel Numerical Techniques Based on Fokas Transforms for the Solution of Initial Boundary Value Problems. Journal of Computational and Applied Mathematics, 227, 75-82.
 Michaeli, A. (1996) Asymptotic Analysis of Edge-Excited Currents on a Convex Face of a Perfectly Conducting Wedge under Overlapping Penumbra Region Conditions. IEEE Transactions on Antennas and Propagation, 44, 97-101.