Fast and high-order approximation of parabolic equations using hierarchical direct solvers and implicit Runge-Kutta methods

TL;DR

This work presents a fast, high‑order solver for parabolic PDEs by coupling the hierarchical Poincaré–Steklov direct elliptic solver (HPS) with high‑order Runge–Kutta time stepping. By treating the stiff linear term with ESDIRK and the nonlinear/nonstiff terms with ERK in an IMEX framework, and leveraging a fixed implicit matrix across stages, the method achieves offline cost $\mathcal{O}(N^{1.5})$ and per‑step cost $\mathcal{O}(N\log N)$ in 2D. The authors prove stability for the backward‑Euler case and provide extensive numerical evidence suggesting stability up to order five for a broad class of problems, including 1D and 2D tests such as diffusion, convection–diffusion, Burgers, and nonlinear Schrödinger dynamics. The approach yields high‑order accuracy and scalable performance on large grids, with detailed comparisons to FEM highlighting computational advantages of the hierarchical direct solver. Together, these results indicate a practical, robust framework for efficient high‑order time integration of parabolic equations with strong spatial discretization via HPS.

Abstract

An additive Runge-Kutta method is used for the time stepping, which integrates the linear stiff terms by an explicit singly diagonally implicit Runge-Kutta (ESDIRK) method and the nonlinear terms by an explicit Runge-Kutta (ERK) method. In each time step, the implicit solve is performed by the recently developed Hierarchical Poincaré-Steklov (HPS) method. This is a fast direct solver for elliptic equations that decomposes the space domain into a hierarchical tree of subdomains and builds spectral collocation solvers locally on the subdomains. These ideas are naturally combined in the presented method since the singly diagonal coefficient in ESDIRK and a fixed time-step ensures that the coefficient matrix in the implicit solve of HPS remains the same for all time stages. This means that the precomputed inverse can be efficiently reused, leading to a scheme with complexity (in two dimensions) $\mathcal{O}(N^{1.5})$ for the precomputation where the solution operator to the elliptic problems is built, and then $\mathcal{O}(N \log N)$ for the solve in each time step. The stability of the method is proved for first order in time and any order in space, and numerical evidence substantiates a claim of stability for a much broader class of time discretization methods. Numerical experiments supporting the accuracy of efficiency of the method in one and two dimensions are presented.

Figures (14)

• Figure 1: The domain $\Omega$ is hierarchically halved and eventually is partitioned into squares. Then each square is discretized with $p\times p$ Chebyshev nodes.
• Figure 2: An illustration of boundary points on the parent domain $\Omega^\tau = \Omega^\alpha \cup \Omega^\beta$.
• Figure 3: The vector $J^\tau$ is partitioned into five blocks.
• Figure 4: Eigenvalues of the time-stepping map for 1D variable coefficient parabolic equation
• Figure 5: time stamps of convection diffusion equation

Theorems & Definitions (3)

• Theorem 1
• proof
• Remark 1