Jingyu Huang, Chenliang Li^*
School of Mathematics and Computing Science, Guilin University of Electronic Technology, Guilin, China.
DOI: 10.4236/ajcm.2026.162008   PDF    HTML   XML   4 Downloads   36 Views   Citations

Abstract

In this paper, based on Tucker product tensor form, we present an economic cascadic tensor multigrid method(ECTMG) for the 3D partial differential equations discretized by the finite element method. Meanwhile, we provide the convergence analysis of the new method. Finally, we verify the effectiveness of the new method through some numerical examples.

Keywords

Partial Differential Equations, Economic Cascadic Tensor Multigrid Method, Tucker Product, Finite Element Method

Share and Cite:

1. Introduction

The finite element method is widely applied in fields such as industrial equipment [1], biomedicine [2] [3], and new energy catalysis [4]. Partial differential equations are discretized into linear equation systems (1) through the finite element method [5].

$Ax=b,$ (1)

$A=\left({\displaystyle \sum }_{k=1}^N{I}_{J_N}\otimes \cdots \otimes I_{J_{k+1}}\otimes A^{\left(k\right)}\otimes I_{J_{k-1}}\otimes \cdots \otimes I_{J_1}\right)$ , $A^{\left(k\right)}\in {ℝ}^{J_k\times J_k}\left(k=1,\cdots ,N\right)$ , $A^{\left(k\right)}\otimes I_{J_{k-1}}$ represent as follows

$A^{\left(k\right)}\otimes I_{J_{k-1}}=\left[\begin{array}{ccc}{a}_{11}{I}_{J_{k-1}}& \cdots & a_{1J_2}{I}_{J_{k-1}}\\ ⋮& \ddots & ⋮\\ a_{1J_1}{I}_{J_{k-1}}& \cdots & a_{J_1{J}_2}{I}_{J_{k-1}}\end{array}\right].$

Chen and Lu [6] equivalently transform the linear system (1) into the following Sylvester tensor equation

$\mathcal{X}{\times }_1{A}^{\left(1\right)}+\mathcal{X}{\times }_2{A}^{\left(2\right)}+\cdots +\mathcal{X}{\times }_N{A}^{\left(N\right)}=\mathcal{D},$ (2)

with $A^{\left(k\right)}\in {ℝ}^{J_k\times J_k}\left(k=1,\cdots ,N\right)$ and $\mathcal{D}\in {ℝ}^{J_1\times \cdots \times J_N}$ , $\mathcal{X}\in {ℝ}^{J_1\times J_2\times \cdots \times J_N}$ . The mode-k product of $\mathcal{X}$ and $\mathcal{A}$ is denoted by $\mathcal{X}{\times }_k{A}^{\left(k\right)}$ , the result is of size $J_1\times \cdots \times J_{k-1}\times J\times J_{k+1}\times \cdots \times J_N$ , and its entries are defined by

${\left(\mathcal{X}{\times }_k{A}^{\left(k\right)}\right)}_{i_1\cdots i_{k-1}j{i}_{k+1}\cdots i_N}={\displaystyle \sum }_{i_k=1}^{J_k}{x}_{i_1\cdots i_{k-1}{i}_k{i}_{k+1}\cdots i_N}{a}_{j{i}_k}^{\left(k\right)}.$

Compared with Equation (1), Equation (2) reduces the storage space required for the solution process, thereby improving the solution efficiency. Grasedyck et al. [5] utilized the finite element method and expressed $A^{\left(i\right)}$ in Equation (2) as follows

$A^{\left(i\right)}={\left(M^{\left(i\right)}\right)}^{-1}{\tilde{A}}^{\left(i\right)}.$

Here, $M^{\left(i\right)}$ represents the mass matrix, and ${\tilde{A}}^{\left(i\right)}$ is the coefficient matrix obtained by the finite element method.

The structure of this paper is as follows. In Section 2, we discretize the high-dimensional partial differential equation by the finite element method, and give the convergence analysis of the new method. Section 3 shows and discusses the numerical results. Finally, in Section 4, we conclude this article with some content.

2. Main Result

2.1. Discretization of 3D PDE by Cubic Element

Let $\Omega \subset {ℝ}^3$ be a bounded Lipschitz domain with boundary $\partial \Omega ={\Gamma }_D\cup {\Gamma }_N$ , ${\Gamma }_D\cap {\Gamma }_N=\varnothing$ , and $\text{meas}\left({\Gamma }_D\right)>0$ . Consider the following elliptic boundary value problem:

$\{\begin{array}{ll}-\nabla \cdot \left(A\nabla u\right)+b\cdot \nabla u+c u=f\hfill & \text{in}\Omega ,\hfill \\ u=0\hfill & \text{on}{\Gamma }_D,\hfill \\ \left(A\nabla u\right)\cdot n=g\hfill & \text{on}{\Gamma }_N,\hfill \end{array}$

where

• $A\left(x\right)\in {ℝ}^{9\times 9}$ is symmetric and positive definite: there exist constants $0<\alpha \le \beta <\infty$ such that

$\alpha {\Vert \xi \Vert }^2\le {\xi }^{\text{T}}A\left(x\right)\xi \le \beta {\Vert \xi \Vert }^2\forall \xi \in {ℝ}^3,\text{a}.\text{e}.x\in \Omega ;$

• $b\left(x\right)\in {ℝ}^3$ is a bounded vector field; $c\left(x\right)\ge 0$ is a bounded scalar;

• $f\in L^2\left(\Omega \right)$ , $g\in L^2\left({\Gamma }_N\right)$ are given data; $n$ is the outward unit normal on ${\Gamma }_N$ .

Define the function space $V=\left\{v\in H^1\left(\Omega \right):{v|}_{{\Gamma }_D}=0\right\}$ . Multiplying the PDE by a test function $v\in V$ and integrating by parts yields the weak problem: find $u\in V$ such that

$a\left(u,v\right)=\ell \left(v\right)\forall v\in V,$

with

$\{\begin{array}{l}a\left(u,v\right)={\displaystyle {\int }_{\Omega }\left(\left(A\nabla u\right)\cdot \nabla v+\left(b\cdot \nabla u\right)v+cuv\right)\text{d}x},\hfill \\ \ell \left(v\right)={\displaystyle {\int }_{\Omega }fv\text{d}x}+{\displaystyle {\int }_{{\Gamma }_N}gv\text{d}s}.\hfill \end{array}$

Then its variational form is as follows

$-{\displaystyle {\iiint }_V\left[\nabla \cdot \left(A\nabla u\right)\right]v\text{d}V}+{\displaystyle {\iiint }_V\left(b\cdot \nabla u\right)v\text{d}V}+{\displaystyle {\iiint }_{V}cuv\text{d}V}={\displaystyle {\iiint }_{V}fv\text{d}V}.$

Thus, the complete weak form is obtained

${\displaystyle {\iiint }_V\left(A\nabla u\cdot \nabla v+\left(b\cdot \nabla u\right)v+cuv\right)\text{d}V}-{\displaystyle {\iint }_{\partial V}v\frac{\partial u}{\partial n_A}\text{d}S}={\displaystyle {\iiint }_{V}fv\text{d}V}.$

Let

$N_x^{\left(l\right)},N_y^{\left(l\right)},N_z^{\left(l\right)}=2^{l+1},$

$h=\frac{1}{N_i^{\left(l\right)}+1}\left(i=x,y,z\right)$ ,

it can be obtained that

$\Delta x=\frac{x_L-x_0}{N_x^{\left(l\right)}},\Delta y=\frac{y_L-y_0}{N_{ya}^{\left(l\right)}},\Delta z=\frac{z_L-z_0}{N_z^{\left(l\right)}}.$

Figure 1 presents the cube element.

Figure 1. Cube element.

We write the following basic functions,

$\begin{array}{ll}{\widehat{L}}_1=\xi \left(1-\eta \right)\left(1-\gamma \right),\hfill & {\widehat{L}}_2=\xi \eta \left(1-\gamma \right),\hfill \\ {\widehat{L}}_3=\xi \left(1-\eta \right)\gamma ,\hfill & {\widehat{L}}_4=\xi \eta \gamma ,\hfill \\ {\widehat{L}}_5=\left(1-\xi \right)\left(1-\eta \right)\left(1-\gamma \right),\hfill & {\widehat{L}}_6=\left(1-\xi \right)\eta \left(1-\gamma \right),\hfill \\ {\widehat{L}}_7=\left(1-\xi \right)\left(1-\eta \right)\gamma ,\hfill & {\widehat{L}}_8=\left(1-\xi \right)\eta \gamma .\hfill \end{array}$

And the coordinates are mapped as

$\begin{array}{l}x={\displaystyle \sum }_{i=1}^8{x}_i{\widehat{L}}_i\left(\xi ,\eta ,\gamma \right),\\ y={\displaystyle \sum }_{i=1}^8{y}_i{\widehat{L}}_i\left(\xi ,\eta ,\gamma \right),\\ z={\displaystyle \sum }_{i=1}^8{z}_i{\widehat{L}}_i\left(\xi ,\eta ,\gamma \right).\end{array}$

Furthermore, it can be concluded that the Jacobian matrix $J$ is

$J=\left|\begin{array}{lll}\frac{\partial x}{\partial \xi }\hfill & \frac{\partial y}{\partial \xi }\hfill & \frac{\partial z}{\partial \xi }\hfill \\ \frac{\partial x}{\partial \eta }\hfill & \frac{\partial y}{\partial \eta }\hfill & \frac{\partial z}{\partial \eta }\hfill \\ \frac{\partial x}{\partial \gamma }\hfill & \frac{\partial y}{\partial \gamma }\hfill & \frac{\partial z}{\partial \gamma }\hfill \end{array}\right|.$

Therefore,

$\begin{array}{c}{K}_{ij}={\displaystyle {\iiint }_{e_1}\nabla L_i\nabla L_j\nabla L_k\text{d}x\text{d}y\text{d}z}\\ ={\displaystyle {\iiint }_{{\widehat{e}}_1}\nabla {\widehat{L}}_i\nabla {\widehat{L}}_j\nabla {\widehat{L}}_k\left|J\right|\text{d}\xi \text{d}\eta \text{d}\gamma }.\end{array}$

After assembling the stiffness matrix of the assembly unit, the linear system (1) can be obtained,

$A_l{x}_l=b_l.$

We have converted it equivalently to

$\mathcal{X}{\times }_1{A}_l^{\left(1\right)}+\mathcal{X}{\times }_2{A}_l^{\left(2\right)}+\mathcal{X}{\times }_3{A}_l^{\left(3\right)}={\mathcal{D}}_l,$ (3)

where, $A_l^{\left(i\right)}\in {ℝ}^{J_k\times J_k}$ ($k=1,2,3$ ) and ${\mathcal{D}}_l$ is the right-hand term, the structure of coefficient matrix $A_l^{\left(i\right)}$ is as follows

$A_l^{\left(i\right)}={\left(M_l^{\left(i\right)}\right)}^{-1}{\tilde{A}}_l^{\left(i\right)},$

Here, $M_l^{\left(i\right)}$ represents the mass matrix, and ${\tilde{A}}_l^{\left(i\right)}$ is the coefficient matrix obtained by the finite element method.

2.2. Inverse of Mass Matrix $M$

According to [7], we can use the Stenger quadrature formula to approximate the inverse of the mass matrix $M$ .

Let $M\in {ℝ}^{n\times n}$ be a symmetric positive definite matrix, with all eigenvalues $\lambda >0$ . Consider the integral

${\displaystyle {\int }_0^{\infty }{\text{e}}^{tM}\text{d}t}.$

For each eigenvalue $\lambda$ , the integral is

${\displaystyle {\int }_0^{\infty }{\text{e}}^{-\lambda t}\text{d}t}.$

Due to the relationship between the matrix exponential and the eigenvalues, we have

${\displaystyle {\int }_0^{\infty }{\text{e}}^{-tM}\text{d}t}=M^{-1}.$

Let $p\in \mathbb{N}$ , Step length is $h_{\text{st}}={\pi }^2/\sqrt{p}$ , the integration nodes and weights are as follows

$t_j=\text{ln}\left({\text{e}}^{j{h}_{\text{st}}}+\sqrt{1+{\text{e}}^{2j{h}_{\text{st}}}}\right),j=-p,\cdots ,p.$

$w_j=\frac{h_{\text{st}}}{\sqrt{1+{\text{e}}^{-2j{h}_{\text{st}}}}}.$

Then the integral is approximately equal to

${\displaystyle {\int }_0^{\infty }{\text{e}}^{-t\lambda }\text{d}t}\approx {\displaystyle \sum }_{j=-p}^p{w}_j{\text{e}}^{-t_j\lambda }.$

2.3. ECTMG Algorithm

Combining Sections 2.1 and 2.2, we present the following iterative method (Algorithm 1).

Algorithm 1. The calculation of $A^{\left(i\right)}$

Similar to [7], we employ the following algorithm for the solution (Algorithm 2).

Algorithm 2. ECTMG algorithm [7].

With is prologation [8], $T_l^{m_l}$ is the smoother. And $m_l$ represents the number of iterations, we similar to [9] provides the following selection rules:

When $d=2$

1) If $l>L_0$ , then $m_l=\left[m_L{\beta }^{L-l}\right]$ .

2) If $l\le L_0$ , $m_l=\left[m_{\text{*}}^{\frac{1}{2}}\left(L-\left(2-{\epsilon }_0\right)l\right){\kappa }_l\right]$ .

When $d\ge 3$

1) If $l>L_0$ , then $m_l=\left[m_L{\beta }^{L-l}\right].$

2) If $l\le L_0$ , there are two cases:

a) If $\left(2-{\epsilon }_0\right)L_0\le L$ , then $m_l=\left[m_{\text{*}}^{\frac{1}{2}}\left(L-\left(2-{\epsilon }_0\right)l\right){\kappa }_l\right]$ .

b) If $\left(2-{\epsilon }_0\right)L_0>L$ , there exists a positive integer$L^{\prime }<L_0$ such that $\left(2-{\epsilon }_0\right)L^{\prime }\le L$ , for all $l\le L^{\prime }$ , choose $m_l=\left[m_{\text{*}}\left(L-\left(2-{\epsilon }_0\right)l\right){\kappa }_l\right]$ .

Remark Let ${\kappa }_l=h_l^{-2}$ , and ${\epsilon }_0$ be a positive constant in the interval [0, 1].

3. Numerical Example

In this section, we verify the effectiveness of the ECTMG algorithm through some numerical examples. All tests will be done with configuration: 11th Gen Intel(R) Core(TM) i5-11400H @ 2.70 GHz 2.69 GHz. Let CPU(s) represent the iteration time. ${\epsilon }_0=0.1$ .

We define the error between the exact solution ${\mathcal{X}}^{\text{*}}$ and the numerical solution ${\mathcal{X}}^k$ as follows

$E\left({\mathcal{X}}^{*}-{\mathcal{X}}^k\right)=\underset{i_1,i_2,\cdots ,i_N}{\max }\left|x_{i_1,i_2,\cdots ,i_N}^{*}-x_{i_1,i_2,\cdots ,i_N}^k\right|,$

Obviously, we have $E\left({\mathcal{X}}^{\text{*}}-{\mathcal{X}}^k\right)={\Vert vec\left({\mathcal{X}}^{\text{*}}-{\mathcal{X}}^k\right)\Vert }_{\infty }$ . In the table of numerical examples in the following text, the symbol $†$ indicates insufficient memory during the algorithm’s execution.

Example 1. ([8]) Consider the following Poisson equation

$\begin{array}{lll}-△u=f\hfill & \text{in}\hfill & V={\left[0,1\right]}^3\hfill \\ u=0\hfill & \text{on}\hfill & \partial V.\hfill \end{array}$

Taking the partition step length as $h=\frac{1}{n+1}$ , using the hexahedral element, we can obtain the Sylvester tensor Equation (1), with $A^{\left(k\right)}$ in Equation (2) as follows

$A^{\left(k\right)}={\left(M^{\left(k\right)}\right)}^{-1}{\tilde{A}}^{\left(k\right)},$

with ${\tilde{A}}^{\left(k\right)}$ as follows

${\tilde{A}}^{\left(k\right)}=\frac{2}{9}\left[\begin{array}{ccc}4& -2& 0\\ -2& 4& -2\\ 0& -2& 4\end{array}\right].$

and $M^{\left(k\right)}$ as follows

$M^{\left(k\right)}=\frac{1}{4}\left[\begin{array}{ccc}4& 1& 0\\ 1& 4& 1\\ 0& 1& 4\end{array}\right].$

Set $V_l=2/{\displaystyle \sum }_{k=1}^3\left({\alpha }_l^k{\beta }_l^k\right)$ , where ${\alpha }_l^k$ and ${\beta }_l^k$ are the maximum and minimum eigenvalues of matrix $A_l^{\left(k\right)}$ , respectively. Select the right-hand side tensor $ℬ$ such that the solution to the corresponding 3-order Poisson equation is

$u=\frac{1}{3}\cdot {\displaystyle \prod }_{k=1}^3\left(x_k-x_k^2\right)$ , where ${\left(x_k\right)}_i=i/\left(n+1\right)$ , $i=1,2,\cdots ,n$ . Table 1 and Table 2

lists the numerical results of BiCG_BTF and ECTMG(BiCG_BTF) algorithms.

Table 1. Numerical results for $n=255\times 255\times 255$ .

Algorithm	CPU (s)	$\mathbb{E}\left({\mathcal{X}}^{\text{*}}-{\mathcal{X}}_L\right)$	IT
BiCG_BTF	663.32	2.20 × 10−7	994
ECTMG (BiCG_BTF)	142.81	4.72 × 10−7	(3, 205, 256, 1024, 4096, 512, 64)

Remark: In the BiCG_BTF algorithm, IT represents the total number of iterations required for the entire algorithm process; in the ECTMG algorithm, it is the number of iterations $m_l$ for the corresponding level.

Table 2. Numerical results for $n=511\times 511\times 511$ .

Algorithm	CPU (s)	$\mathbb{E}\left({\mathcal{X}}^{\text{*}}-{\mathcal{X}}_L\right)$	IT
BiCG_BTF	14076.29	4.96 × 10−7	1958
ECTMG (BiCG_BTF)	790.07	5.78 × 10−7	(3, 269, 256, 1024, 15,360, 1920, 240, 30)

As can be seen from Table 1 and Table 2, through the linear system (2) discretized by the finite element method, the ECTMG algorithm can still efficiently solve the problem and maintain the same accuracy as BiCG_BTF, which verifies the feasibility of the finite element discretized equation for this solution framework.

4. Conclusion

Based on the finite element method, we discretize the three-dimensional partial differential equations and solves it using the economic cascadic tensor multigrid method. The numerical results show that the new method is effective.

Acknowledgements

We acknowledge the financial support from the Natural Science Foundation of China under Grant 12161027 and the Science and Technology Project of Guangxi (Guike AD25069086).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Zhao, Y., Chang, J.C. and Liang, J.L. (2024) Simulation of Temperature Field of Rotary Kiln Burner Based on Finite Element Method. Applied Mathematics Journal of Chinese Universities, 39, 51-63. (In Chinese)
[2]	Wu, H.X., Liu, X.Y., Wang, T.Y., et al. (2026) Finite Element Simulation of Scoliosis with Muscle Unit Introduction: Verification of Correction Effect under Bidirectional Load. Chinese Journal of Tissue Engineering Research, 30, 2172-2181. (In Chinese)
[3]	Gao, X. and Xing, W.H. (2023) Application of Finite Element Analysis in Spine Surgery. Chinese Journal of Tissue Engineering Research, 27, 2921-2927. (In Chinese)
[4]	Nie, Q.H., Song, X.D., Liu, S., et al. (2025) Advances in Application of Finite Element Numerical Simulations to Electrocatalytic Reduction Systems of Energy-Related Small Molecules. Clean Coal Technology, 31, 1-13. (In Chinese)
[5]	Grasedyck, L. (2004) Existence and Computation of Low Kronecker-Rank Approximations for Large Linear Systems of Tensor Product Structure. Computing, 72, 247-265.[CrossRef]
[6]	Chen, Z. and Lu, L. (2012) A Projection Method and Kronecker Product Preconditioner for Solving Sylvester Tensor Equations. Science China Mathematics, 55, 1281-1292.[CrossRef]
[7]	Huang, J.Y. and Li, C.L. (2026) An Economic Cascadic Tensor Multigrid Method for Solving High Dimensional Elliptic Linear Partial Differential Problems. arXiv preprint arXiv: 2606. 25643.
[8]	Ballani, J. and Grasedyck, L. (2013) A Projection Method to Solve Linear Systems in Tensor Format. Numerical Linear Algebra with Applications, 20, 27-43.[CrossRef]
[9]	Shi, Z., Xu, X. and Huang, Y. (2007) Economical cascadic multigrid method (ECMG). Science in China Series A: Mathematics, 50, 1765-1780. [Google Scholar] [CrossRef]