A Fully Discrete Fast Fourier–Galerkin Method Solving a Boundary Integral Equation for the Biharmonic Equation

Abstract

We develop a fully discrete fast Fourier–Galerkin method for solving a boundary integral equation for the biharmonic equation by introducing a quadrature scheme for computing the integrals of non-smooth functions that appear in the Fourier–Galerkin method. A key step in developing the fully discrete fast Fourier–Galerkin method is the design of a fast quadrature scheme for computing the Fourier coefficients of the non-smooth kernel function involved in the boundary integral equation. We prove that with the proposed quadrature algorithm, the total number of additions and multiplications used in generating the compressed coefficient matrix for the proposed method is quasi-linear (linear with a logarithmic factor), and the resulting numerical solution of the equation preserves the optimal convergence order. Numerical examples are presented to demonstrate the approximation accuracy and computational efficiency of the proposed method.

Appendices

Appendix

Algorithm FEcd

We present here a fast algorithm for evaluating Fourier integrals

$$\begin{aligned} \int _{I_{\mathfrak {ab}}\otimes I^{d-1}} h(\mathbf {x}) e_{-\mathbf {l}}(\mathbf {x})d\mathbf {x}, \end{aligned}$$

(78)

where $\mathbf {l}\in \mathbb {L}_{N,2}$, $I_{{\mathfrak {a}}{\mathfrak {b}}}:=[{\mathfrak {a}},{\mathfrak {b}}]$, ${\mathfrak {a}},{\mathfrak {b}}\in I$ and $d\in \{2,3\}$ and h is a continuous function on $I_{\mathfrak {ab}}\otimes I^{d-1}$. The algorithm is a modification of the fast Fourier transform on sparse grids, proposed in [42], based on the multiscale Lagrange interpolation of h on sparse grids in [40].

We recall the multiscale Lagrange interpolation on sparse grids in $I^d$ [40]. The multiscale Lagrange interpolation on sparse grids is builded upon the multiscale Lagrange interpolation on I introduced in [20], see also [21]. Following [20], we choose two contractive mappings $\phi _\rho : I\rightarrow I,~\rho \in \mathbb {Z}_2$, defined by $\phi _0(x) := \frac{x}{2}$ and $\phi _1(x) := \frac{x}{2}+\pi ,~x \in I$, and let $\psi :=\{\phi _\rho :\rho \in \mathbb {Z}_2\}$. Fix $m\in \mathbb {N}$ and let $V_0:=\{0\leqslant \zeta _0<\zeta _1<\cdots<\zeta _{m-1}<2\pi \}$ be a refinable set in the sense that $V_0\subset \phi _0(V_0)\cup \phi _1(V_0)$. Let $v_{0,r}:=\zeta _r$, $r\in \mathbb {Z}_m$. For each $r\in \mathbb {Z}_m$, we define

$$\begin{aligned} l_{0,r}(x):=\prod _{q=0,q\ne r}^{m-1} \frac{x-v_{0,q}}{v_{0,r}-v_{0,q}}, \quad x\in I. \end{aligned}$$

For all $\tau \in \mathbb {N}$ and ${\mathbf {p}}_\tau :=[\rho _\gamma :\gamma \in \mathbb {Z}_{\tau }]\in \mathbb {Z}_2^{\tau }$, we let $\phi _{{\mathbf {p}}_\tau }:=\phi _{\rho _{\tau -1}}\circ \cdots \circ \phi _{\rho _{0}}$ and $\mu ({\mathbf {p}}_\tau ):=\sum _{\gamma \in \mathbb {Z}_\tau } \rho _\gamma 2^\gamma$. The points $v_{\tau ,r}:=\phi _{{\mathbf {p}}_\tau }(\zeta _{r'})$, with ${\mathbf {p}}_\tau \in \mathbb {Z}_2^\tau ,~r'\in \mathbb {Z}_m$ and $r=m\mu ({\mathbf {p}}_\tau )+r'$, will be used as the interpolation points at higher levels. Specifically, for all $\tau \in \mathbb {N}$, we define

$$\begin{aligned} l_{\tau ,r}(x):=\left\{ \begin{array}{ll} \prod \nolimits _{q=m\mu ({\mathbf {p}}_\tau ),q\ne r}^{m\mu ({\mathbf {p}}_\tau )+m-1} \frac{x-v_{\tau ,q}}{v_{\tau ,r}-v_{\tau ,q}}, &{} \quad x\in \phi _{{\mathbf {p}}_\tau }(I), \\ 0, &{} \quad x\in I \setminus \phi _{{\mathbf {p}}_\tau }(I). \end{array} \right. \end{aligned}$$

We now describe the multiscale interpolation functionals. Let $\mathbb {W}_{0,m}:=\mathbb {Z}_m$. For all $\tau \in \mathbb {N}$, denote $V_\tau := \{v_{\tau ,r}, r\in \mathbb {Z}_{2^\tau m}\}$ and $\mathbb {W}_{\tau ,m}:=\{r\in \mathbb {Z}_{2^\tau m}: v_{\tau ,r}\in V_\tau \setminus V_{\tau - 1}\}$. Let $a_{q,r}:=l_{0,q}(v_{1,r})$, for $q\in \mathbb {Z}_m$ and $r\in \mathbb {Z}_{2m}$. For $k,l\in \mathbb {N}$, $k\mod l$ is the remainder of the division of k by l. Define linear functionals $\eta _{\tau ,r}: C(I)\rightarrow \mathbb {R}$, for $\tau \in \mathbb {N}_0$ and $r\in \mathbb {W}_{\tau , m}$, by

$$\begin{aligned} \eta _{\tau ,r}(f) :=\left\{ \begin{array}{ll} f(v_{\tau ,r}) - \sum \nolimits _{q\in \mathbb {Z}_m} a_{q,r\mod 2m}f(v_{\tau -1, m \lfloor \frac{r}{2m}\rfloor +q}),&{} \quad \tau >0,\\ f(v_{0,r}), &{} \quad \tau =0. \end{array} \right. \end{aligned}$$

Next, we describe the multiscale Lagrange interpolation of f on sparse grids. For all ${{\varvec{\tau }}}:=[\tau _k: k\in \mathbb {Z}_d]\in \mathbb {N}_0^d$, we define index set $\mathbb {W}_{{\varvec{\tau }},m}:=\mathbb {W}_{\tau _0,m}\otimes \mathbb {W}_{\tau _1,m}\otimes \cdots \otimes \mathbb {W}_{\tau _{d-1},m}.$ For all ${{\varvec{\tau }}}:=[\tau _k: k\in \mathbb {Z}_d]\in \mathbb {N}_0^d$ and ${\varvec{r}}:=[r_k: k\in \mathbb {Z}_d]\in \mathbb {W}_{{\varvec{\tau }},m}$, we define $\eta _{{\varvec{\tau }},{\varvec{r}}}:= \eta _{\tau _0, r_0}\otimes \cdots \otimes \eta _{\tau _{d-1},r_{d-1}}$ and $l_{{\varvec{\tau }},{\varvec{r}}}:= l_{\tau _0, r_0}\otimes \cdots \otimes l_{\tau _{d-1},r_{d-1}}.$ For $f\in C(I^d)$ and $N\in \mathbb {N}$, the multiscale Lagrange interpolation of f on sparse grids is defined by

$$\begin{aligned} \mathcal {S}_{N,m}^d f=\sum _{{\varvec{\tau }}\in \mathbb {S}_{N,d}}\sum _{{\varvec{r}}\in \mathbb {W}_{{\varvec{\tau }},m}} \eta _{{\varvec{\tau }},{\varvec{r}}}(f) l_{{\varvec{\tau }},{\varvec{r}}}, \end{aligned}$$

(79)

where $\mathbb {S}_{N,d}:= \Big \{[\tau _k:k\in \mathbb {Z}_d]\in \mathbb {Z}_{N+1}^d: \sum _{k\in \mathbb {Z}_d}\tau _k \leqslant N \Big \}$.

Below, we describe a fast algorithm for evaluating the Fourier integral (78) by employing the multiscale Lagrange interpolation on sparse grids (79). Let $\alpha :=\frac{{\mathfrak {b-a}}}{2\pi }$. For $\mathbf {x}:=[x_0,\mathbf {x}_1]\in I^{d}$ with $x_0\in I$ and $\varsigma \in \mathbb {R}$, we define $\varsigma \cdot \mathbf {x}:=[\varsigma x_0, \mathbf {x}_1]$ and $h_{I_{\mathfrak {ab}}}(\mathbf {x}):=h([\alpha x_0+{\mathfrak {a}},\mathbf {x}_1])$. By a change of variables, the integral in (78) can be reformulated as the following integral on $I^d$

$$\begin{aligned} \int _{I_{\mathfrak {ab}}\times I^{d-1}} h(\mathbf {x}) e_{-\mathbf {l}}(\mathbf {x})d\mathbf {x}=\alpha e_{-l_0}({\mathfrak {a}}) \int _{I^{d}}h_{I_{\mathfrak {ab}}}(\mathbf {x}) e_{-\mathbf {l}}(\alpha \cdot \mathbf {x}) d\mathbf {x}. \end{aligned}$$

(80)

By replacing $h_{I_{\mathfrak {ab}}}$ with the multiscale Lagrange interpolation $\mathcal {S}_{N,m}^d h_{I_{\mathfrak {ab}}}$ on sparse grids, we obtain the quadrature formula for the integral (78)

$$\begin{aligned} Q^d_{N,m,{\mathfrak {a,b}}}(h,\mathbf {l}):=\alpha e_{-l_0}({\mathfrak {a}})\int _{I^{d}} {\mathcal {S}}^d_{N,m}h_{I_{\mathfrak {ab}}}(\mathbf {x}) e_{-\mathbf {l}}(\alpha \cdot \mathbf {x}) d\mathbf {x},\quad \mathbf {l}\in \mathbb {L}_{N,2}. \end{aligned}$$

(81)

To simplify this formula, we introduce necessary notation. For all $N\in \mathbb {N}$, $\tau \in \mathbb {Z}_{N+1}$, $r\in \mathbb {Z}_{2^\tau m}$, $\varsigma \in \mathbb {R}$ and $l\in \mathbb {L}_{N,1}$, we define

$$\begin{aligned} A_{\tau ,r, \varsigma }(l):= \int _{I} l_{\tau , r}(x) e_{- l}(\varsigma x) dx. \end{aligned}$$

For all $N\in \mathbb {N}$, $\tau \in \mathbb {Z}_{N+1}$, $r\in \mathbb {W}_{\tau ,m}$ and $l\in \mathbb {L}_{N,1}$, we define

$$\begin{aligned} A_{\tau , r}(l):= \int _{I} l_{\tau , r}(x) e_{-l}(x) dx. \end{aligned}$$

For all ${{\varvec{\tau }}}:=[\tau _k:k\in \mathbb {Z}_d]\in \mathbb {S}_{N,d}$, ${\varvec{r}}:=[r_k:k\in \mathbb {Z}_d]\in \mathbb {W}_{{\varvec{\tau }},m}$, $\varsigma \in \mathbb {R}$ and $\mathbf {l}:=[l_0,l_1]\in \mathbb {L}_{N,2}$, we define

$$\begin{aligned} A_{{\varvec{\tau }},{\varvec{r}},\varsigma }(\mathbf {l}) := \left\{ \begin{array}{ll} A_{\tau _0, r_0,\varsigma }(l_0)A_{\tau _1, r_1}(l_1),&{} \quad d=2,\\ A_{\tau _0, r_0,\varsigma }(l_0)A_{\tau _1, r_1}(l_1)A_{\tau _2, r_2}(0),&{} \quad d=3. \end{array} \right. \end{aligned}$$

Substituting the interpolation (79) into formula (81), the quadrature formula for the integral (78) can be rewritten as

$$\begin{aligned} Q^d_{N,m,{\mathfrak {a,b}}}(h,\mathbf {l})= \alpha e_{-l_0}({\mathfrak {a}})\sum _{{\varvec{\tau }}\in \mathbb {S}_{N,d}}\sum _{{\varvec{r}}\in \mathbb {W}_{{\varvec{\tau }},m}} \eta _{{\varvec{\tau }},{\varvec{r}}}(h_{I_{\mathfrak {ab}}}) A_{{\varvec{\tau }},{\varvec{r}}, \alpha }(\mathbf {l}), \ \ \text{ for }\ \ \mathbf {l}\in \mathbb {L}_{N,2}. \end{aligned}$$

(82)

From Section 5 of [42], we see that (82) can be rewritten as a dimension-reducible sum. This observation leads us to developing a fast algorithm for computing (82) by employing Algorithm 4.5 in [42], designed for fast computing dimension-reducible sum.

In order to propose a fast algorithm for computing sum (82), we next rewrite it as a dimension-reducible sum. For $\upsilon \in \mathbb {Z}_m$, we let $r_{0,\upsilon }=\upsilon$. For each $\tau \in \mathbb {N}$ and $\upsilon \in \mathbb {Z}_{2^{\tau -1} m }$, we also let $r_{\tau , \upsilon }\in \mathbb {W}_{\tau , m}$ be the number such that there exactly $\upsilon$ elements in $\mathbb {W}_{\tau ,m}$ less than $r_{\tau , \upsilon }$. For ${{\varvec{\tau }}}\in \mathbb {S}_{N,2}$ and $\mathbf {r}\in \mathbb {I}_{{{\varvec{\tau }}}, m}$, we define ${\varvec{\upsilon }}:=[ r_k - \lfloor 2^{\tau _k-1} \rfloor m: k\in \mathbb {Z}_{2}]$ ,

$$\begin{aligned} \eta _\mathbf {r}(h_{I_{\mathfrak {ab}}}):=\left\{ \begin{array}{ll} \eta _{{{\varvec{\tau }}},\mathbf {r}_{{{\varvec{\tau }}},{\varvec{\upsilon }}}}(h_{I_{\mathfrak {ab}}}),&{} \quad d=2,\\ \sum \nolimits _{\tau \in \mathbb {Z}_{N+1-|{{\varvec{\tau }}}|}}\sum \nolimits _{\upsilon \in \mathbb {Z}_{2^{\tau -1} m }}\eta _{[\tau ,{{\varvec{\tau }}}],[r_{\tau ,\upsilon },\mathbf {r}_{{{\varvec{\tau }}},{\varvec{\upsilon }}}]}(h_{I_{\mathfrak {ab}}})A_{\tau , r_{\tau ,\upsilon }}(0),&{} \quad d=3, \end{array} \right. \end{aligned}$$

(83)

and for $\varsigma \in \mathbb {R}$,

$$\begin{aligned} A_{\mathbf {r}, \varsigma }:= A_{{{\varvec{\tau }}},\mathbf {r}_{{{\varvec{\tau }}},{\varvec{\upsilon }}}, \varsigma }, \end{aligned}$$

(84)

where ${\mathbf {r}}_{{{\varvec{\tau }}},{\varvec{\upsilon }}}:=[r_{\tau _k,\upsilon _k}:k\in \mathbb {Z}_2]$. For all $\tau \in \mathbb {N}_0$, we define $\mathbb {V}_{\tau }:=\{r\in \mathbb {N}_0: m\lfloor 2^{\tau -1}\rfloor \leqslant r < m2^{\tau } \},$ where m is the order of polynomial $l_{\tau ,r}$. For all ${{\varvec{\tau }}}:=[\tau _0,\tau _1]\in \mathbb {N}_0^2$, we also define $\mathbb {V}_{{{\varvec{\tau }}}}:=\mathbb {V}_{\tau _0}\otimes \mathbb {V}_{\tau _1}$. For $N\in \mathbb {N}_0$, let

$$\begin{aligned} \mathbb {H}_{N,1}:= \bigcup \limits _{\tau \in \mathbb {Z}_{N+1}} \mathbb {V}_{\tau } \ \ \mathrm{and}~\ \ \mathbb {H}_{N,2}:= \bigcup \limits _{{{\varvec{\tau }}}\in \mathbb {S}_{N,2}} \mathbb {V}_{{{\varvec{\tau }}}}. \end{aligned}$$

With the notation introduced above, for all $\mathbf {l}:=[l_0,l_1]\in \mathbb {L}_{N,2}$ Eq. (82) is rewritten as

$$\begin{aligned} Q^d_{N, m,{\mathfrak {a,b}}} (h,\mathbf {l}) =\alpha e_{-l_0}({\mathfrak {a}})\sum _{\mathbf {r}\in \mathbb {H}_{N, 2}} \eta _\mathbf {r}(h_{I_{\mathfrak {ab}}})A_{\mathbf {r},\alpha }(\mathbf {l}). \end{aligned}$$

(85)

Lemmas 5.2 and 5.3 in [42] show that (85) is a dimension-reducible sum.

We now describe a fast algorithm for computing the dimension-reducible sum in (85). It is shown in [42] that the dimension-reducible sum can be split recursively into one-dimension sums. This idea is employed to design Algorithm 4.5 in [42] for the fast Fourier transform on sparse grids. Here, we modify Algorithm 4.5 in [42] to fit the dimension-reducible sum (85). To this end, we consider a generalized dimension-reducible sum as follows. For all $\varsigma \in \mathbb {R}$, $N\in \mathbb {N}$, ${\varvec{\xi }}_N:=[\xi _{\mathbf {r}}\in \mathbb {R}:{\mathbf {r}}\in \mathbb {L}_{N,2}]$ and each $\mathbf {l}\in \mathbb {L}_{N, 2}$, we define

$$\begin{aligned} \varepsilon _{{\varvec{\xi }}_N,\varsigma }(\mathbf {l}):=\sum _{\mathbf {r}\in \mathbb {H}_{N, 2}} \xi _{\mathbf {r}}A_{\mathbf {r},\varsigma }(\mathbf {l}). \end{aligned}$$

(86)

For all $N\in \mathbb {N}$, $\tau \in \mathbb {Z}_{N+1}$ and $r\in \mathbb {V}_{\tau }$, let ${\varvec{\xi }}_{N, \tau , r}:= [\xi _{[r,r']}: r'\in \mathbb {H}_{N-\tau , 1}].$ For each $\tau \in \mathbb {Z}_{N+1}$, $r\in \mathbb {H}_{N-\tau , 1}$, we also let

$$\begin{aligned} \mathbf {d}_{{\varvec{\xi }}_{N}, \tau , r}:= [d_{r'}:\tau '\in \mathbb {Z}_{\tau +1},r'\in \mathbb {V}_{\tau '}], \end{aligned}$$

where $d_{r'}:=\xi _{[r',r]}$ for $r'\in \mathbb {V}_{\tau }$ and $d_{r'}:=0$ for $r'\in \mathbb {V}_{\tau '}$ with $\tau '< \tau$. For each $\varsigma \in \mathbb {R}$, $\tau \in \mathbb {Z}_{N+1}$ and $r\in \mathbb {V}_{N-\tau }$, we define two one dimensional sums by

$$\begin{aligned} \varepsilon _{{\varvec{\xi }}_{N}, \tau , r}: = \sum _{r'\in \mathbb {H}_{N-\tau , 1}} \xi _{[r,r']} A_{r',1}. \end{aligned}$$

(87)

and

$$\begin{aligned} {\widetilde{\varepsilon }}_{{\varvec{\xi }}_{N}, \tau , r,\varsigma }: = \sum _{\tau '\in \mathbb {Z}_{\tau +1}}\sum _{r'\in \mathbb {V}_{\tau '}} d_{r'} A_{r',\varsigma }. \end{aligned}$$

(88)

For each $\varsigma \in \mathbb {R}$, $\tau \in \mathbb {Z}_{N+1}$ and $\mathbf {l}:=[l_0,l_1]\in \mathbb {L}_{N-\tau ,2}$, define

$$\begin{aligned} \nu _{{\varvec{\xi }}_{N},\tau ,\varsigma }(\mathbf {l}):=\sum _{\tau '\in \mathbb {Z}_{\tau +1}}\sum _{r\in \mathbb {V}_{\tau '}} A_{r,\varsigma }(l_0)\varepsilon _{{\varvec{\xi }}_{N}, \tau ', r}(l_1), \end{aligned}$$

(89)

and

$$\begin{aligned} \displaystyle {\widetilde{\nu }}_{{\varvec{\xi }}_{N},\tau ,\varsigma }(\mathbf {l}):=\sum _{\tau '\in \mathbb {Z}_{N+1}\setminus \mathbb {Z}_{\tau +1}}\sum _{r\in \mathbb {H}_{N-\tau , d-1}} {\widetilde{\varepsilon }}_{{\varvec{\xi }}_{N},\tau ', r,\varsigma }(l_0)A_{r,1}(l_1). \end{aligned}$$

With the notation, we have the following equality

$$\begin{aligned} \varepsilon _{{\varvec{\xi }}_{N},\varsigma }=\nu _{{\varvec{\xi }}_N,\tau ,\varsigma }+{\widetilde{\nu }}_{{\varvec{\xi }}_N,\tau ,\varsigma }. \end{aligned}$$

(90)

The formula (90) shows that the values of $\varepsilon _{{\varvec{\xi }}_{N},\varsigma }$ can be obtained by computing $\nu _{{\varvec{\xi }}_N,\tau ,\varsigma }$ and ${\widetilde{\nu }}_{{\varvec{\xi }}_N,\tau ,\varsigma }$. The definition of $\nu _{{\varvec{\xi }}_{N},\tau ,\varsigma }$ shows that we can obtain the values of $\nu _{{\varvec{\xi }}_{N},\tau ,\varsigma }$ by a one-dimensional algorithm when we have the values of $\varepsilon _{{\varvec{\xi }}_{N},\tau ', r}$. Thus, it remains to develop a method for computing ${\widetilde{\nu }}_{{\varvec{\xi }}_N,\tau ,\varsigma }$. For $\varsigma \in \mathbb {R}$, $\tau \in \mathbb {Z}_{N+1}$, $r\in \mathbb {H}_{N-\tau , 1}$ and $l\in \mathbb {I}_{\tau }$, we let

$$\begin{aligned} {\tilde{h}} _{{\varvec{\xi }}_{N}, \tau , r,\varsigma , l} := \left\{ \begin{array}{ll} {\widetilde{\varepsilon }}_{{\varvec{\xi }}_N,\tau +1, r,\varsigma }(l) + {\tilde{h}} _{{\varvec{\xi }}_{N}, \tau +1, r, \varsigma , l} , &{} \tau \leqslant N-1 \hbox { and } r \in \mathbb {H}_{N-\tau -1,1},\\ 0, &{} \tau \geqslant N \hbox { or } r\in \mathbb {H}_{N-\tau ,1}\setminus \mathbb {H}_{N-\tau -1,1}. \end{array} \right. \end{aligned}$$

(91)

From the definitions of $\displaystyle {\widetilde{\nu }}_{{\varvec{\xi }}_{N},\tau ,\varsigma }$ and ${\tilde{h}} _{{\varvec{\xi }}_{N}, \tau , r,\varsigma , l}$, we have the following equality

$$\begin{aligned} {\widetilde{\nu }}_{{\varvec{\xi }}_{N},\tau ,\varsigma }(\mathbf {l})=\sum _{r\in \mathbb {H}_{N-\tau ,1}}{\tilde{h}} _{{\varvec{\xi }}_{N}, \tau , r,\varsigma , l_0} A_{r,1}(l_1). \end{aligned}$$

(92)

We are ready to present the fast algorithm for computing (86). In this algorithm, we call the one dimensional algorithm “FE1d” described in “Appendix B”, for computing one dimensional sum (88) and (89). For given $\varsigma \in \mathbb {R}$ and ${\varvec{\xi }}_N:=[\xi _{\mathbf {r}}\in \mathbb {R}:{\mathbf {r}}\in \mathbb {H}_{N,2}]$, we define $E_{{\varvec{\xi }}_{N},\varsigma }:=[\varepsilon _{{\varvec{\xi }}_N,\varsigma }(\mathbf {l}): \mathbf {l}\in \mathbb {L}_{N,2}]$. For each $\varsigma \in \mathbb {R}$, $\tau \in \mathbb {Z}_{N+1}$ and $l\in \mathbb {L}_{N-\tau , 1}$, we define $V_{{\varvec{\xi }}_{N},\tau , l,\varsigma }:= [\nu _{{\varvec{\xi }}_N,\tau ,\varsigma }([l',l]): l'\in \mathbb {I}_{\tau } ].$ For each $\varsigma \in \mathbb {R}$, $\tau \in \mathbb {Z}_{N+1}$ and $l\in \mathbb {I}_{\tau }$, we define ${\widetilde{V}}_{{\varvec{\xi }}_{N},\tau , l,\varsigma }:= [{\widetilde{\nu }}_{{\varvec{\xi }}_N,\tau ,\varsigma }([l,l']): l'\in \mathbb {L}_{N-\tau , 1} ].$ For each $\varsigma \in \mathbb {R}$, $\tau \in \mathbb {Z}_{N+1}$ and $l\in \mathbb {I}_{\tau }$, we define ${\widetilde{\mathbf {h}}}_{{\varvec{\xi }}_{N},\tau ,\varsigma , l}:= [{\tilde{h}}_{{\varvec{\xi }}_N,\tau , r,\varsigma , l}: r\in \mathbb {H}_{N-\tau , 1}].$ For each $\tau \in \mathbb {Z}_{N+1}$ and $l\in \mathbb {L}_{N-\tau , 1}$, we let $\mathbf {c}_{{\varvec{\xi }}_{N}, \tau , l}:= [\varepsilon _{{\varvec{\xi }}_{N}, \tau , r}(l): \tau '\in \mathbb {Z}_{\tau +1}, r\in \mathbb {V}_{\tau '} ].$ For all $\tau \in \mathbb {Z}_{N+1}$, $r\in \mathbb {V}_{\tau }$, we define $E_{{\varvec{\xi }}_{N}, \tau ,r}:=[\varepsilon _{{\varvec{\xi }}_{N}, \tau , r}(l): l\in \mathbb {L}_{N-\tau ,1}]$.

Algorithm A.1

$\,\mathrm{FEd}\,(N, {\varvec{\xi }}_{N}, \varsigma )$

Input: $N\in \mathbb {N}$, ${\varvec{\xi }}_{N}:=[\xi _{\mathbf {r}}:\mathbf {r}\in \mathbb {H}_{N,2}]$ and $\varsigma \in \mathbb {R}$

Output: $E_{{\varvec{\xi }}_{N},\varsigma }$.

Step 1. For all $\tau \in \mathbb {Z}_{N+1}$, $r\in \mathbb {V}_{\tau }$, compute $E_{{\varvec{\xi }}_{N}, \tau , r}$ by using $\,\mathrm{FE1d}\,(1, N-\tau , {\varvec{\xi }}_{N, \tau , r}, 1)$.

Step 2. For $\tau \in \mathbb {Z}_{N+1}$, $l\in \mathbb {L}_{N-\tau , 1}$, compute $V_{{\varvec{\xi }}_{N},\tau , l,\varsigma }$ by using $\,\mathrm{FE1d}\,(\tau , \mathbf {c}_{{\varvec{\xi }}_{N}, \tau , l}, \varsigma )$.

Step 3. For all $\tau \in \mathbb {Z}_{N+1}$, $r\in \mathbb {H}_{N-\tau , 1}$, compute $[{\widetilde{\varepsilon }}_{{\varvec{\xi }}_N,{\tau }, r,\varsigma }(l): l\in \mathbb {I}_{\tau } ]$ by using $\,\mathrm{FE1d}\,(\tau , \mathbf {d}_{{\varvec{\xi }}_{N}, \tau , r}, \varsigma )$.

Step 4. Compute $[{\tilde{h}} _{{\varvec{\xi }}_{r}, \tau , r, \varsigma , l}: \tau \in \mathbb {Z}_{N+1}, l\in \mathbb {I}_\tau , r\in \mathbb {H}_{N-\tau , d-1}],$ according to (91).

Step 5. For $\tau \in \mathbb {Z}_{N+1}$, $l\in \mathbb {I}_{\tau }$, compute ${\widetilde{V}}_{{\varvec{\xi }}_{N},\tau , l, \varsigma }$ by using $\,\mathrm{FE1d}\,(1, N-\tau , {\widetilde{\mathbf {h}}}_{{\varvec{\xi }}_{N},\tau ,\varsigma , l}, 1)$.

Step 6. Compute $E_{{\varvec{\xi }}_{N},\varsigma }$ according to (90). Return $E_{{\varvec{\xi }}_{N},\varsigma }$.

Let ${\varvec{\eta }}_{N}:=[\eta _{\mathbf {r}}(h_{I_{\mathfrak {ab}}}): \mathbf {r}\in \mathbb {H}_{N, 2}]$. From the definitions of $Q^d_{N, m,{\mathfrak {a,b}}} (h,\mathbf {l})$ and $\varepsilon _{{\varvec{\xi }}_N,\varsigma }(\mathbf {l})$, it is easy to see that for all $\mathbf {l}\in \mathbb {L}_{N,2}$,

$$\begin{aligned} Q^d_{N, m,{\mathfrak {a,b}}} (h,\mathbf {l})=\alpha e_{-l_0}\varepsilon _{{\varvec{\eta }}_N,\alpha }(\mathbf {l}). \end{aligned}$$

(93)

Thus, applying Algorithm A.1 to ${\varvec{\eta }}_{N}$ leads to the algorithm for computing (85). Let ${\mathbf {Q}}^d_{N, m,{\mathfrak {a,b}}}(h):=[Q^d_{N, m,{\mathfrak {a,b}}} (h,\mathbf {l}):\mathbf {l}\in \mathbb {L}_{N,2}]$ and ${\varvec{\eta }}_{N}:=[\eta _{\mathbf {r}}(h_{I_{\mathfrak {ab}}}): \mathbf {r}\in \mathbb {H}_{N, 2}]$.

Algorithm A.2

${\mathrm {FECd}}(d,N, {\mathfrak {a,b}}, h)$

Input: $d\in \{2,3\}$, $N\in \mathbb {N}$, ${\mathfrak {a}}\in I,~{\mathfrak {b}}\in I$ and $h\in C(I_{\mathfrak {ab}}\otimes I^{d-1})$

Output: ${\mathbf {Q}}_{N, m,{\mathfrak {a,b}}}(h)$.

Step 1. Compute the coefficient ${\varvec{\eta }}_{N}$.

Step 2. Compute $E_{{\varvec{\eta }}_N,\alpha }=[\varepsilon _{{\varvec{\eta }}_N,\alpha }(\mathbf {l}): \mathbf {l}\in \mathbb {L}_{N,2}]$ by calling $\,\mathrm{FEd}\,(N, {\varvec{\eta }}_N, \alpha )$.

Step 3. Compute ${\mathbf {Q}}^d_{N, m,{\mathfrak {a,b}}}(h)$ by (93).

We next estimate the number $\mathcal {M}({\mathrm {FECd}}, N)$ of additions and multiplications used in Algorithm A.2.

Lemma A.3

There exists a positive constant c such that for all $d\in \{2,3\}$ and $N\in \mathbb {N}$,

$$\begin{aligned} \mathcal {M}({\mathrm {FECd}},N) \leqslant cN^22^N. \end{aligned}$$

(94)

Proof

According to Theorem 2.13 in [40], the number of additions and multiplications used in Step 1 is ${\mathcal {O}}(N2^N )$ when $d=2$, and ${\mathcal {O}}(N^22^N )$ when $d=3$. Since the number of additions and multiplications used in the Fourier transform “FE1d” is ${\mathcal {O}}(N2^N)$, and the cardinality of ${\varvec{\eta }}_N$ is ${\mathcal {O}}(N2^N)$, by Theorem 4.6 in [42], the number of additions and multiplications used in Step 2 is ${\mathcal {O}}(N^22^N)$. Since the cardinality of $\mathbb {L}_{N,2}$ is ${\mathcal {O}}(N2^N)$, from Eq. (93) we know that the number of additions and multiplications used in Step 3 is ${\mathcal {O}}(N 2^N)$. Thus, we obtain the estimation (94). $\square$

Lemma A.3 is used in the proof of Theorem 3.4.

Algorithm:  FE1d 

In this subsection, we present the fast Fourier transform “FE1d”, used in the Algorithm A.1, in computing the quantity

$$\begin{aligned} \varepsilon _{{\varvec{\xi }}_N,\varsigma }(l):= \sum _{\tau \in \mathbb {Z}_{N+1}}\sum _{r'\in \mathbb {I}_{\tau }} \xi _{r'} A_{r',\varsigma }(l),\quad l\in \mathbb {L}_{N,1}, \end{aligned}$$

(95)

where $m\in \mathbb {N}$, $N\in \mathbb {N}$, ${\varvec{\xi }}_N:=[\xi _{r'}:r'\in \mathbb {I}_{\tau }, \tau \in \mathbb {Z}_{N+1}]$ and $\varsigma >0$. To this end, we first rewrite $\varepsilon _{{\varvec{\xi }}_N,\varsigma }(l)$ as a combination of the values of $e^{-i2\pi \varsigma l u /2^{N}}$, $u\in \mathbb {Z}_{2^N}$. From the definition of $A_{r',\varsigma }$ in (84), we have for each $r'\in \mathbb {I}_{\tau }$ that $A_{r',\varsigma }=A_{\tau , r, \varsigma }$ where $r=r_{\tau ,r'-\lfloor 2^{\tau -1}\rfloor m}$. Thus, we can rewrite (95) as

$$\begin{aligned} \varepsilon _{{\varvec{\xi }}_{N},\varsigma }(l)&= \sum _{\tau \in \mathbb {Z}_{N+1}} \sum _{r\in \mathbb {W}_{\tau ,m}}\xi _{\tau , r} A_{\tau , r,\varsigma }(l), \end{aligned}$$

where $\xi _{\tau , r}:=\xi _{r'}$. For $q\in \mathbb {Z}_m$ and $\kappa \in \mathbb {Z}_{2m}$, we let $a_{q,\kappa }:=\ell _{0,q}(v_{1,\kappa })$. From Lemma 2.1 in [40], the Lagrange basis functions $l_{\tau , r}$,$\tau \in \mathbb {N}_0, r\in \mathbb {W}_{\tau ,m}$ satisfy the following refinement equation

$$\begin{aligned} l_{\tau , r} =\sum _{q\in \mathbb {Z}_{2m}} a_{r \mod m, q} l_{\tau +1, 2m\lfloor r/m\rfloor +q}. \end{aligned}$$

Thus, we have for all $\tau \in \mathbb {N}_0, r\in \mathbb {W}_{\tau ,m}$ that

$$\begin{aligned} A_{\tau , r, \varsigma } =\sum _{q\in \mathbb {Z}_{2m}} a_{r \mod m, q} A_{\tau +1, 2m\lfloor r/m\rfloor +q, \varsigma }. \end{aligned}$$

By a direct computation, we rewrite $\varepsilon _{{\varvec{\xi }}_{N},\varsigma }(l)$ as

$$\begin{aligned} \varepsilon _{{\varvec{\xi }}_{N},\varsigma }(l)= \sum _{r\in \mathbb {Z}_{2^N m}} {\bar{\xi }}_{N, r} A_{N, r,\varsigma }(l),\quad l\in \mathbb {L}_{{N},1}, \end{aligned}$$

(96)

where ${\bar{\xi }}_{N, r}$ is defined recursively by for $\tau \in \mathbb {Z}_{N+1}$,

$$\begin{aligned} {\bar{\xi }}_{\tau , r}:= \left\{ \begin{array}{ll} \xi _{0, r}, &{} \quad \tau =0, \\ \xi _{\tau , r} + \sum \nolimits _{q\in \mathbb {Z}_m} a_{q, r\mod 2m} {\bar{\xi }}_{\tau -1, m\lfloor r/2m\rfloor +q}, &{} \quad \tau>0, r\in \mathbb {W}_{\tau ,m}, \\ \sum \nolimits _{q\in \mathbb {Z}_m} a_{q, r\mod 2m} {\bar{\xi }}_{N-1, m\lfloor r/2m\rfloor +q} , &{} \quad \tau >0, r\in \mathbb {Z}_{2^\tau m}\setminus \mathbb {W}_{\tau , m}. \end{array} \right. \end{aligned}$$

(97)

For all $\varsigma >0$ and $q\in \mathbb {Z}_{m}$ and $l\in \mathbb {Z}$, we let

$$\begin{aligned} {\tilde{t}}_{N,q,\varsigma }(l) := \frac{1}{2^N} \frac{1}{\sqrt{2\pi }}\int _I l_{0,q}(x)e^{-i\varsigma lx/2^{N}}dx. \end{aligned}$$

(98)

By a simple computation, we obtain for all $r\in \mathbb {Z}_{2^Nm}$, $\varsigma >0$ and $l\in \mathbb {Z}$ that

$$\begin{aligned} A_{N, r,\varsigma }(l)={\tilde{t}}_{N,r\mod m,\varsigma }(l)e^{-i2\pi \varsigma l u/2^N}, \end{aligned}$$

(99)

where $u=\lfloor r/m \rfloor$. Let $\left( {\hat{\xi }}_{N,q,\varsigma } \right) _l:=\sum \limits _{u\in \mathbb {Z}_{2^N}} {\bar{\xi }}_{N,mu+q} e^{-i2\pi \varsigma l u /2^{N}}$. Then, from (96) and (99), we have that

$$\begin{aligned} \varepsilon _{{\varvec{\xi }}_{N},\varsigma }(l)=\sum _{q\in \mathbb {Z}_{m}} {\tilde{t}}_{N,q,\varsigma }(l) \left( {\hat{\xi }}_{N,q,\varsigma }\right) _l, \quad l\in \mathbb {L}_{{N},1}. \end{aligned}$$

(100)

Equation (100) and the definition of $\left( {\hat{\xi }}_{N,q,\varsigma } \right) _l$ show that we can construct a fast algorithm for computing $\varepsilon _{{\varvec{\xi }}_{N},\varsigma }(l)$ by the fast Fourier transform in [25] when $\varsigma =1$. In the case $\varsigma >0$ and $\varsigma \ne 1$, $\left( {\hat{\xi }}_{N,q,\varsigma } \right) _l$ is the fractional Fourier transform of $\left[ {\bar{\xi }}_{N,mu+q}:u\in \mathbb {Z}_{2^N}\right]$. In this case, we employ the fast algorithm for fractional Fourier transform in [9] to construct a fast algorithm for computing $\varepsilon _{{\varvec{\xi }}_{N},\varsigma }(l)$.

Below, we recall the fast algorithm for fractional Fourier transform for ready reference. Let $\varpi , {\tilde{\varpi }}\in \mathbb {N}$ and ${\mathbf {x}}:=[x_u\in \mathbb {R}:u\in \mathbb {Z}_{\varpi }]$. We consider computing

$$\begin{aligned} G_\varsigma (\mathbf {x}, l+s):=\sum _{u=0}^{\varpi -1} x_u e^{-i 2\pi (l+s)u \varsigma /\varpi },\qquad l\in \mathbb {N}~\mathrm{with}~ 0\leqslant l<{\tilde{\varpi }}, \end{aligned}$$

where $s\in \mathbb {Z}$ and $\varsigma >0$ is a fractional number. First we define the sequences $\mathbf {y}:=[y_u\in \mathbb {C}:u\in \mathbb {Z}_{\varpi +{\tilde{\varpi }}}]$ and $\mathbf {z}:=[z_u\in \mathbb {C}:u\in \mathbb {Z}_{\varpi +{\tilde{\varpi }}}]$ by

$$\begin{aligned} \begin{array}{lll} y_u:=&{} x_u e^{-i\pi u^2\varsigma /\varpi },&{} 0\leqslant u<\varpi ,\\ y_u:=&{} 0, &{} \varpi \leqslant u<\varpi +{\tilde{\varpi }},\\ z_u:=&{} e^{ i\pi (-u+s)^2\varsigma /\varpi }, &{} 0\leqslant u<\varpi ,\\ z_u:=&{} e^{ i\pi (-u+s+\varpi +{\tilde{\varpi }})^2\varsigma /\varpi },&{} \varpi \leqslant u <{\tilde{\varpi }}+\varpi . \end{array} \end{aligned}$$

(101)

It can be seen that

$$\begin{aligned} G_{\varsigma }(\mathbf {x}, l+s)=e^{-i\pi (l+s)^2\varsigma /\varpi }\sum _{u=0}^{\varpi -1} y_u z_{(u-l+\varpi +{\tilde{\varpi }})\mod (\varpi +{\tilde{\varpi }})}, \quad 0\leqslant l< {\tilde{\varpi }}. \end{aligned}$$

(102)

It can be verified that the sequence $\mathbf {z}$ now satisfies the required property for a $(\varpi +{\tilde{\varpi }})$-point circular convolution. Thus it follows that $(\varpi +{\tilde{\varpi }})$-point discrete Fourier transform may be used to evaluate (102). By $F_l({\mathbf {x}})$ we denote the discrete Fourier transform of ${\mathbf {x}}:=[x_u\in \mathbb {R}:u\in \mathbb {Z}_{\varpi +{\tilde{\varpi }}}]$. That is,

$$\begin{aligned} F_l({\mathbf {x}}):=\sum _{u\in \mathbb {Z}_{\varpi +{\tilde{\varpi }}}} x_{u}e^{i\pi lu/(\varpi +{\tilde{\varpi }})}, \qquad l\in \mathbb {Z}_{\varpi +{\tilde{\varpi }}}, \end{aligned}$$

(103)

which can be computed by the fast Fourier transform. Let $F^{-1}_l({\mathbf {x}})$ denote the inverse discrete Fourier transform of ${\mathbf {x}}$. By the convolution property of the discrete Fourier transform, Eq. (102) can be rewritten as

$$\begin{aligned} G_{\varsigma }(\mathbf {x}, l+s)=e^{-i\pi (l+s)^2\varsigma } F_l^{-1}(\mathbf {w}), \quad 0\leqslant l< {\tilde{\varpi }}, \end{aligned}$$

(104)

where $\mathbf {w}$ is defined by $\mathbf {w}:=[w_k:w_k = F_k(\mathbf {y})F_k(\mathbf {z}), 0\leqslant k<{\tilde{\varpi }}+\varpi ]$.

We now present the algorithm “FrFT” for computing $G_{\varsigma }(\mathbf {x}, l+s)$ for $s\in \mathbb {Z}$ and $0 \leqslant l < {\tilde{\varpi }}$. Let $\mathbf {G}_{\mathbf {x}, \varsigma , s}:=\{G_{\varsigma }(\mathbf {x}, l): s\leqslant l \leqslant {\tilde{\varpi }} +s\}$.

Algorithm B.1

FrFT $(\mathbf {x}, \varsigma , s, {\tilde{\varpi }})$

Input: $\mathbf {x}:=\{x_u:u\in \mathbb {Z}_\varpi \}, s\in \mathbb {Z}$, $\varsigma \in \mathbb {R}$ and ${\tilde{\varpi }}\in \mathbb {N}$.

Output: $G_{\mathbf {x}, \varsigma , s}$.

Step 1. Compute $\mathbf {y}$ and $\mathbf {z}$ according to formula (101).

Step 2. Compute $\mathbf {w}:=\{w_l=F_l(\mathbf {y})F_l(\mathbf {z}): l\in \mathbb {Z}_{{\tilde{\varpi }}}\}$ by applying the fast Fourier transform [25] to $\mathbf {y}$ and $\mathbf {z}$.

Step 3. Compute $\mathbf {G}_{\mathbf {x}, \varsigma , s}$ according to (104) by applying the fast inverse Fourier transform to $\mathbf {w}$.

The total computational cost of the above algorithm is $\mathcal {O}(\varpi \log _2 \varpi )$ number of floating point operations. With the help of the fractional Fourier transform, we obtain the following algorithm for computing the summation (100). Let ${\bar{{\varvec{\xi }}}}_N:=[{\bar{\xi }}_{N, r}:r\in \mathbb {Z}_{2^Nm}]$, ${\bar{{\varvec{\xi }}}}_{N,q}:=[{\bar{\xi }}_{N, mu+q}:u\in \mathbb {Z}_{2^N}]$ and ${\varvec{\varepsilon }}_{{\varvec{\xi }}_{N},\varsigma }:=[\varepsilon _{{\varvec{\xi }}_{N},\varsigma }(l):l\in \mathbb {L}_{N,1}]$.

Algorithm B.2

FE1d $(N,{\varvec{\xi }}_{N}, \varsigma )$

Input: $N\in \mathbb {N}, {\varvec{\xi }}_N:=\{x_u\in \mathbb {C}:u\in \mathbb {Z}_{2^Nm}\}$ and $\varsigma \in \mathbb {R}$.

Output: ${\varvec{\varepsilon }}_{{\varvec{\xi }}_{N},\varsigma }$.

Step 1. Compute ${\bar{{\varvec{\xi }}}}_N$ by (97).

Step 2. For each $q\in \mathbb {Z}_m$, compute $\left( {\hat{\xi }}_{N,q,\varsigma }\right) _l$, $l\in \mathbb {L}_{{N},1}$, by using FrFT$({\bar{{\varvec{\xi }}}}_{N, q}, \varsigma , -2^N{\tilde{m}}+1,2^{N+1}{\tilde{m}}-1)$.

Step 3. For all $q \in \mathbb {Z}_m$ and $l\in \mathbb {L}_{N,1}$, compute ${\tilde{t}}_{N,q,\varsigma }(l)$ according to formula (98).

Step 3. For all $l\in \mathbb {L}_{N,1}$, compute $\varepsilon _{{\varvec{\xi }}_{N},\varsigma }(l)$ according to formula (100).

The total computational cost of the above algorithm is $\mathcal {O}(n\log _2 n)$ number of floating point operations, where $n:=2^N$.