On fast Fourier solvers for the tensor product high-order FEM for a generalized Poisson equation

TL;DR

This paper introduces efficient Fourier-based algorithms for high-order finite element methods on rectangular domains, enabling fast solutions to Poisson and related PDEs with optimal theoretical and practical performance.

Contribution

It develops novel FFT-based algorithms for eigenvector expansions in high-order FEM, improving computational efficiency for multi-dimensional PDEs.

Findings

Algorithms are logarithmically optimal in theory.

Numerical experiments confirm practical efficiency in 2D and 3D.

Applicable to various time-dependent PDEs.

Abstract

We present direct logarithmically optimal in theory and fast in practice algorithms to implement the tensor product high order finite element method on multi-dimensional rectangular parallelepipeds for solving PDEs of the Poisson kind. They are based on the well-known Fourier approaches. The key new points are the fast direct and inverse FFT-based algorithms for expansion in eigenvectors of the 1D eigenvalue problems for the high order FEM. The algorithms can further be used for numerous applications, in particular, to implement the tensor product high order finite element methods for various time-dependent PDEs. Results of numerical experiments in 2D and 3D cases are presented.

Tables

Table 1. Table 1: The 2D case. The errors in the mesh uniform norm and their ratios R C subscript 𝑅 𝐶 R_{C} in dependence on K = 2 , 4 , … , 1024 𝐾 2 4 … 1024 K=2,4,...,1024 and n = 1 , 9 ¯ 𝑛 ¯ 1 9 n=\overline{1,9} for algorithm (a)

$K$		$n=1$		$R_C$		$n=2$		$R_C$		$n=3$		$R_C$		$n=4$		$R_C$		$n=5$		$R_C$
$2$		$5$	$.1\cdot {10}^{-2}$	–		$2$	$.4\cdot {10}^{-1}$	–		$8$	$.7\cdot {10}^{-2}$	–		$3$	$.7\cdot {10}^{-2}$	–		$6$	$.1\cdot {10}^{-3}$	–
$4$		$3$	$.8\cdot {10}^{-1}$	–		$2$	$.5\cdot {10}^{-2}$	$9$	$.6$	$8$	$.4\cdot {10}^{-3}$	$10$	$.3$	$1$	$.2\cdot {10}^{-3}$	$32$	$.2$	$2$	$.1\cdot {10}^{-4}$	$29$	$.3$
$8$		$1$	$.0\cdot {10}^{-1}$	$3$	$.8$	$1$	$.6\cdot {10}^{-3}$	$16$	$.1$	$5$	$.9\cdot {10}^{-4}$	$14$	$.2$	$4$	$.7\cdot {10}^{-5}$	$24$	$.5$	$3$	$.3\cdot {10}^{-6}$	$63$	$.8$
$16$		$2$	$.6\cdot {10}^{-2}$	$3$	$.9$	$1$	$.0\cdot {10}^{-4}$	$15$	$.7$	$4$	$.1\cdot {10}^{-5}$	$14$	$.6$	$1$	$.6\cdot {10}^{-6}$	$29$	$.3$	$5$	$.4\cdot {10}^{-8}$	$60$	$.8$
$32$		$6$	$.6\cdot {10}^{-3}$	$3$	$.9$	$6$	$.2\cdot {10}^{-6}$	$16$	$.1$	$2$	$.6\cdot {10}^{-6}$	$15$	$.6$	$5$	$.2\cdot {10}^{-8}$	$31$	$.1$	$8$	$.5\cdot {10}^{-10}$	$63$	$.2$
$64$		$1$	$.6\cdot {10}^{-3}$	$4$	$.0$	$3$	$.9\cdot {10}^{-7}$	$15$	$.9$	$1$	$.6\cdot {10}^{-7}$	$15$	$.9$	$1$	$.7\cdot {10}^{-9}$	$30$	$.6$	$1$	$.3\cdot {10}^{-11}$	$64$	$.0$
$128$		$4$	$.1\cdot {10}^{-4}$	$4$	$.0$	$2$	$.4\cdot {10}^{-8}$	$16$	$.0$	$1$	$.0\cdot {10}^{-8}$	$16$	$.0$	$5$	$.4\cdot {10}^{-11}$	$31$	$.4$	$2$	$.1\cdot {10}^{-13}$	$63$	$.1$
$256$		$1$	$.0\cdot {10}^{-4}$	$4$	$.0$	$1$	$.5\cdot {10}^{-9}$	$16$	$.0$	$6$	$.4\cdot {10}^{-10}$	$16$	$.0$	$1$	$.7\cdot {10}^{-12}$	$31$	$.7$	$6$	$.4\cdot {10}^{-15}$	$32$	$.8$
$512$		$2$	$.6\cdot {10}^{-5}$	$4$	$.0$	$9$	$.6\cdot {10}^{-11}$	$16$	$.0$	$4$	$.0\cdot {10}^{-11}$	$16$	$.0$	$5$	$.4\cdot {10}^{-14}$	$31$	$.7$	$4$	$.7\cdot {10}^{-15}$	–
$\mathrm{1\hspace{0.17em}024}$		$6$	$.4\cdot {10}^{-6}$	$4$	$.0$	$6$	$.0\cdot {10}^{-12}$	$16$	$.0$	$2$	$.5\cdot {10}^{-12}$	$16$	$.0$	$2$	$.9\cdot {10}^{-15}$	$18$	$.6$	$3$	$.6\cdot {10}^{-15}$	–

Table 2. Table 2: The 2D case. The ratios of the execution times for algorithm (a) in dependence on K = 2 , 4 , … , 1024 𝐾 2 4 … 1024 K=2,4,...,1024 and n = 1 , 9 ¯ 𝑛 ¯ 1 9 n=\overline{1,9}

$K$		$n=1$		$n=2$		$n=3$		$n=4$		$n=5$		$n=6$		$n=7$		$n=8$		$n=9$
$8$		–		–		–		$1$	$.46$	$1$	$.6$	$1$	$.97$	$2$		$2$	$.02$	$2$	$.01$
$16$		$1$	$.15$	$1$	$.55$	$1$	$.8$	$2$	$.07$	$2$	$.1$	$2$	$.1$	$2$	$.13$	$2$	$.14$	$2$	$.15$
$32$		$1$	$.74$	$2$	$.24$	$2$	$.31$	$2$	$.35$	$2$	$.35$	$2$	$.35$	$2$	$.37$	$2$	$.42$	$2$	$.49$
$64$		$1$	$.99$	$2$	$.29$	$2$	$.37$	$2$	$.38$	$2$	$.41$	$2$	$.45$	$2$	$.55$	$2$	$.57$	$2$	$.81$
$128$		$2$	$.37$	$2$	$.44$	$2$	$.51$	$2$	$.57$	$2$	$.62$	$2$	$.69$	$3$	$.24$	$3$	$.19$	$3$	$.1$
$256$		$2$	$.46$	$2$	$.49$	$2$	$.46$	$2$	$.53$	$3$	$.07$	$3$	$.08$	$2$	$.9$	$3$	$.08$	$3$	$.05$
$512$		$2$	$.49$	$2$	$.66$	$2$	$.66$	$3$	$.33$	$3$	$.02$	$3$	$.15$	$2$	$.99$	$3$	$.08$	$3$	$.11$
$\mathrm{1\hspace{0.17em}024}$		$2$	$.52$	$2$	$.69$	$3$	$.33$	$3$	$.04$	$3$	$.07$	$3$	$.13$	$3$	$.14$	$3$	$.22$	$3$	$.34$

Table 3. Table 3: The 2D case. The ratios of the execution times for algorithm (b) in dependence on K = 2 , 4 , … , 1024 𝐾 2 4 … 1024 K=2,4,...,1024 and n = 1 , 9 ¯ 𝑛 ¯ 1 9 n=\overline{1,9}

$K$		$n=1$		$n=2$		$n=3$		$n=4$		$n=5$		$n=6$		$n=7$		$n=8$		$n=9$
$8$		–		–		–		$1$	$.17$	$1$	$.28$	$1$	$.97$	$1$	$.99$	$1$	$.99$	$2$	$.01$
$16$		–		$1$	$.39$	$1$	$.64$	$2$	$.02$	$2$	$.09$	$2$	$.08$	$2$	$.13$	$2$	$.2$	$2$	$.22$
$32$		$2$	$.04$	$1$	$.84$	$2$	$.17$	$2$	$.38$	$2$	$.47$	$2$	$.54$	$2$	$.49$	$2$	$.53$	$2$	$.61$
$64$		$1$	$.57$	$2$	$.26$	$2$	$.42$	$2$	$.42$	$2$	$.44$	$2$	$.52$	$2$	$.67$	$2$	$.74$	$2$	$.88$
$128$		$2$	$.34$	$2$	$.49$	$2$	$.52$	$2$	$.67$	$2$	$.77$	$2$	$.85$	$3$	$.21$	$3$	$.21$	$3$	$.27$
$256$		$2$	$.5$	$2$	$.6$	$2$	$.64$	$2$	$.8$	$3$	$.13$	$3$	$.27$	$3$	$.25$	$3$	$.34$	$3$	$.44$
$512$		$2$	$.58$	$2$	$.76$	$2$	$.91$	$3$	$.31$	$3$	$.37$	$3$	$.29$	$3$	$.31$	$3$	$.4$	$3$	$.43$
$\mathrm{1\hspace{0.17em}024}$		$2$	$.68$	$2$	$.98$	$3$	$.36$	$3$	$.36$	$3$	$.33$	$3$	$.53$	$3$	$.58$	$3$	$.84$	$3$	$.94$

Table 4. Table 4: The 3D case. The errors in the mesh uniform norm and their ratios R C subscript 𝑅 𝐶 R_{C} in dependence on K = 2 , 4 , … , 64 𝐾 2 4 … 64 K=2,4,...,64 and n = 1 , 9 ¯ 𝑛 ¯ 1 9 n=\overline{1,9} for algorithm (a)

$K$		$n=1$		$R_C$		$n=2$		$R_C$		$n=3$		$R_C$		$n=4$		$R_C$		$n=5$		$R_C$
$2$		$1$	$.3\cdot {10}^{-2}$	–		$2$	$.6\cdot {10}^{-2}$	–		$2$	$.8\cdot {10}^{-1}$	–		$2$	$.1\cdot {10}^{-1}$	–		$3$	$.1\cdot {10}^{-2}$	–
$4$		$3$	$.1\cdot {10}^{-2}$	–		$6$	$.9\cdot {10}^{-2}$	–		$3$	$.7\cdot {10}^{-2}$	$7$	$.6$	$3$	$.8\cdot {10}^{-3}$	$54$	$.6$	$1$	$.7\cdot {10}^{-3}$	$18$	$.5$
$8$		$5$	$.0\cdot {10}^{-1}$	–		$1$	$.5\cdot {10}^{-2}$	$4$	$.7$	$3$	$.1\cdot {10}^{-3}$	$12$	$.2$	$3$	$.0\cdot {10}^{-4}$	$12$	$.5$	$2$	$.9\cdot {10}^{-5}$	$57$	$.8$
$16$		$1$	$.2\cdot {10}^{-1}$	$4$	$.2$	$8$	$.4\cdot {10}^{-4}$	$17$	$.4$	$2$	$.3\cdot {10}^{-4}$	$13$	$.4$	$1$	$.1\cdot {10}^{-5}$	$27$	$.7$	$5$	$.1\cdot {10}^{-7}$	$56$	$.8$
$32$		$3$	$.0\cdot {10}^{-2}$	$4$	$.0$	$5$	$.1\cdot {10}^{-5}$	$16$	$.4$	$1$	$.5\cdot {10}^{-5}$	$15$	$.5$	$3$	$.6\cdot {10}^{-7}$	$30$	$.6$	$8$	$.3\cdot {10}^{-9}$	$62$	$.0$
$64$		$7$	$.5\cdot {10}^{-3}$	$4$	$.0$	$3$	$.2\cdot {10}^{-6}$	$16$	$.0$	$9$	$.2\cdot {10}^{-7}$	$16$	$.0$	$1$	$.2\cdot {10}^{-8}$	$31$	$.1$	$1$	$.3\cdot {10}^{-10}$	$64$	$.0$

Table 5. Table 5: The 3D case. The ratios of the execution times for algorithm (a) in dependence on K = 2 , 4 , … , 64 𝐾 2 4 … 64 K=2,4,...,64 and n = 1 , 9 ¯ 𝑛 ¯ 1 9 n=\overline{1,9}

$K$		$n=1$		$n=2$		$n=3$		$n=4$		$n=5$		$n=6$		$n=7$		$n=8$		$n=9$
$4$		$1$	$.68$	$2$	$.05$	$2$	$.34$	$4$	$.51$	$4$	$.35$	$4$	$.62$	$4$	$.58$	$4$	$.52$	$4$	$.49$
$8$		$1$	$.15$	$2$	$.6$	$3$	$.82$	$4$	$.23$	$4$	$.23$	$4$	$.25$	$4$	$.22$	$4$	$.29$	$4$	$.24$
$16$		$4$	$.21$	$4$	$.47$	$4$	$.31$	$4$	$.22$	$4$	$.29$	$4$	$.27$	$4$	$.29$	$4$	$.3$	$4$	$.31$
$32$		$4$	$.85$	$4$	$.79$	$4$	$.8$	$4$	$.89$	$4$	$.78$	$4$	$.84$	$4$	$.83$	$4$	$.86$	$5$	$.03$
$64$		$4$	$.32$	$4$	$.46$	$4$	$.78$	$4$	$.77$	$4$	$.88$	$4$	$.93$	$5$	$.05$	$5$	$.16$	$5$	$.52$

Table 6. Table 6: The 3D case. The ratios of the execution times for algorithm (b) in dependence on K = 2 , 4 , … , 64 𝐾 2 4 … 64 K=2,4,...,64 and n = 1 , 9 ¯ 𝑛 ¯ 1 9 n=\overline{1,9}

$K$		$n=1$		$n=2$		$n=3$		$n=4$		$n=5$		$n=6$		$n=7$		$n=8$		$n=9$
$4$		–		$1$	$.9$	$2$	$.09$	$4$	$.24$	$4$	$.32$	$4$	$.51$	$4$	$.47$	$4$	$.51$	$4$	$.41$
$8$		–		$2$	$.29$	$3$	$.6$	$4$	$.19$	$4$	$.25$	$4$	$.23$	$4$	$.23$	$4$	$.26$	$4$	$.25$
$16$		$3$	$.65$	$4$	$.28$	$4$	$.29$	$4$	$.3$	$4$	$.32$	$4$	$.35$	$4$	$.36$	$4$	$.44$	$4$	$.46$
$32$		$4$	$.89$	$4$	$.75$	$4$	$.78$	$4$	$.81$	$4$	$.81$	$4$	$.96$	$4$	$.94$	$5$		$5$	$.22$
$64$		$4$	$.43$	$4$	$.68$	$4$	$.87$	$4$	$.89$	$5$	$.01$	$5$	$.08$	$5$	$.31$	$5$	$.47$	$5$	$.71$

Equations

$$-u^{\prime \prime }(x)=\lambda u(x)\text{on}[0,X],u(0)=u(X)=0,u(x)\not\equiv 0.$$

$$\mathcal{A}v=\lambda \mathcal{C}v,v\in S_K^{(n)},v\ne 0.$$

$$A_{kl}={\int }_{{\sigma }_0}{e}_k^{\prime }(x)e_l^{\prime }(x)dx,C_{kl}={\int }_{{\sigma }_0}{e}_k(x)e_l(x)dx,$$

$$A=\left(\begin{array}[]{ccc}a_{0}&a^{T}&a_{n}\\ a&\widetilde{A}&\check{a}\\ a_{n}&\check{a}^{T}&a_{0}\end{array}\right),\ C=\left(\begin{array}[]{ccc}c_{0}&c^{T}&c_{n}\\ c&\widetilde{C}&\check{c}\\ c_{n}&\check{c}^{T}&c_{0}\end{array}\right),\ G(\lambda)=\left(\begin{array}[]{ccc}g_{0}(\lambda)&g^{T}(\lambda)&g_{n}(\lambda)\\ g(\lambda)&\widetilde{G}(\lambda)&\check{g}(\lambda)\\ g_{n}(\lambda)&\check{g}^{T}(\lambda)&g_{0}(\lambda)\end{array}\right).$$

$$p=p_e+p_o,p_e:=0.5(p+\check{p}),p_o:=0.5(p-\check{p}).$$

$$\check{p}\cdot q=p\cdot \check{q},\check{p}\cdot \check{q}=p\cdot q\text{for any}p,q\in {\mathbb{R}}^{n-1}.$$

$$g_n(\lambda )v_{j-1}+\check{g}(\lambda )\cdot v_{j-1/2}+2g_0(\lambda )v_j+g(\lambda )\cdot v_{j+1/2}+g_n(\lambda )v_{j+1}=0,j=\overline{1,K-1},$$

$$g(\lambda )v_{j-1}+\tilde{G}(\lambda )v_{j-1/2}+\check{g}(\lambda )v_j=0,j=\overline{1,K},$$

$$Ae=\lambda Ce,e\in {\mathbb{R}}^{n+1},e\ne 0;$$

$$\tilde{A}e=\lambda \tilde{C}e,e\in {\mathbb{R}}^{n-1},e\ne 0,$$

$$BP=PB,$$

$$a=\sum _{l=1}^{n-1}{a}^{(l)}\tilde{C}{e}^{(l)},c=\sum _{l=1}^{n-1}{c}^{(l)}\tilde{C}{e}^{(l)},\text{with}{a}^{(l)}=a\cdot e^{(l)},c^{(l)}=c\cdot e^{(l)}.$$

$$p=\sum _{l=1}^{n-1}\frac{a^{(l)}-\lambda c^{(l)}}{\lambda -{\lambda }_0^{(l)}}{e}^{(l)}=\sum _{l=1}^{n-1}\frac{a^{(l)}-{\lambda }_0^{(l)}{c}^{(l)}}{\lambda -{\lambda }_0^{(l)}}{e}^{(l)}-{\tilde{C}}^{-1}c.$$

$$\widetilde{G}(\lambda)p=\sum_{l=1}^{n-1}\big{(}\lambda_{0}^{(l)}-\lambda\big{)}p_{l}\widetilde{C}e^{(l)},\ \ g(\lambda)=\sum_{l=1}^{n-1}\big{(}a^{(l)}-\lambda c^{(l)}\big{)}\widetilde{C}e^{(l)},$$

$$(A)\text{the eigenvalues in}{S}_n\text{and}{\tilde{S}}_n\text{are simple and}{S}_n\cap {\tilde{S}}_n=\varnothing .$$

$$\hat{\gamma }(\lambda )\equiv -(g_0-g\cdot {\tilde{G}}^{-1}g)(\lambda )/(g_n-\check{g}\cdot {\tilde{G}}^{-1}g)(\lambda )=\theta ,$$

$$\displaystyle a_{0}-\lambda c_{0}+\sum_{l=1}^{n-1}\frac{(a^{(l)}-\lambda c^{(l)})^{2}}{\lambda-\lambda_{0}^{(l)}}=-\theta\Big{(}a_{n}-\lambda c_{n}+\sum_{l=1}^{n-1}\frac{(\check{a}^{(l)}-\lambda\check{c}^{(l)})(a^{(l)}-\lambda c^{(l)})}{\lambda-\lambda_{0}^{(l)}}\Big{)}.$$

$$(y,v{)}_{S_K^{(n)}}:=\sum _{j=1}^{K-1}{y}_j{v}_j+\sum _{j=1}^K{y}_{j-1/2}\cdot v_{j-1/2},\Vert v{\Vert }_{\mathcal{C}}^2:=(\mathcal{C}v,v{)}_{S_K^{(n)}}.$$

$$s_{0,j}^{(l)}=0,j=\overline{1,K-1},s_{0,j-1/2}^{(l)}=(-P{)}^{j-1}{e}^{(l)},j=\overline{1,K},$$

$$s_{k,j}^{(l)}=s_{k,j},j=\overline{1,K-1},s_{k,j-1/2}^{(l)}=p_k^{(l)}\sin \frac{\pi k(j-1)}{K}+{\check{p}}_k^{(l)}\sin \frac{\pi kj}{K},j=\overline{1,K},$$

$$\displaystyle\big{(}\mathcal{C}s_{k}^{(l)},s_{\tilde{k}}^{(\tilde{l})}\big{)}_{S_{K}^{(n)}}=0$$

$$w=\sum _{l=1}^{n-1}{w}_{0l}{s}_0^{(l)}+\sum _{k=1}^{K-1}\sum _{l=1}^n{w}_{kl}{s}_k^{(l)}.$$

$$0=v_{j-1/2}\cdot\widetilde{G}(\lambda)e=\widetilde{G}(\lambda)v_{j-1/2}\cdot e=-\big{[}(g(\lambda)\cdot e)v_{j-1}+(\check{g}(\lambda)\cdot e)v_{j}\big{]}.$$

$$\displaystyle G(\lambda)\left(\begin{array}[]{c}0\\[2.84526pt] e\\[2.84526pt] 0\end{array}\right)=\left(\begin{array}[]{c}g(\lambda)\cdot e\\[2.84526pt] 0\\[2.84526pt] \check{g}(\lambda)\cdot e\end{array}\right).$$

$$v_j=0,j=\overline{0,K}.$$

$$c_{j-1/2}\check{g}(\lambda )\cdot e+c_{j+1/2}g(\lambda )\cdot e=0,j=\overline{1,K-1}.$$

$$v_{j-1/2}=(-1{)}^{j-1}e\text{for even}e,v_{j-1/2}=e\text{for odd}e,j=\overline{1,K}$$

$$v_{j-1/2}=-G^{-1}(\lambda )[v_{j-1}g(\lambda )+v_j\check{g}(\lambda )],j=\overline{1,K}.$$

$${\hat{g}}_n(\lambda )v_{j-1}+2{\hat{g}}_0(\lambda )v_j+{\hat{g}}_n(\lambda )v_{j+1}=0,j=\overline{1,K-1},$$

$${\hat{g}}_0(\lambda )=(g_0-g\cdot {\tilde{G}}^{-1}g)(\lambda ),{\hat{g}}_n(\lambda )=(g_n-\check{g}\cdot {\tilde{G}}^{-1}g)(\lambda ).$$

$$\displaystyle G(\lambda)\left(\begin{array}[]{c}1\\[2.84526pt] -(\widetilde{G}^{-1}g)(\lambda)\\[2.84526pt] 0\end{array}\right)=\left(\begin{array}[]{c}\hat{g}_{0}(\lambda)\\[2.84526pt] 0\\[2.84526pt] \hat{g}_{n}(\lambda)\end{array}\right),\ \ G(\lambda)\left(\begin{array}[]{c}0\\[2.84526pt] -(\widetilde{G}^{-1}\check{g})(\lambda)\\[2.84526pt] 1\end{array}\right)=\left(\begin{array}[]{c}\hat{g}_{n}(\lambda)\\[2.84526pt] 0\\[2.84526pt] \hat{g}_{0}(\lambda)\end{array}\right).$$

Full text

Keywords. Fast direct algorithm, high order finite element method, FFT, Poisson equation.

MSC subject classifications. 65F05, 65F15, 65M60, 65T99.

1 Introduction

We present direct fast algorithms to implement $n$th order ($n\geqslant 2$) finite element method (FEM) on rectangular parallelepipeds [4] for solving a $N$-dimensional generalized Poisson equation ($N\geqslant 2$) with the Dirichlet boundary condition. The algorithms are based on the well-known Fourier approaches, e.g., see [2, 10, 11, 12] and references therein. The key new points are the fast direct and inverse algorithms for expansion in eigenvectors of the 1D eigenvalue problems for the high order FEM utilizing several versions of the discrete fast Fourier transform (FFT) [3]. This solves the old known problem, see [2, p. 271], and makes the full algorithms logarithmically optimal with respect to the number of elements as in the case of the bilinear elements ($n=1$) or standard finite-difference schemes. The algorithms are fast in practice (faster than the theoretical expectations) and demonstrate only a mild growth in $n$ starting from the standard case $n=1$. For example, in the 9th order case, the 2D FEM system for $2^{20}$ elements containing almost $85\cdot 10^{6}$ unknowns and the 3D FEM system for $2^{18}$ elements containing more than $190\cdot 10^{6}$ unknowns are solved respectively in less than 2 and 15 min on an ordinary laptop using Matlab R2016a code, see details below.

The algorithms can further serve for a variety of applications including general 2nd order elliptic equations (as preconditioners), for the $N$-dimensional heat, wave or time-dependent Schrödinger PDEs, etc. They can be applied for some non-rectangular domains, in particular, by involving meshes topologically equivalent to rectangular ones [7]. Other standard boundary conditions can be covered as well, see a brief description in [14]. Moreover, the Fourier structure of algorithms is especially valuable for solving some wave physics problems, in particular, involving non-local boundary conditions, e.g., see [2, 6, 13], whence our own interest arose. Clearly the algorithms are also highly parallelizable.

The paper is organized as follows. In Section 2, the statement of the 1D FEM eigenvalue problem together with auxiliary FEM eigenvalue problems on and inside the reference element are given. The basic Section 3 is devoted to a description of eigenvalues and eigenvectors of the 1D FEM eigenvalue problem and the fast direct and inverse algorithms for expansion in these eigenvectors. Applications to the generalized Poisson equation in a $N$-dimensional rectangular parallelepiped with the Dirichlet boundary condition are described in Section 4. Results of numerical experiments for $N=2$ and 3 are presented in detail in Section 5; all of them include the standard case $n=1$ for comparison.

2 The statement of 1D FEM eigenvalue problem

We first consider in detail the FEM for the simplest 1D eigenvalue ODE problem

[TABLE]

We take the uniform mesh $\bar{\omega}_{h}$ with the nodes $x_{j}=jh$, $j=\overline{0,K}$ (i.e., $0\leqslant j\leqslant K$) and the step $h=X/K$. Let $H_{h}^{(n)}[0,X]$ be the FEM space of the piecewise-polynomial functions $\varphi\in C[0,X]$ such that $\varphi(x)\in\mathcal{P}_{n}$ for $x\in[x_{j-1},x_{j}]$, $j=\overline{1,K}$, with $\varphi(0)=\varphi(X)=0$; here $\mathcal{P}_{n}$ is the space of polynomials having at most $n$th degree, $n\geqslant 2$.

Let $S_{K}^{(n)}$ be the space of vector functions $w$ such that $w_{j}\in\mathbb{R}$ for $j=\overline{0,K}$ with $w_{0}=w_{K}=0$ and $w_{j-1/2}\in\mathbb{R}^{n-1}$, $j=\overline{1,K}$. Clearly $\dim S_{K}^{(n)}=nK-1$. A function $\varphi\in H_{h}^{(n)}[0,X]$ is uniquely defined by its values at the mesh nodes $\varphi_{j}=\varphi(x_{j})$, $j=\overline{0,K}$, with $\varphi_{0}=\varphi_{K}=0$, and inside the elements $\varphi_{j-1/2}=\{\varphi(x_{j-1}+(l/n)h)\}_{l=1}^{n-1}$, $j=\overline{1,K}$, that form the element in $S_{K}^{(n)}$.

We use the following scaled operator form of the standard FEM discretization for problem (2.1)

[TABLE]

Here $\mathcal{A}=\mathcal{A}^{T}>0$ and $\mathcal{C}=\mathcal{C}^{T}>0$ are the global (scaled) stiffness and mass operators (matrices) acting in $S_{K}^{(n)}$ and together with $\lambda$ independent on $h$; the true approximate eigenvalues are $\lambda_{h}=4h^{-2}\lambda$.

Let $A=\{A_{kl}\}_{k,l=0}^{n}$ and $C=\{C_{kl}\}_{k,l=0}^{n}$ be the local stiffness and mass matrices related to the reference element $\sigma_{0}=[-1,1]$ with the following entries

[TABLE]

where $\{e_{l}\}_{l=0}^{n}$ is the Lagrange basis in $\mathcal{P}_{n}$ such that e_{l}\big{(}-1+(2k)/n\big{)}=\delta_{kl}, for $k,l=\overline{0,n}$, and $\delta_{kl}$ is the Kronecker delta. The matrices $A$, $C$ and the related matrix pencil $G(\lambda):=A-\lambda C$ have the following $3\times 3$–block form

[TABLE]

Here $\widetilde{A}$, $\widetilde{C}$ and $\widetilde{G}(\lambda)=\widetilde{A}-\lambda\widetilde{C}$ are square matrices of order $n-1$ with the column vectors $a,c,g(\lambda)=a-\lambda c\in\mathbb{R}^{n-1}$ whereas $\check{p}_{l}\equiv(Pp)_{l}=p_{n-l}$, $l=\overline{1,n-1}$, for $p\in\mathbb{R}^{n-1}$. The matrices $A$, $C$ and $G(\lambda)$ are bisymmetric (i.e. symmetric with respect to the main and secondary diagonals).

Notice that $P_{ij}=\delta_{i(n-j)}$ and $P^{T}=P^{-1}=P$. Let $\mathbb{R}_{e}^{n-1}$ and $\mathbb{R}_{o}^{n-1}$ be the subspaces of even and odd vectors in $\mathbb{R}^{n-1}$, i.e. such that respectively $Pp=p$ and $Pp=-p$. The decomposition $\mathbb{R}^{n-1}=\mathbb{R}_{e}^{n-1}\oplus\mathbb{R}_{o}^{n-1}$ (for $n\geqslant 3$) is implemented by the formulas

[TABLE]

Notice that $\dim\mathbb{R}_{e}^{n-1}=[n/2]$ and $\dim\mathbb{R}_{o}^{n-1}=[(n-1)/2]$, with $\mathbb{R}_{o}^{n-1}=\{0\}$ for $n=2$. Clearly

[TABLE]

Hereafter the symbol $\cdot$ denotes the inner product of vectors in $\mathbb{R}^{n-1}$.

Then problem (2.2) can be written in the following explicit form

[TABLE]

with $v_{0}=v_{K}=0$, $v\not\equiv 0$.

We also consider the auxiliary eigenvalue problems on and inside the reference element $\sigma_{0}$

[TABLE]

where clearly $A\geqslant 0$, $C>0$ and $\widetilde{A}=\widetilde{A}^{T}>0$, $\widetilde{C}=\widetilde{C}^{T}>0$; see some their properties in [13] (where the problem similar to (2.6), (2.7) on the uniform mesh on $[0,\infty)$ for $\lambda\in\mathbb{C}$ was studied). Denote by $S_{n}$ and $\tilde{S}_{n}$ their spectra. Let $\{\lambda_{0}^{(l)},e^{(l)}\}_{l=1}^{n-1}$ be eigenpairs of problem (2.9).

Lemma 2.1**.**

1. The subspaces $\mathbb{R}_{e}^{n-1}$ and $\mathbb{R}_{e}^{n-1}$ are invariant with respect to $\widetilde{A}$ and $\widetilde{C}$. Thus each eigenvector in $\{e^{(l)}\}_{l=1}^{n-1}$ can be chosen either even or odd. Also $\lambda_{0}^{(l)}>0$, $l=\overline{1,n-1}$.

2. Similar properties are valid for problem (2.8) with the exception of one simple zero eigenvalue.

Proof.

Any bisymmetric matrix $B$ of the order $n-1$ commutes with $P$, i.e.

[TABLE]

that implies the main result of Item 1. The property $\lambda_{0}^{(l)}>0$, $l=\overline{1,n-1}$ is well-known.

For Item 2, the argument is similar taking into account that $A(1\ldots 1)^{T}=0$ (concerning simplicity of $\lambda=0$, see Proposition 5 in [13]). ∎

One can check by the direct computation that all the eigenvalues in $S_{n}$ and $\tilde{S}_{n}$ are simple and $S_{n}\cap\tilde{S}_{n}=\emptyset$ at least for $1\leqslant n\leqslant 9$, see [13].

For low $n$, one can find $S_{n}$ and $\tilde{S}_{n}$ analytically (with the help of Mathematica), in particular, $\tilde{S}_{2}=\{2.5\}$, $\tilde{S}_{3}=\{2.5,10.5\}$, $\tilde{S}_{4}=\{14\pm\sqrt{133},10.5\}$ and $\tilde{S}_{5}=\{14\pm\sqrt{133},30\pm 9\sqrt{5}\}$ (repeatability of the eigenvalues is not occasional, see [13]).

We choose $\{e^{(l)}\}_{l=1}^{n-1}$ as in Lemma 2.1 using scaling $\widetilde{C}e^{(l)}\cdot e^{(l)}=1$.

Lemma 2.2**.**

Let $\widetilde{G}(\lambda)p=-g(\lambda)$, where $\lambda\not\in\tilde{S}_{n}$. Let the vectors $a$ and $c$ be expanded as

[TABLE]

See $\widetilde{G}(\lambda)$, $g(\lambda)$, $a$ and $c$ in (2.3). Then the following formulas hold

[TABLE]

Proof.

For the expansions $p=\sum_{l=1}^{n-1}p_{l}e^{(l)}$ and (2.11), we have

[TABLE]

and the result easily follows. ∎

3 Solving of the 1D FEM eigenvalue problem and related FFT-based algorithms

Below we impose the following assumption:

[TABLE]

Recall that it is valid at least for $2\leqslant n\leqslant 9$ (below in Section 5 we verify it up to $n=21$).

We introduce the auxiliary equation

[TABLE]

where $\lambda\not\in\tilde{S}_{n}$, with the parameter $\theta$. Its solving is equivalent to finding the roots of a polynomial having at most $n$th degree, see [13]. In particular, owing to Lemma 2.2 this equation can be rewritten as

[TABLE]

Here $\check{a}^{(l)}=\check{a}\cdot e^{(l)}$ and $\check{c}^{(l)}=\check{c}\cdot e^{(l)}$. Moreover, for $2\leqslant n\leqslant 9$ computations help to confirm that the vectors $e^{(l)}$ are even and odd respectively for odd and even $l$ provided that $\lambda_{0}^{(1)}<\ldots<\lambda_{0}^{(n-1)}$; therefore $\check{a}^{(l)}=(-1)^{l}a^{(l)}$ and $\check{c}^{(l)}=(-1)^{l}c^{(l)}$, $l=\overline{1,n-1}$.

We define the simplest inner product in $S_{K}^{(n)}$ and the corresponding squared $\mathcal{C}$-norm

[TABLE]

Next theorem describes eigenvalues and eigenvectors of problem (2.2).

Theorem 3.1**.**

1. The spectrum of problem (2.2) consists in $\tilde{S}_{n}$ and the numbers \big{\{}\lambda_{k}^{(l)}\big{\}}_{l=1}^{n} that are all $n$ (and all positive real) solutions to equation (3.2) with $\theta=\theta_{k}:=\cos\frac{\pi k}{K}$ for $k=\overline{1,K-1}$. The numbers \big{\{}\lambda_{k}^{(l)}\big{\}}_{l=1}^{n} differ from $\{\lambda_{0}^{(l)}\}_{l=1}^{n-1}$ and are different for fixed $k$ .

2. The following eigenvector corresponds to the eigenvalue $\lambda_{0}^{(l)}$:

[TABLE]

for $l=\overline{1,n-1}$. Here $(-P)^{j-1}e=(-1)^{j-1}e$ for even $e$, $(-P)^{j-1}e=e$ for odd $e$.

3. The following eigenvector corresponds to the eigenvalue $\lambda_{k}^{(l)}$:

[TABLE]

where $s_{k,j}:=\sin\frac{\pi kj}{K}$ and $p_{k}^{(l)}\in\mathbb{R}^{n-1}$ solves the non-degenerate algebraic system \widetilde{G}\big{(}\lambda^{(l)}_{k}\big{)}p_{k}^{(l)}=-g\big{(}\lambda^{(l)}_{k}\big{)}, for $k=\overline{1,K-1}$, $l=\overline{1,n}$.

4. The introduced eigenvectors are $\mathcal{C}$-orthogonal, i.e.

[TABLE]

for any $k,\tilde{k}\in\overline{0,K-1}$, $l\in\overline{1,n-\delta_{k0}}$ and $\tilde{l}\in\overline{1,n-\delta_{\tilde{k}0}}$ such that $k\neq\tilde{k}$ and/or $l\neq\tilde{l}$.

Consequently they form the basis in $S_{K}^{(n)}$, i.e. any $w\in S_{K}^{(n)}$ can be uniquely expanded as

[TABLE]

Proof.

1. We distinguish between two cases. Let first $\lambda\in\tilde{S}_{n}$ and $e$ satisfy (2.9). Then, for any $j=\overline{1,K}$, using equation (2.7) we get

[TABLE]

Clearly

[TABLE]

Since $\lambda\not\in S_{n}$ by assumption (A), we have $g(\lambda)\cdot e\neq 0$ or $\check{g}(\lambda)\cdot e\neq 0$. Owing to assumption (A) and Lemma 2.1, $\lambda$ is simple in $\tilde{S}_{n}$ and $e$ is either even or odd. Correspondingly either $\check{g}(\lambda)\cdot e=g(\lambda)\cdot e\neq 0$ and $v_{j}=-v_{j-1}$, or $\check{g}(\lambda)\cdot e=-g(\lambda)\cdot e\neq 0$ and $v_{j}=v_{j-1}$. Since $v_{0}=0$, in both cases we get

[TABLE]

Thus equation (2.7) is reduced to $\widetilde{G}v_{j-1/2}=0$ and implies that $v_{j-1/2}=c_{j-1/2}e$.

Now equation (2.6) is reduced to

[TABLE]

Therefore $c_{j+1/2}=-c_{j-1/2}$ for even $e$ or $c_{j+1/2}=c_{j-1/2}$ for odd $e$. Consequently the sought eigenvector $v$ satisfies (3.7) together with

[TABLE]

($v$ is defined up to a non-zero constant multiplier). Thus we come to eigenvectors (3.3).

2. Next let $\lambda\not\in\tilde{S}_{n}$. Then from equation (2.7) we get

[TABLE]

Inserting this into equation (2.6), we find the three-point equation

[TABLE]

where

[TABLE]

Straightforwardly the following equalities hold

[TABLE]

If $\hat{g}_{n}(\lambda)=\hat{g}_{0}(\lambda)=0$, then the equalities mean that $\lambda$ is at least double eigenvalue for problem (2.8) that contradicts assumption (A).

If $\hat{g}_{n}(\lambda)=0$ and $\hat{g}_{0}(\lambda)\neq 0$, then equation (3.9) together with (3.8) lead to $v=0$ thus such $\lambda$ does not satisfy (2.2).

Therefore $\hat{g}_{n}(\lambda)\neq 0$ and equation (3.9) is simplified to

[TABLE]

with the function $\hat{\gamma}(\lambda)=\hat{g}_{0}(\lambda)/\hat{g}_{n}(\lambda)$ (see it also in (3.1)). Since $v_{0}=v_{n}=0$, we can use the expansion

[TABLE]

and define the vector $\tilde{{\bf v}}:=(\tilde{v}_{1},\ldots,\tilde{v}_{K-1})$ of its coefficients. Using the expansion in (3.10) gives

[TABLE]

Clearly this equality is valid for some $\tilde{{\bf v}}\neq 0$ if and only if

[TABLE]

Notice that $\tilde{{\bf v}}=0$ is equivalent to $v=0$ in $S_{K}^{(n)}$ (taking into account formula (3.8)). Therefore $\lambda$ satisfies (3.12); moreover, $\tilde{v}_{j}=\delta_{kj}$ and consequently $v_{j}=s_{k,j}$, $j=\overline{0,K}$, together with

[TABLE]

see (3.8) (all last three equalities are valid up to the same non-zero multilplier). Thus we come to eigenvectors (3.4).

The total amount of eigenvalues $\lambda\not\in\tilde{S}_{n}$ (taking into account their possible multiplicity) is $\dim S_{K}^{(n)}-(n-1)=n(K-1)$. The maximal amount of roots algebraic equations (3.12) for all $k$ is the same so that each equation (3.12) has to possess exactly $n$ distinct roots (for fixed $k$, the written eigenvector $v$ is defined by $\lambda$ uniquely).

3. Property (3.5) is knowingly valid for eigenvectors $s_{k}^{(l)}$ and $s_{\tilde{k}}^{(\tilde{l})}$ corresponding to different eigenvalues of problem (2.2), in particular, for $k=0$ and $\tilde{k}\neq 0$, or $k=\tilde{k}$ and $l\neq\tilde{l}$. The remaining case will be covered below in Corollary 3.1 of the related Lemma 3.1. ∎

Notice that: (1) the vectors $s_{0}^{(l)}$ are used only to describe the algorithm, and only the vectors $e^{(l)}$ are applied in its implementation; (2) $s_{k,j}^{(l)}$ are independent on $l$; (3) the vectors $p_{k}^{(l)}$ are independent on $j$ and can also be computed owing to Lemma 2.2.

Lemma 3.1**.**

Let $w\in S_{K}^{(n)}$ and $w_{j-1/2}=qw_{j-1}+\check{q}w_{j}$, $j=\overline{1,K}$, for some $q\in\mathbb{R}^{n-1}$. Then

[TABLE]

for $k=\overline{1,K-1}$, $l=\overline{1,n}$, where $(s_{k},w)_{\omega_{h}}:=\sum_{j=1}^{K-1}s_{k,j}w_{j}$.

Proof.

1. For any $v,w\in S_{K}^{(n)}$, recalling notation (2.3) we have

[TABLE]

2. According to formulas (3.15) and (3.3), for even $e^{(l)}$, we get

[TABLE]

since $(c-\check{c})\cdot e^{(l)}=0$ and $\widetilde{C}e^{(l)}\cdot(q-\check{q})=0$ for any $c,q\in\mathbb{R}^{n-1}$ as well as $w_{0}=w_{K}=0$.

For odd $e^{(l)}$, we similarly get

[TABLE]

since $(c+\check{c})\cdot e^{(l)}=0$ and $\widetilde{C}e^{(l)}\cdot(q+\check{q})=0$ for any $c,q\in\mathbb{R}^{n-1}$ as well as $w_{0}=w_{K}=0$. Equality (3.13) is proved.

3. Formula (3.15) together with (3.4) and (2.5) imply that

[TABLE]

Owing to formula

[TABLE]

we first derive

[TABLE]

Second, using formulas (2.5) and (2.10), we get

[TABLE]

Owing to $s_{k,0}=s_{k,K}=0$, $w_{0}=w_{K}=0$ as well as formulas (3.17), we further derive

[TABLE]

Adding (3.18) and (3.19), we prove (3.14). ∎

Corollary 3.1**.**

The orthogonality property (3.5) from Theorem 3.1, Item 4 is valid.

Proof.

It remains to consider the case $k,\tilde{k}\in\overline{1,K-1}$ and $k\neq\tilde{k}$. Since then $(s_{k},s_{\tilde{k}})_{\omega_{h}}=0$, the result follows from (3.14). ∎

We call the calculation of $w\in S_{K}^{(n)}$ by the coefficients $w_{kl}$ of the expansion (3.6) as the inverse $F_{n}$-transform and the calculation of the coefficients $w_{kl}$ by $w\in S_{K}^{(n)}$ as the direct $F_{n}$-transform. We also consider related expansion of $y\in S_{K}^{(n)}$

[TABLE]

and the calculation of the coefficients $\widetilde{y}_{kl}$ by $y\in S_{K}^{(n)}$ that we call as the direct $FC_{n}$-transform.

Let us describe their fast FFT-based implementation.

Theorem 3.2**.**

1. The inverse $F_{n}$-transform can be implemented according to the following formulas

[TABLE]

where $d_{k,e}$ and $d_{k,o}$ are respectively even and odd components of the vector $d_{k}:=\sum_{l=1}^{n}w_{kl}p_{k}^{(l)}$. Note that $(-P)^{j-1}e=e$ for odd $j$ and $(-P)^{j-1}e=-\check{e}$ for even $j$ for any $e\in\mathbb{R}^{n-1}$.

The collection $\{w_{j}\}_{j=1}^{K-1}$ can be computed by the standard inverse FFT with respect to sines. The collection $\{w_{j-1/2}\}_{j=1}^{K}$ can be computed by $n-1$ modified inverse FFT related to the centers of elements in the amount of $[n/2]$ with respect to sines and $[(n-1)/2]$ with respect to cosines using extensions $d_{K,e}:=0$ and $d_{0,o}:=0$, see algorithms DST-I, DST-III and DCT-III in [3].

2. The direct $FC_{n}$-transform can be implemented basing on the standard formula

[TABLE]

Here, first, for $k=0$, $l=\overline{1,n-1}$, the following formulas hold

[TABLE]

Second, for $k=\overline{1,K-1}$, $l=\overline{1,n}$, the following formulas hold

[TABLE]

Notice that the sums in formula (3.25) are independent on $l$. The collection of all these coefficients can be computed using $n$ standard direct FFTs with respect to sines.

3. Similarly to Item 2, the direct $F_{n}$-transform can be implemented basing on the standard formula

[TABLE]

Here, for $k=0$, $l=\overline{1,n-1}$, the following formula holds

[TABLE]

For $k=\overline{1,K-1}$, $l=\overline{1,n}$, formula (3.25) with $y:=\mathcal{C}w$ is applicable. Alternatively, the following formula holds as well

[TABLE]

where $q_{k,e}^{(l)}$ and $q_{k,o}^{(l)}$ are respectively even and odd components of the vector $q_{k}^{(l)}:=\tilde{C}p_{k}^{(l)}+c$. Once again all these coefficients can be computed using $n$ standard direct FFTs with respect to sines.

Proof.

1. Let the coefficients $w_{kl}$ of expansion (3.6) be known. According to the first formulas (3.3) and (3.4), the values of $w$ for integer indices in (3.6) are reduced to (3.21).

To compute $w$ for half-integer indices, we transform the second sum in (3.6). Owing to decomposition (2.4) we rewrite the second formula (3.4) in the form

[TABLE]

Then using also the second formula (3.3), we obtain formula (3.22).

2. Now we consider the computation of the coefficients in expansion (3.20) for given $y\in S_{K}$. Owing to the orthogonality property (3.5), they first can be expressed in the form (3.23) for $k=0$, $l=\overline{1,n-1}$ and $k=\overline{1,K-1}$, $l=\overline{1,n}$.

Formulas (3.15) and (3.3) imply that

[TABLE]

Lemma 3.1 immediately implies formula (3.26) since $(s_{k},s_{k})_{\omega_{h}}=K/2$.

By virtue of formulas (3.3) for the numerator of formula (3.23) for $k=0$ we can write

[TABLE]

By virtue of formulas (3.4) for the same numerator for $k=\overline{1,K-1}$ we get

[TABLE]

Therefore shifting by 1 the index in the second of these sums, and applying the identity $a_{1}b_{1}+a_{2}b_{2}=0.5(a_{1}+a_{2})(b_{1}+b_{2})+0.5(a_{1}-a_{2})(b_{1}-b_{2})$ and recalling decomposition (2.4), we derive

[TABLE]

Since also

[TABLE]

we obtain formula (3.25).

3. Owing to the orthogonality property (3.5), formula (3.27) is valid.

By virtue of formulas (3.3), (3.15) as well as $\widetilde{C}^{T}=\widetilde{C}$ and $P=P^{T}$ for its numerator for $k=0$ we can write

[TABLE]

By virtue of formulas (3.16), (3.18) for the same numerator for $k=\overline{1,K-1}$ we get

[TABLE]

where $q_{k}^{(l)}=\tilde{C}p_{k}^{(l)}+c$. Transforming the last sum in the same manner as above the second and third terms in (3.30), we obtain (3.29). ∎

4 Applications to the generalized Poisson equation

To begin with, we turn to the simple 1D ODE boundary value problem

[TABLE]

where for simplicity $\alpha=\textrm{const}>-(\pi/X)^{2}$. Its FEM discretization has the operator form

[TABLE]

where $f^{h}\in S_{K}^{(n)}$ is the FEM average of $f$. Its solution can be written in the form

[TABLE]

of the expansion like (3.6), where $\widetilde{f}^{h}_{kl}$ are the coefficients of the expansion like (3.20) for the vector $f^{h}$; recall that $\delta_{k0}$ is the Kronecker delta.

Next we consider in detail solving of the $N$-dimensional ($N\geqslant 2$) boundary value problem

[TABLE]

where $\Delta$ is the Laplace operator and \alpha=\textrm{const}>-\pi^{2}\big{(}X_{1}^{-2}+\ldots+X_{N}^{-2}\big{)} (for simplicity, in order to treat the positive definite operator that actually is not so necessary).

We define the space $H_{h_{1}}^{(n_{1})}[0,X_{1}]\otimes\ldots\otimes H_{h_{N}}^{(n_{N})}[0,X_{N}]$ of the piecewise-polynomial in $\overline{\Omega}$ functions, where $h_{i}=X_{i}/K_{i}$ and $n_{i}\geqslant 2$, $i=\overline{1,N}$. Let $\mathbf{K}=(K_{1},\ldots,K_{N})$ and $\mathbf{n}=(n_{1},\ldots,n_{N})$.

We also define the space $S_{\mathbf{K}}^{(\mathbf{n})}=S_{K_{1}}^{(n_{1})}\otimes\ldots\otimes S_{K_{N}}^{(n_{N})}$ of vector functions. For example, for $N=2$, these functions are numbers for the indices $(j_{1},j_{2})$, $j_{1}=\overline{0,K_{1}}$, $j_{2}=\overline{0,K_{2}}$, and vectors from $\mathbb{R}^{n_{1}-1}$, $\mathbb{R}^{n_{2}-1}$ and $\mathbb{R}^{(n_{1}-1)\times(n_{2}-1)}$ respectively for the indices

[TABLE]

as well as zero vectors for $j_{1}=0,K_{1}$ and $j_{2}=0,K_{2}$. Similarly to the 1D case, there is the natural isomorphism between functions in $H_{h_{1}}^{(n_{1})}[0,X_{1}]\otimes\ldots\otimes H_{h_{N}}^{(n_{N})}[0,X_{N}]$ and vectors in $S_{\mathbf{K}}^{(\mathbf{n})}$.

The FEM dicretization of problem (4.4) can be written in the following operator form

[TABLE]

where $\mathcal{A}_{i}$ and $\mathcal{C}_{i}$ are versions of the above defined operators $\mathcal{A}$ and $\mathcal{C}$ acting in variable $x_{i}$ (depending on $K_{i}$ and $n_{i}$), $i=\overline{1,N}$, and $f^{h}\in S_{\mathbf{K}}^{(\mathbf{n})}$ is the FEM average of $f$. Recall that the case of the non-homogeneous Dirichlet boundary condition $u(x)=b(x)$ on $\partial\Omega$ in (4.4) could be easily covered by reducing to (4.5) with the modified $f^{h}$ depending on an approximation $b^{h}$ of $b$.

To compute its solution, the $F_{n}$-transforms from Theorem 3.2 can be applied twofold.

(a) We consider the multiple expansion of $f^{h}\in S_{\mathbf{K}}^{(\mathbf{n})}$ like (3.20)

[TABLE]

Then the expansion of the solution has the following form

[TABLE]

Here \big{\{}\lambda_{i,k_{i}}^{(l_{i})},s_{i,k_{i}}^{(l_{i})}\big{\}} are versions of the above defined eigenpairs \big{\{}\lambda_{k}^{(l)},s_{k}^{(l)}\big{\}} with respect to $x_{i}$.

Algorithm (a) comprises two rather standard steps:

(1) finding the coefficients of expansion (4.6) for $f^{h}$ (by the direct $FC_{n}$-transforms in $x_{1}$,…, $x_{N}$);

(2) finding $v$ by the coefficients of its expansion (4.7) (by the inverse $F_{n}$-transforms in $x_{1}$,…, $x_{N}$).

(b) We consider the expansion of $f^{h}$ like (3.20) in $x_{2}$,…, $x_{N}$, i.e.

[TABLE]

now with the coefficients $\widetilde{f}^{h}_{k_{2}l_{2},\ldots,k_{N}l_{N}}\in S_{K_{1}}^{(n_{1})}$. Then the coefficients $v_{kl}\in S_{K_{1}}^{(n_{1})}$ in the similar expansion of the solution $v\in S_{\mathbf{K}}^{(\mathbf{n})}$

[TABLE]

serve as the solutions to 1D problems in $x_{1}$

[TABLE]

Their matrices are symmetric and positive definite.

Algorithm (b) comprises three rather standard steps:

(1) finding the coefficients of the expansion (4.8) for $f^{h}$ (by the direct $FC_{n}$-transforms in $x_{2}$,…, $x_{N}$);

(2) solving the collection of the independent 1D problems (4.10) for the coefficients of the expansion of $v$;

(3) finding $v$ by the coefficients of its expansion (4.9) (by the inverse $F_{n}$-transforms in $x_{2}$,…, $x_{N}$).

Implementing algorithms (a) and (b) costs respectively O\big{(}K_{1}\ldots K_{N}\log_{2}(K_{1}\ldots K_{N})\big{)} and O\big{(}K_{1}\ldots K_{N}\log_{2}(K_{2}\ldots K_{N})\big{)} arithmetic operations.

Importantly, they can be applied to solve various time-dependent PDEs such as the heat, wave or Schrödinger’s equations since for their implicit time discretizations one usually gets problems like (4.5) at the upper time level.

Moreover, algorithm (b) is directly extended to the case of more general equations than in (4.4) with the coefficients depending on $x_{1}$ (that is essential, in particular, in the polar and cylindrical coordinates), various boundary conditions for $x_{1}=0,X_{1}$ and the nonuniform mesh in $x_{1}$ [11]. It can also be applied for reducing 3D problems in a cylindrical domain to a collection of independent 2D problems in the cylinder base.

5 Numerical experiments

1. We first check that the eigenvalues of each of problems (2.8) and (2.9) are well separated. We define their spectral gaps as

[TABLE]

where $S_{n}=:\{\lambda_{0(l)}\}_{l=1}^{n+1}$ and present $\tilde{\delta}_{n}$ and $\tilde{\delta}_{n}$ in Fig. 1 (left). The terms $\pi^{2}/4$ and $3\pi^{2}/4$ are the spectral gaps (in fact, the gaps between two minimal eigenvalues) of the corresponding ODE problems, see [13]. We observe that both $\delta_{n}$ and $\tilde{\delta}_{n}$ are decreasing and rapidly tend to 0 as $n$ increases. We also checked that $S_{n}\cap\tilde{S}_{n}=\emptyset$ for all $2\leqslant n\leqslant 21$.

Also we give the spectral condition numbers $\operatorname*{cond}\widetilde{A}$ and $\operatorname*{cond}\widetilde{C}$ in Fig. 1 (right) and remark their rapid growth as $n$ increases (unfortunately).

Notice that all our computations are accomplished on an ordinary notebook ASUS-U36S with Intel Core i3-2350M CPU 2.3 GHz, 8 Gb, Win 10 x64. The codes in Matlab R2016a were developed to implement the algorithms, and we emphasize that several basic and advanced [1] code vectorisation techniques were applied to speed up them notably.

In the case of the 1D ODE problem (4.2) for $\alpha=1$, we provide the condition numbers of the global FEM matrix $4h^{-2}\mathcal{A}+\alpha\mathcal{C}$ for $\alpha=1$ and its local version in dependence on $K=2,4,...,1024$ for $n=\overline{1,9}$ in Fig. 2 for a further reference. Notice their rather rapid growth as $K$ increases.

Next we analyze the execution time for $F_{n}$-transforms in dependence on the choice of $K$ as it grows up to very high values $2^{20}$. We consider three choices of $K$: powers of two ($K=2^{p}$), primes or other composite numbers (not powers of two). Primes and composite values of $K$ are randomly chosen as different ten numbers between each two consecutive powers of two, and the comparison is accomplished in the spirit of [8]. We use the external function comp_dst for DST-I and DST-III and comp_dct for DCT-III from Large Time-Frequency Analysis Toolbox [15] based on FFTW library [9], which is also used in the Matlab function fft for FFT.

The execution times for the DST-I and DST-III as well as our transforms in the case $n=5$ (for definiteness) are given respectively in Fig. 3 and Fig. 4. We do not take into account the execution time for the pairs $\{\lambda^{(l)}_{k},p_{k}^{(l)}\}$ since it becomes negligible for multiple using of the transforms for fixed $K$ as required below. For our transforms the execution times are respectively larger (than for the DST-I and DST-III) due to multiple using of FFTs and some additional computations. The inverse $F_{n}$-transform takes less time than the direct one. Moreover, the best results are mainly for $K=2^{p}$ though this is not the case for DST-I (fortunately, we apply it in the optimal case $K=2^{p}-1$). For the two other choices of $K$ the execution times are worse but close, and the difference between all of them is less than in the case of both DST-I and DST-III. These results are attractive. Notice that the data in Fig. 4 can be approximated by linear functions for $2^{5}\leqslant K\leqslant 2^{10}$ and $2^{10}\leqslant K\leqslant 2^{20}$ but with rather flat slope in the former case and a visibly higher slope in the latter one.

Remark 5.1**.**

The last phenomenon is explained by the advanced architecture of modern processors involving cache memory, streaming SIMD extensions and advanced vector extensions of the instruction set, etc. Also the above mentioned high-quality implementations of FFTs are applied.

2. Our main computational results concern solving problem (4.4) for both $N=2$ and 3, $\alpha=1$ and $X_{i}=1$. We take $K_{i}=K$ and $n_{i}=n$ for simplicity. We apply the multiple Gauss quadrature formulas with $n+1$ nodes in $x_{i}$ to compute $f^{h}$. The eigenpairs of the 1D problems are computed with the quadruple precision (using Mathematica) to improve the stability with respect to round-off errors. We notice that system (4.5) contains $(Kn-1)^{N}$ unknowns. For comparison, we include the simplest known case $n=1$ implemented in the code in the unified manner.

We first consider the 2D case ($N=2$) and take the exact solution $u(x):=\sin(2\pi x_{1})\sin(3\pi x_{2})$ \cosh\big{(}\sqrt{2}x_{1}-x_{2}\big{)}. We compute the FEM solution for different values of $n$ and $K$ in order to study the (absolute) error behaviour, see Fig. 5 and Table 1 (where values of $R_{C}$ less than 3 are omitted). To reduce the possible round-off errors, hereafter we compute the pairs $\{\lambda^{(l)}_{k},p_{k}^{(l)}\}$ with the quadruple precision. For algorithm (a), the error behaviour in the uniform mesh norm is standard: the rate of its decreasing $R_{C}$ is mainly proportional to $(n+1)^{2}$ except for $n=2$ when it is faster and similar to $n=3$ (the exception has previously been noted in [5]). Curiously, for $n=9$ and the minimal $K=2$ the error is already less than for $n=1$ and the maximal $K=1024$. Of course, the error cannot be better than the level of round-off errors that is achieved the faster, the higher $5\leqslant n\leqslant 9$. We observe that impact of the round-off errors is almost absent. For algorithm (b), the situation is similar only up to the error level \sim$$10^{-11}. Once this value is reached for fixed $3\leqslant n\leqslant 9$, we further see the error growth as $K$ increases that means the perceptible impact of the round-off errors. This is due to the respective growth of condition numbers for matrices in system (4.10) as $K$ or $n$ increases, see above Fig. 2. Consequently algorithm (a) is preferable than (b) provided that very high precision is required. Notice that if we use the double precision olny, then the level of the best error becomes at \sim$$10^{-12}-10^{-13} and the results remain stable for algorithm (a) but they remain practically unchanged for algorithm (b). Thus only the double precision computations are possible provided that the mentioned accuracy is sufficient (that is the case in a lot of applications).

We also analyze the execution time for both algorithms (using multiple program runs and their median execution time), see Fig. 6 and Tables 2 and 3 for the same $K$ and $n$. Clearly this time is independent of the above specific choice of $u$. We do not take into account the execution time for the pairs $\{\lambda^{(l)}_{k},p_{k}^{(l)}\}$ considering the case where they are computed in advance and stored (recall that they are independent of the data of PDE problem (4.4)). We call attention to the weakly superlinear behaviour of time in $K$ and its mild monotone growth in $n$. Notice that all the ratios of the consecutive execution times in the both tables are even less than the lower bound 4 for the theoretical ratios (the ratios less than one are omitted); see Remark 5.1 in this respect. In contrast to theoretical expectations, almost all ratios for algorithm (a) are less than for (b). The ratios grows as $K$ and $n$ increase, and for algorithm (b) the highest value is already close to 4. For the maximal $K=1024$ and $n=9$, system (4.5) contains almost $85\cdot 10^{6}$ unknowns but only less than 2 min is required to solve it that is the excellent result.

Finally we consider the most interesting 3D case ($N=3$) and take the exact solution $u(x):=\sin(2\pi x_{1})\sin(3\pi x_{2})$ \sin(4\pi x_{3})\cosh\big{(}\sqrt{2}x_{1}-x_{2}+x_{3}/\sqrt{3}\big{)}. Once again we compute the FEM solution for different values of $K=2,4,\ldots,64$ and $n=\overline{1,9}$ and study the error behaviour, see Fig. 7 and Table 4 (where values of $R_{C}$ less than 3 are omitted). Conclusions are generally the same as in the 2D case. Notice that now the worse stability properties of algorithm (b) are visible only for $n=8,9$ since much less maximal value of $K$ is taken.

The execution times in the 3D case are presented in Fig. 8 and Tables 5 and 6. Once more conclusions are similar to the 2D case. All the ratios in the both tables are notably less than the lower bound 8 for the theoretical ratios. Importantly, for the maximal $K=64$ and $n=9$, system (4.5) contains more than $190\cdot 10^{6}$ unknowns but only less than 15 min is required to solve it that is the nice result (especially taking into account the Matlab implementation of loops).

The study has been funded within the framework of the Academic Fund Program at the National Research University Higher School of Economics in 2016-2017 (grant no. 16-01-0054) and by the Russian Academic Excellence Project ‘5-100’ as well as by the RFBR, grant № 16-01-00048.