On the numerical solution of the heat equation I: Fast solvers in free space

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

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(NMlogN), 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.

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 artificial (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], [15], [16], [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 Rd:Ut(x,t)-2U(x,t)=f(x,t),t>0,in the absence of physical boundaries, subject to the initial conditionU(x,0)=U0(x),xRd.The functions f(x,t) and U0(x) are assumed to be compactly supported in the box B=[-R/2,R/2]d, centered at the origin. We also assume that f(x,t) and U0(x) are k-times differentiable: f(x,t)Ck(B×[0,T]) and U0(x)Ck(B). In subsequent works we will consider complex geometry and discontinuous data.

From standard potential theory [13], [23], the solution can be written asU(x,t)=Rdkx-y,tU0(y)dy+0tRdk(x-y,t-τ)f(y,τ)dydτ,where the fundamental solution for the heat equation in Rd isk(x,t):=e-x2/4t(4πt)d2,t>0.

We will refer to the first integral in (3) as an initial potential and the second integral as a volume potential. There is a substantial literature on Green’s function methods for problems of diffusion (see for example [3]). However, straightforward discretization of the above integrals leads to an enormously expensive numerical scheme – the solution is dependent on the full space–time history of the diffusion process. With N points in the discretization of the domain and M time steps, it is easy to see that O(N2M+N2M2) work is required. Thus, the obvious advantages of the approach (stability, robustness, and the correctness of the far-field behavior) appear to be overwhelmed by the problems of cost.

In recent years, however, several fast algorithms have been developed [10], [11], [12], [26] that lead to nearly optimal schemes, for which the work required is O(MNlogN). Related methods can be found in [18], [20], [27]. They involve a fairly substantial amount of numerical and analytic machinery. In the present context, where we need to evaluate volume potentials with smooth data, a simpler method can be developed based entirely on the continuous Fourier transform. After outlining the algorithm itself, we illustrate its performance with some numerical examples from materials science, simulating dendritic solidification.

Section snippets

Fourier representation of the solution

While (3) describes the solution to the heat equations (1), (2), significant advantage can be obtained by considering its Fourier transform. For this, we letUˆ(s,t)=Rdeis·xU(x,t)dx,U(x,t)=1(2π)dRde-is·xUˆ(s,t)ds.It is obvious from (1), (2), and well known that Uˆ(s,t) satisfies the ordinary differential equationdUˆdt(s,t)=-s2Uˆ(s,t)+fˆ(s,t),wherefˆ(s,t)=Rdeis·xf(x,t)dx.An elementary calculation shows thatUˆ(s,t)=e-s2ΔtUˆ(s,t-Δt)+Φ(s,t,Δt),whereΦ(s,t,Δt)=t-Δtte-s2(t-τ)fˆ(s,τ)dτ.Thus,

The numerical scheme

With the full complement of tools in place, we can now provide an informal description of the fast recursive marching (FRM) method. For kth order accuracy in time, we refer to the method as FRM(k).
The FRM(2) method

  • Step 1:

    Initialization

    • (a) Select time step Δt and number of time steps M.

    • (b) Select precision ϵ for quadrature in Fourier domain.

    • (c) Obtain N1 quadrature nodes and weights in Fourier space for (18) according to Corollary 1.

    • (d) Compute weights W0, W1 from (10).

  • Step 2:

    Transform initial data

    • (a)

Numerical results

We first test the performance of the algorithm in two dimensions with a known exact solution Uexa(x1,x2,t), constructed so that Uexa(x1,x2,t)=V(x1,x2,t)+W(x1,x2,t), where V corresponds to an initial potential and W corresponds to a volume potential. For this, we letV(x1,x2,t)=e-(x1-0.05)2+(x2-0.15)24(t+0.01)4π(t+0.01).Clearly V satisfies the homogeneous heat equation with initial dataV(x1,x2,0)=e-(x1-0.05)2+(x2-0.15)24(0.01)4π(0.01).We let W satisfy the inhomogeneous heat equationWt=2W+f,W(x1,x

Conclusions and generalizations

We have described the fast recursive marching method, a simple Fourier-based method for the solution of the heat equation in free space with smooth initial data and a smooth source term. It allows for efficient and accurate long-time simulations without the need for artificial boundary conditions on a finite computational domain. The convergence theory can be stated trivially – the error in the solution is the quadrature error in computing the space–time integral (3). The CPU time of the method

Acknowledgment

This work was supported by the Applied Mathematical Sciences Program of the U.S. Department of Energy under Contract DEFG0288ER25053.

References (27)

  • Gregory Beylkin

    On the fast Fourier transform of functions with singularities

    Appl. Comput. Harmonic Anal.

    (1995)
  • C.A. Brebbia et al.

    Boundary Element Techniques

    (1984)
  • P.J. Davis et al.

    Methods of Numerical Integration

    (1984)
  • Cited by (40)

    View all citing articles on Scopus
    View full text