Approximate Analytical Solutions to the Heat and Stokes Equations on the Half-Line Obtained by Fokas’ Transform

In this paper, we evaluate the integrals that are solutions of the heat and Stokes’ equations obtained by Fokas’ transform method by deriving exact formulas. Our method is more accurate and efficient than the contour de-formation and parametrization used by Fokas to compute these integrals. In fact, for the heat equation, our solution is exact up to the imaginary error function and for the Stokes equation, our solution is exact up to the incomplete Airy function. In addition, our solutions extend to the lateral boundary without convergence issues, allow for asymptotic expansions, and are much faster than those obtained by other methods.


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.  respectively, this method constructs the solutions ( ) , q x t as integrals (2.1) and (2.2) in the complex plane involving an x-transform of the initial condition 0 q and a t-transform of the boundary condition 0 g . 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 ( ) 0 g t 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 ( )   (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 0 x = , 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 0 x = without this issue of convergence (they are still ill-behaved as 0 t → , 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.

Fokas' Integral Solutions
In this section, we present the solutions of Equations (1.1) and (1.2) given by 1 e e e e , e d 2 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, [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 D + ∂ . After that, they mapped the hyperbola from the complex plane to the real line using the following parametrisation:

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 ( ) , q x t of the heat Equation (1.1). Our starting point is identity (2.1). Define the following four quantities i z , 1 4 i ≤ ≤ , which will appear in the computation: We consider the contour 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 iλ . Using the residue theorem and Jordan's lemma, we get: In this integral, we need not consider the residue at the other pole iλ − 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 iλ − − . Using the residue theorem, we get: In this integral as well, we need not consider the residue at the other pole iλ − 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 4 2 k k λ + as: We get: We next compute these four integrals. After changing the variable k by k k t =  and completing the square in k  , we have:  The first integral The integrand in the right hand side of (3.10) has 1 z − as a pole. Deforming the integral back to the real line will result in a residue at there is a residue. If 2 x t λ < , then the pole is not inside the contour. Therefore, the first integral is: where the integral is interpreted in the principal value sense when 1 z ∈  and ( ) ( ) Note that in this case there is no residue associated with the pole 2 z − since it is outside the contour when deforming the integral to the real line.  The third integral . Note also that deforming the integral on the real line in this case will not result in a residue since the pole    lim e e e .
The difference of the residues at

Res
Res e e e 2 e sin , 2 which cancels in the boundary case 2 x t λ = . Therefore, the sum of (3.7) and

The Stokes Equation of the First Kind
We convert the integral along the infinite contour to an integral on the real line. Consider the contour Using the analyticity of the integrand, we obtain: show that. The proof is available upon request).
The contributions of the integral vanish along  For the second integral, we have: Also, as R → ∞ , the contour 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: . 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 1 k and 2 k . We do not need to consider the residue at the other pole 3 k since it is not in the considered contour. Therefore: where (after some computations) ( ) ( ) Res k in the integral gives For the second integral, the roots of the denominator of the integrand are: . 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 1 k and 2 k . We do not need to consider the residue at the other pole 3 k since it is not in the considered contour.
For the third integral, we use the partial fractions decomposition of   When restricted to the positive half-line [ ) 0, +∞ , 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. Lemma 3.3. Consider two entire functions f and g, such that f ig + is bounded

The Hilbert Transform
In particular, consider two entire functions f and g such that f and g are real-valued on the real line and f ig + is bounded on +  . Then f g =  .
Writing the Cauchy kernel as a combination of the Poisson and Hilbert kernels, we get: where PV denotes the Cauchy principal value,  the Hilbert transform, The above formula is valid for all ( ) not defined when 0 t = ). This agrees with the fact that, although the expression This is the Cauchy integral of  , which is an analytic function bounded on the real line-so, for fixed x and y ∈  , ( ) Proposition 3.5. The exact solution of Equation (1.2) is given by: Proof. By analogy with (3.18) and following Constandinides-Marhefka [8], let the incomplete Airy functions be given by: Ai e d , For the integral to converge, the upper integration limit in (3.26) can be taken to be e iθ ∞ , for 2 4 5 0, , , Over each range, the integral is constant.
Due to the integrand's rapid decay, the integral (3.26) converges for any , x k ∈  , defining an analytic function of both.
Since (by a change of variable t it ′ = ) Ci e d e d , x ∈  for k ∈  and x ∈  . Journal of Applied Mathematics and Physics By the Schwartz reflection principle, it follows that , is also bounded. Therefore: Now Let's apply the same change of variables to the derivative of the above function with respect to x. We have: , which are both bounded. Now, we have the general formula of the Hilbert transform for the function Ci .
k k Using Equation (3.28), the above system becomes: Therefore, the solution for α , β , and γ : The first line can also be used to obtain other bounds, e.g. x t λ = , but covers both the case when x is fixed and t is large and the case when t is fixed and x is large.

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 simpli-Journal of Applied Mathematics and Physics

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 Table A1.
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.

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 Figure B1.