Approximate Analytical Solutions to the Heat and Stokes Equations on the Half-Line Obtained by Fokas’ Transform ()
1. Introduction
Fokas’ method (Fokas, 1997 [1] and 2002 [2] ) 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
(1.1)
and
(1.2)
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, [3] ), Fokas et al. (2009, [4] ), and Papatheodorou-Kandili (2009, [5] ) 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 [3] 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 [6] ). Hence, our numerical solution is more accurate than Fokas’ solution and even than that of Papatheodorou-Kandili (2009, [5] ), 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:
(2.1)
For the Stokes’ equation of the first kind with the same initial and boundary conditions as above the solution is:
(2.2)
To evaluate numerically these integrals, Flyer-Fokas (2008, [3] ) and Fokas et al., (2009, [4] ) 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:
(2.3)
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):
(2.4)
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
.
(2.5)
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:
(3.1)
Proof. First, we convert the integral along the infinite contour to an integral on the real line. From (2.1), we have
(3.2)
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:
(3.3)
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
,
and,
(3.4)
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:
(3.5)
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:
(3.6)
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:
(3.7)
For the third integral in (3.4), we use the partial fractions decomposition of
as:
(3.8)
We get:
(3.9)
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:
(3.10)
where
.
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
(3.11)
• 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:
(3.12)
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
(3.13)
Therefore Equation (3.1) becomes:
(3.14)
3.2. The Stokes Equation of the First Kind
Lemma 3.2. The solution of Equation (1.2) is given by:
(3.15)
where
,
;
,
;
,
.
Proof. From Equation (2.2), we have that
(3.16)
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
Therefore,
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:
(3.17)
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
Therefore,
For the third integral, we use the partial fractions decomposition of
as follows:
Therefore,
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:
and
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 [7]; 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 [7].
Writing the Cauchy kernel as a combination of the Poisson and Hilbert kernels, we get:
(3.18)
Therefore:
(3.19)
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
(3.20)
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
(3.21)
Here we used the fact that
and
, so
.
• If
, then
and
are in
; and
and
are in
. Therefore,
and
(3.22)
• If
, then all the roots
are in
. Therefore,
and
(3.23)
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:
(3.24)
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:
(3.25)
where
are computed later.
Proof. By analogy with (3.18) and following Constandinides-Marhefka [8], let the incomplete Airy functions be given by:
(3.26)
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, [9] ). Also see Michaeli (1996, [10] ) 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
)
(3.27)
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
.
Proof.
• 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.
Therefore:
(3.28)
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:
(3.29)
Therefore, the solution for
,
, and
:
(3.30)
Using the formulas for
,
,
,
,
,
and the Wronskian
, we get after simplification:
(3.31)
Now, we have:
Thus,
and finally
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 [7], that the imaginary error function
has the following asymptotic expansion:
(4.1)
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
Therefore,
(4.2)
where
(4.3)
Replacing (4.1) in (3.20), we formally get the asymptotic expansion of the solution
(4.4)
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
,
:
More precisely,
(4.5)
where the error is bounded by
(4.6)
which converges to 0 for fixed N as
.
Proof. The fraction
can be written as:
Therefore,
(4.7)
where
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,
(4.8)
Therefore
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:
(4.9)
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.
5. Discussion
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)
(5.1)
The leading term
and the error are both of size
.
The next term in this expansion is (after simplification)
(5.2)
Explicitly, as
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
Acknowledgements
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 [3] )