An integral equation–based numerical method for the forced heat equation on complex domains

Integral equation–based numerical methods are directly applicable to homogeneous elliptic PDEs and offer the ability to solve these with high accuracy and speed on complex domains. In this paper, such a method is extended to the heat equation with inhomogeneous source terms. First, the heat equation is discretised in time, then in each time step we solve a sequence of so-called modified Helmholtz equations with a parameter depending on the time step size. The modified Helmholtz equation is then split into two: a homogeneous part solved with a boundary integral method and a particular part, where the solution is obtained by evaluating a volume potential over the inhomogeneous source term over a simple domain. In this work, we introduce two components which are critical for the success of this approach: a method to efficiently compute a high-regularity extension of a function outside the domain where it is defined, and a special quadrature method to accurately evaluate singular and nearly singular integrals in the integral formulation of the modified Helmholtz equation for all time step sizes.

We are grateful for the support from the Natural Science and Engineering Research Council of Canada.

Open access funding provided by Royal Institute of Technology. We recevied support of the Swedish Research Council under Grant No. 2015-04998 and funding from the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine.

Authors and Affiliations

1. Department of Mathematics, KTH Royal Institute of Technology, Stockholm, Sweden

   Fredrik Fryklund & Anna-Karin Tornberg

2. Department of Mathematics, Simon Fraser University, Burnaby, Canada

   Mary Catherine A. Kropinski

Corresponding author

Correspondence to Fredrik Fryklund.

Communicated by: Leslie Greengard

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix A. Adaptive time stepping with IMEX Runge-Kutta methods

1.1 A.1 Adaptive discretisation in time

This appendix shows how applying implicit-explicit Runge-Kutta (IMEXRK) schemes from [5] to the heat equation reduces it to a sequence of modified Helmhotlz equations to solve at each time step. Formulate the heat equation (1)–(3) as:

where the superscripts denote implicit and explicit, referring to the term being classified as stiff or nonstiff, respectively.

Let t_N denote an instance in time that is the sum of previous discrete time steps \(\{\delta t_{i}\}_{i = 1}^{N-1}\) that may be of different size:

for some initial time t₀. Let U_N be the approximation of U(t_N), then at time t_N+ 1 it is:

where N_S is the number of stages for k^σ, σ ∈ {I, E}, computed as

The second argument of F^σ in (A.5) is defined as:

and \(\bar {U}_{1} = U_{N}\). The coefficients \(\{a_{i,j}^{\sigma }\}_{i,j = 1}^{N_{S}}\), \(\{b_{j}^{\sigma }\}_{j = 1}^{N_{S}}\) and \(\{c_{j}^{\sigma }\}_{j= 1}^{N_{S}}\) are tabulated in the two associated Butcher tableaus for σ = I and σ = E; see Table 3 for a general IMEXRK scheme. The principal difference between the coefficients for implicit and explicit methods is that \(a_{i,j}^{E} = 0\) for i ≤ j while \(a_{i,j}^{I} \neq 0\) for i = j, excluding i = 1. The quantity \(\bar {U}_{i}\) is unknown for every i = 2, … , N_S, since the corresponding implicit stage \({k^{I}_{i}}\) is unknown.

The implicit stage at i is \({k^{I}_{i}} ={\Delta } \bar {U}_{i}\) by definition (A.2). To avoid approximating the differential operator, replace \({k^{I}_{i}}\) in (A.6) with \({\Delta } \bar {U}_{i}\) and reformulate as:

The (A.7) has the form of the modified Helmholtz equation (4)–(5): f(x) corresponds to the right-hand side, \(u(\mathbf {x}) = \bar {U}_{i}(\mathbf {x})\) and \(\alpha ^{2} = (\delta t_{N}a_{i,i}^{I})^{-1}\). We stress that the larger α² is the harder (4)–(5) is to solve accurately in terms of numerics (see Section 3.2.1). The associated boundary condition g is (3) evaluated at time \(t_{N}+\delta t_{N}{c_{i}^{I}}\).

To obtain the next stage \({k^{I}_{i}}\), equation (A.7) must be solved for \(\bar {U}_{i}\) in Ω. Once \(\bar {U}_{i}\) is known, reformulate (A.7) and compute:

With \({k^{I}_{i}}\) known the stage \({k^{E}_{i}}\), that is F^E, can be computed explicitly. Note that for (1)–(3) F^E = F(t, x), so the explicit stage \({k^{E}_{i}}\) is independent of the implicit stages; thus, it is computed directly.

To summarise, the approximate solution U^N+ 1 at time t^N+ 1 is given by (A.4). The stages \({k^{I}_{i}}\), for i = 1, … , N_S are obtained by solving (4)–(5), corresponding to (A.7), and explicit computation of (A.8). Once \(\bar {U}_{i}\) is known, \({k_{i}^{E}} = F^{E}\left (t_{N}+\delta t_{N}{c_{i}^{E}},\bar {U}_{i}\right )\) is computed explicitly. See the flowchart in Appendix B for a graphical overview.

1.1.1 A.1.1 IMEXRK2

This scheme is never used in this paper, but serves as a simple example of applying an IMEX Runge-Kutta scheme. The stencil for IMEXRK2, with coefficients tabulated in Table 5, involves taking a half time step δt_N/2 and solving for \(\bar {U}_{2}\) satisfying:

By (A.4) the solution at the next time step, t_N+ 1 = δt_N + t_N for every x ∈Ω is:

An important aspect of IMEXRK2 is that we obtain a second-order method by only solving (4)–(5) once, i.e. only one intermediate stage is required.

An adaptive time stepper can be constructed by coupling IMEXRK2 with a method of lower order. A simple IMEX scheme of first order is the Forward-Backward Euler scheme, with coefficients given in Table 4. Applied to the heat equation (1), we have:

1.1.2 A.1.2 IMEXRK34

The IMEKRK34 scheme is a coupled third- and fourth-order scheme; see Tables 6 and 7 for the associated Butcher tableaus. It has two sets of six stages \(\{k_{i}^{\sigma }\}_{i = }\) for σ = I, E, but only five implicit stages need to be solved for every iterate in time [5]. This is due to \({k^{I}_{6}}\) at t_N is equal to \({k^{I}_{1}}\) at t_N+ 1 for N > 1, a property sometimes referred to as first same as last, or FSAL. Note that the explicit stages do not have this property.

For N = 0, the first stage must be given by supplementary initial data ΔU₀. Otherwise, the procedure is exactly as described in Appendix A.1: for a given i solve (A.7) for \(\bar {U}_{i}\). Once known extract \({k^{I}_{i}} = {\Delta } \bar {U}_{i}\) from (A.7) and compute \({k_{i}^{E}} = F^{E}\left (t_{N}+\delta t_{N}{c_{i}^{E}},\bar {U}_{i}\right )\) explicitly and start over for i + 1 until all six stages are known. An approximate solution U_N+ 1 at t_N+ 1 is given by (A.4), which is a fourth-order approximation. The third-order approximation \(\tilde {U}_{N+1}\) is given by (A.4) as well, but with the coefficients \(\{\tilde {b}_{j}^{\sigma }\}_{j = 1}^{N_{S}}\) instead of \(\{b_{j}^{\sigma }\}_{j = 1}^{N_{S}}\) (Table 6).

1.1.3 A.1.3 Adaptivity

Denote the solution given by Forward-Backward Euler or the third-order method in IMEXRK34 as \(\tilde {U}(\mathbf {x})\). At each discrete time instance t_N+ 1 = δt_N + t_N for some δt_N, we compute U_N+ 1(x) and \(\tilde {U}_{N+1}(\mathbf {x})\). The relative temporal error is approximated by:

where ∥⋅∥ is the standard discrete ℓ₂-norm (50). If r is less than some tolerance TOL, then U_N+ 1(x) is accepted as solution at time t_N+ 1. The quantity r is used to produce a new time step δt_{N, NEW} by:

where p = 2 from the order of the IMEXRK2 scheme and p = 4 for IMEXRK34. The value 0.9 is a safety factor. If the solution is accepted, then δt_N+ 1 = δt_{N, NEW}; otherwise, the computations start over at t_N with δt_N = δt_{N, NEW}. Thus, even if the solution is accepted, the step size is updated by the scheme (A.15), meaning growth is possible if appropriate. See the flowchart in Appendix B for a graphical overview.

Appendix B. Flowchart over solution procedure

Flowchart over the procedure for updating the approximate solution U_N at t_N for the heat equation (1)–(3). Note that the grey block corresponds to the flowchart in Fig. 16

Full size image

Flowchart over the procedure for solving the modified Helmholtz equation (4)–(5). The numbers associated with the yellow boxes denote the major computational steps; see Section 3.3

Full size image

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions