A fast method for solving the heat equation by layer potentials
Introduction
Boundary integral formulations of parabolic problems involve time convolutions in addition to the usual boundary integral operators. While the finite element or finite difference method are highly popular for this type of problems, the method of layer potentials has distinct advantages, in particular, stability of explicit time stepping and reduction of unknowns to the boundary surface. Boundary element methods that evaluate the time convolution have been discussed in the past decades in several places, such as [4], [9], [11].
If the time convolution is computed directly the cost quickly becomes prohibitively expensive with increasing problem size which makes the integral equation approach less attractive than competing approaches. However, the situation changes if fast methods for the computation of heat potentials are considered. One such an algorithm is the fast Gauss transform [7]. In this approach the heat kernel is approximated by a Hermite expansion, which leads to an efficient scheme to compute the spatial convolution efficiently. This approach has been combined with re-starts to avoid time convolutions [12]. The original version of the fast Gauss Transform is based on Hermite expansions of the heat kernel. Since the spatial variables of the kernel separate, the translation operator appear in tensor product form, which can be exploited to reduce the computational cost associated with translation operators [13]. Besides Hermite expansions, exponential expansions have been considered as well, in this case translation operators are diagonal [8].
Another technique that exploits the convolution form of the time integral is to evaluate heat potentials in Fourier domain. This was described for bounded domains in [6], where the Green’s function is expanded in a Fourier series. The method must be combined with nonuniform FFTs since the heat sources are located on a boundary surface. In unbounded domains the kernel appears as a continuous Fourier transform which results in some nontrivial complications, cf. [5].
The approach described in this article is more close to the fast Gauss Transform in that the computation is performed in ‘physical’ and not in Fourier domain. Hence, the method is somewhat simpler when the solution must be evaluated in physical space in every time step, which is the case, for instance, when the heat equation is solved via the Green’s representation theorem. The main features of the discussed approach are as follows:
- •
The heat kernel is expanded in both space and time. The time convolution is computed using a hierarchical scheme, similar to the FMM for one-dimensional source distributions. However, certain modifications are necessary to account for the Volterra form of heat potentials. The cost per time step in this algorithm is almost independent of the length of the time interval.
- •
The heat kernel is interpolated at the Chebyshev nodes in space and time. The translation operators of the multipole method can be easily derived and preserve the tensor product form of the kernel.
After discussing heat potentials and their discretization in Section 2 we introduce the Chebyshev interpolation and the resulting translation operators in Sections 3 Chebyshev expansion of the heat kernel, 4 Translation operators. This leads to the fast algorithm described in Section 5. Section 6 concludes with numerical results.
Section snippets
Heat potentials
It is well known that the solution of the heat equation satisfies Green’s Representation formula. If ut = Δu in the exterior of a domain with boundary S and if u has vanishing initial conditions, then Green’s formula assumes the formIn three spatial variables, the single- and double layer potential are given bywhere S is the boundary surface, and G is the heat kernel which is given by
Chebyshev expansion of the heat kernel
Chebyshev polynomials are commonly used to approximate functions in the interval [−1, 1]. In this article, we use the definitionNote that the inclusion of the factor cn is not standard, but will simplify the formulas below. The roots of Tn+1 are given byand the Lagrange polynomials with respect to these nodes areThe interpolation of a function f at the Chebyshev nodes can be expressed in terms of
Translation operators
The FMM combines sources and evaluation points in a tree of cubes. To accelerate the computation of cube interactions the kernel is replaced by a truncated series expansion. The efficient evaluation of such a series involves the moments of the source cube and the expansion coefficients of the destination cube. In the course of the computation, the moments and coefficients must be translated between the different levels of the cube hierarchy.
In addition to the space variable, heat potentials
A fast algorithm for the time convolution
This section describes in detail the fast method for the smooth potential (10).
Numerical example
To illustrate the convergence behavior we solve the heat equation in the exterior of the unit sphere. The solution approach is based on the Green’s Representation formula (1), desingularized Nyström discretization, and the fast method to evaluate the heat potentials. To obtain a problem with known solution we place the heat sourceat the point x0 = [0.5, 0, 0], supply the Neumann data on the sphere and compare the numerical solution on the sphere with the exact solution
Conclusion
We have discussed a fast method for solving heat conduction problems with layer potentials. The numerical results are in good agreement with the O(p4q2NM) complexity and the quadrature error. It is possible to construct higher-order quadrature schemes that evaluate more terms of the expansion of the local part. The difficulty is that the higher terms involve derivatives of the surface curvature and get very complicated.
For the simplicity of exposition, we have restricted the discussion
References (13)
The Numerical Solution of Integral Equations of the Second Kind
(1997)- et al.
Spectral Methods in Fluid Dynamics
(1988) Numerical evaluation of surface integrals in three dimensions
Math. Comp.
(1995)- et al.
Application of the boundary element method to transient heat conduction
Int. J. Numer. Meth. Eng.
(1991) - et al.
Spectral approximation of the free-space heat kernel
Appl. Comput. Harmonic Anal.
(1999) - et al.
A fast algorithm for the evaluation of heat potentials
Comm. Pure Appl. Math.
(1990)