A fast method for solving the heat equation by layer potentials

https://doi.org/10.1016/j.jcp.2006.11.001Get rights and content

Abstract

Boundary integral formulations of the heat equation involve time convolutions in addition to surface potentials. If M is the number of time steps and N is the number of degrees of freedom of the spatial discretization then the direct computation of a heat potential involves order N2M2 operations. This article describes a fast method to compute three-dimensional heat potentials which is based on Chebyshev interpolation of the heat kernel in both space and time. The computational complexity is order p4q2NM operations, where p and q are the orders of the polynomial approximation in space and time.

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 form12u(x,t)=Ku(x,t)-Vun(x,t),xS,t>0.In three spatial variables, the single- and double layer potential are given byVgs(x,t)=0tSG(x-y,t-τ)gs(y,τ)dSydτ,Kgd(x,t)=0tSnyG(x-y,t-τ)gd(y,τ)dSydτ,where S is the boundary surface, and G is the heat kernel which is given byG

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 definitionTn(x)=cncos(narccos(x)),cn=12,n=0,1,n>0.Note that the inclusion of the factor cn is not standard, but will simplify the formulas below. The roots of Tn+1 are given byωkn=cosπ22k+1n+1,0knand the Lagrange polynomials with respect to these nodes areLi(x)=k=0kinx-ωknωin-ωkn.The interpolation Tnf 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 sourceu(x,t)=1(4πt)3/2exp-x-x024tat 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 O(Δ3/2) 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)

  • Kendall E. Atkinson

    The Numerical Solution of Integral Equations of the Second Kind

    (1997)
  • C. Canuto et al.

    Spectral Methods in Fluid Dynamics

    (1988)
  • D. Chien

    Numerical evaluation of surface integrals in three dimensions

    Math. Comp.

    (1995)
  • G.F. Dargush et al.

    Application of the boundary element method to transient heat conduction

    Int. J. Numer. Meth. Eng.

    (1991)
  • L. Greengard et al.

    Spectral approximation of the free-space heat kernel

    Appl. Comput. Harmonic Anal.

    (1999)
  • L. Greengard et al.

    A fast algorithm for the evaluation of heat potentials

    Comm. Pure Appl. Math.

    (1990)
There are more references available in the full text version of this article.

Cited by (0)

View full text