Efficient highorder integral equation methods have been developed for solving the boundary value problems of the heat equation with complex geometries in two and three dimensions. First of all, the classical heat potential theory is applied to convert such problems to Volterra integral equations of the second kind via the heat layer potentials. Some advantages of the integral formulation as compared with standard finite difference and finite element methods include reduction of the dimension of the problem by one, high order accuracy, unconditional stability, insensitivity to different geometries, and elimination of truncating the computational domain and the need of artificial boundary conditions for exterior problems. However, the heat layer potentials contains convolution integrals in both space and time whose direct evaluation requires O(NS2NT2) work and O(NSNT) storage, where NS is the total number of discretization points in the spatial boundary and NT is the total number of time steps. This is excessively expensive even for problems of modest size, especially for threedimensional problems.
In order to evaluate the heat layer potentials accurately and efficiently, they are split into two parts  the local part which contains the temporal integration from t  δ to t and the history part which contains the temporal integration from 0 to t  δ. For the local part, Product integration is applied on the temporal integral to convert it to a sum of several spatial convolution integrals where the socalled local kernels have logarithmic singularity in two dimensions and 1r singularity in three dimensions. These weakly singular integrals are discretized via highorder quadratures and the resulting discrete summations can then be evaluated via fast algorithms such as the fast multipole method and its descendants.
For the history part, efficient separated sumofexponentials approximations can be constructed for the heat kernel in any dimension. Specifically, in one space dimension, the heat kernel admits an approximation involving a number of terms that is of the order O(log(T δ )(log(1 ) + log log(T δ ))) for any x Î R and δ ≤ t ≤ T, where E is the desired precision. In all higher dimensions, the corresponding heat kernel admits an approximation involving only O (log^{2}(Tδ )) terms for fixed accuracy E. These approximations can be used to accelerate the evaluation of the history part of the heat layer potentials for stationary geometries.
For twodimensional problems with complex stationary geometries, the sumofexponentials approximation is used for the heat kernel and all local and history kernels are compressed only once. The resulting algorithm is very efficient with quasilinear complexity in both space and time for both interior and exterior problems. For twodimensional problems with complex moving geometries, the spectral Fourier approximation is applied for the heat kernel and NUFFT is used to speed up the evaluation of the history part of the heat potentials. The complexity of the algorithm is again quasilinear in both space and time, albeit only for the interior problem. For threedimensional problems, the sumofexponentials approximation is applied to speed up the evaluation of the history part. The singular surface integrals in the local kernels are treated with a spectrally accurate integrator. The algorithm is applicable for both interior and exterior problems and has quasilinear complexity with respect to the temporal variable. All these algorithms can be parallelized in a straightforward manner and their performance is demonstrated with extensive numerical experiments.
