On the numerical solution of the heat equation I: Fast solvers in free space Jing-Rebecca Li * , Leslie Greengard INRIA-Rocquencourt, Projet POEMS, Domaine de Voluceau – Rocquencour, 78153 Le Chesnay Cedex, France Received 11 April 2006; received in revised form 19 June 2007; accepted 19 June 2007 Available online 7 July 2007 Abstract We describe a fast solver for the inhomogeneous heat equation in free space, following the time evolution of the solution in the Fourier domain. It relies on a recently developed spectral approximation of the free-space heat kernel coupled with the non-uniform fast Fourier transform. Unlike finite difference and finite element techniques, there is no need for artificial boundary conditions on a finite computational domain. The method is explicit, unconditionally stable, and requires an amount of work of the order OðNM log N Þ, where N is the number of discretization points in physical space and M is the number of time steps. We refer to the approach as the fast recursive marching (FRM) method. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Heat equation; Free space; Unbounded domain; Integral representation; Spectral approximation 1. Introduction The solution of the heat equation (the diffusion equation) in free space or in unbounded regions arises as a modeling task in a variety of engineering, scientific, and financial applications. While the most commonly used approaches are based on finite difference (FD) and finite element (FE) methods, these must be coupled to arti- ficial (non-reflecting) boundary conditions imposed on a finite computational domain in order to simulate the effect of diffusion into an infinite medium. These boundary conditions are discussed, for example, in [5,14–17]. Here, we describe a mathematically much more straightforward approach, which we will refer to as the fast recursive marching (FRM) method. It is based on evaluating the exact solution of the governing equation, using convolution in space and time with the free-space Green’s function. One advantage of this approach is that essentially no convergence theory is required. The error in the solution is simply the quadrature error in evaluating the solution. In the present paper, we restrict our attention to the simplest setting, namely the isotropic inhomogeneous heat equation in R d : U t ðx; tÞr 2 U ðx; tÞ¼ f ðx; tÞ; t > 0; ð1Þ 0021-9991/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2007.06.021 * Corresponding author. Tel.: +33 1 39 63 53 55; fax: +33 1 39 63 58 84. E-mail address: jingrebecca.li@inria.fr (J.-R. Li). Journal of Computational Physics 226 (2007) 1891–1901 www.elsevier.com/locate/jcp