An efficient multigrid solver for 3D biharmonic equation with a discretization by 25-point difference scheme

TL;DR

This paper introduces an efficient multigrid solver combining extrapolation techniques and a 25-point difference scheme to solve the 3D biharmonic equation with high accuracy and computational efficiency.

Contribution

The paper develops a novel extrapolation cascadic multigrid method with 25-point discretization and demonstrates its effectiveness for 3D biharmonic problems.

Findings

The method achieves higher-order accuracy on the finest grid.

Numerical experiments confirm the efficiency and accuracy of the solver.

The approach outperforms traditional methods in solving 3D biharmonic equations.

Abstract

In this paper, we propose an efficient extrapolation cascadic multigrid (EXCMG) method combined with 25-point difference approximation to solve the three-dimensional biharmonic equation. First, through applying Richardson extrapolation and quadratic interpolation on numerical solutions on current and previous grids, a third-order approximation to the finite difference solution can be obtained and used as the iterative initial guess on the next finer grid. Then we adopt the bi-conjugate gradient (Bi-CG) method to solve the large linear system resulting from the 25-point difference approximation. In addition, an extrapolation method based on midpoint extrapolation formula is used to achieve higher-order accuracy on the entire finest grid. Finally, some numerical experiments are performed to show that the EXCMG method is an efficient solver for the 3D biharmonic equation.

Tables

Table 1. Table 1: Errors and convergence rates using EXCMG bi-cg for Problem 1.

Mesh	Iters	${\Vert u_h-u\Vert }_2$		${\Vert u_h-u\Vert }_{\infty }$		${\Vert w_h-u_h\Vert }_2$
		Error	Order	Error	Order	Error	Order
$32\times 32\times 32$	474	$1.13\left(-2\right)$		$5.16\left(-2\right)$		$5.14\left(-3\right)$
$64\times 64\times 64$	512	$2.89\left(-3\right)$	1.97	$1.29\left(-2\right)$	2.00	$6.15\left(-4\right)$	3.06
$128\times 128\times 128$	64	$7.27\left(-4\right)$	1.99	$3.21\left(-3\right)$	2.00	$7.63\left(-5\right)$	3.01
$256\times 256\times 256$	8	$1.80\left(-4\right)$	2.01	$8.01\left(-4\right)$	2.00	$9.56\left(-6\right)$	3.00
$512\times 512\times 512$	1	$4.25\left(-5\right)$	2.08	$1.99\left(-4\right)$	2.01	$1.19\left(-6\right)$	3.00
	4.12 WU			$1.46\left(-5\right)$		$5.81\left(-6\right)$

Table 2. Table 2: Errors and convergence rates using EXCMG bi-cg for Problem 2.

Mesh	Iters	${\Vert u_h-u\Vert }_2$		${\Vert u_h-u\Vert }_{\infty }$		${\Vert w_h-u_h\Vert }_2$
		Error	Order	Error	Order	Error	Order
$32\times 32\times 32$	259	$8.96\left(-7\right)$		$8.06\left(-6\right)$		$4.59\left(-6\right)$
$64\times 64\times 64$	470	$2.30\left(-7\right)$	1.96	$2.06\left(-6\right)$	1.97	$5.48\left(-7\right)$	3.07
$128\times 128\times 128$	384	$5.80\left(-8\right)$	1.99	$5.15\left(-7\right)$	2.00	$6.66\left(-8\right)$	3.04
$256\times 256\times 256$	48	$1.46\left(-8\right)$	1.99	$1.28\left(-7\right)$	2.01	$8.19\left(-9\right)$	3.02
$512\times 512\times 512$	6	$3.67\left(-9\right)$	1.99	$3.22\left(-8\right)$	1.99	$1.00\left(-9\right)$	3.03
	18.98 WU			$9.90\left(-9\right)$		$2.93\left(-10\right)$

Table 3. Table 3: Errors and convergence rates using EXCMG bi-cg for Problem 3.

Mesh	Iters	${\Vert u_h-u\Vert }_2$		${\Vert u_h-u\Vert }_{\infty }$		${\Vert w_h-u_h\Vert }_2$
		Error	Order	Error	Order	Error	Order
$32\times 32\times 32$	285	$4.10\left(-6\right)$		$1.75\left(-5\right)$		$9.49\left(-6\right)$
$64\times 64\times 64$	533	$1.05\left(-6\right)$	1.96	$4.36\left(-6\right)$	2.00	$1.15\left(-6\right)$	3.04
$128\times 128\times 128$	384	$2.66\left(-7\right)$	1.98	$1.09\left(-6\right)$	2.00	$1.42\left(-7\right)$	3.02
$256\times 256\times 256$	48	$6.71\left(-8\right)$	1.99	$2.73\left(-7\right)$	2.00	$1.77\left(-8\right)$	3.01
$512\times 512\times 512$	6	$1.71\left(-8\right)$	1.97	$6.90\left(-8\right)$	1.98	$2.18\left(-9\right)$	3.02
	19.11 WU			$8.68\left(-9\right)$		$9.16\left(-10\right)$

Table 4. Table 4: Errors and convergence rates using EXCMG bi-cg for Problem 4.

Mesh	Iters	${\Vert u_h-u\Vert }_2$		${\Vert u_h-u\Vert }_{\infty }$		${\Vert w_h-u_h\Vert }_2$
		Error	Order	Error	Order	Error	Order
$32\times 32\times 32$	275	$1.35\left(-6\right)$		$3.47\left(-6\right)$		$2.19\left(-6\right)$
$64\times 64\times 64$	513	$3.47\left(-7\right)$	1.96	$8.69\left(-7\right)$	2.00	$2.68\left(-7\right)$	3.03
$128\times 128\times 128$	512	$8.77\left(-8\right)$	1.98	$2.17\left(-7\right)$	2.00	$3.30\left(-8\right)$	3.02
$256\times 256\times 256$	64	$2.22\left(-8\right)$	1.98	$5.46\left(-8\right)$	1.99	$4.11\left(-9\right)$	3.01
$512\times 512\times 512$	8	$5.70\left(-9\right)$	1.96	$1.39\left(-8\right)$	1.98	$5.11\left(-10\right)$	3.01
	25.07 WU			$1.50\left(-9\right)$		$2.72\left(-10\right)$

Table 5. Table 5: Errors and convergence rates using EXCMG bi-cg for Problem 5.

Mesh	Iters	${\Vert u_h-u\Vert }_2$		${\Vert u_h-u\Vert }_{\infty }$		${\Vert w_h-u_h\Vert }_2$
		Error	Order	Error	Order	Error	Order
$32\times 32\times 32$	432	$8.86\left(-2\right)$		$3.76\left(-1\right)$		$1.91\left(-1\right)$
$64\times 64\times 64$	873	$2.42\left(-2\right)$	1.87	$1.01\left(-1\right)$	1.90	$2.55\left(-2\right)$	2.91
$128\times 128\times 128$	1913	$6.22\left(-3\right)$	1.96	$2.57\left(-2\right)$	1.97	$2.33\left(-3\right)$	3.45
$256\times 256\times 256$	256	$1.60\left(-3\right)$	1.96	$6.45\left(-3\right)$	1.99	$2.16\left(-4\right)$	3.43
$512\times 512\times 512$	32	$4.45\left(-4\right)$	1.85	$1.65\left(-3\right)$	1.97	$2.45\left(-5\right)$	3.14
	95.70 WU			$2.84\left(-4\right)$		$7.89\left(-5\right)$

Equations

$${\Delta }^{2}u(x,y,z)=f(x,y,z),(x,y,z)\in \Omega =(0,1{)}^3,$$

$$u(x,y,z)=g_1(x,y,z),\frac{\partial u}{\partial n}=g_2(x,y,z),(x,y,z)\in \partial \Omega ,$$

$$u(x,y,z)=g_1(x,y,z),\frac{{\partial }^{2}u}{\partial n^2}=g_2(x,y,z),(x,y,z)\in \partial \Omega .$$

$${\Delta }^{2}u(x,y,z)=\frac{{\partial }^{4}u}{\partial x^4}+\frac{{\partial }^{4}u}{\partial y^4}+\frac{{\partial }^{4}u}{\partial z^4}+2\frac{{\partial }^{4}u}{\partial x^2{y}^2}+2\frac{{\partial }^{4}u}{\partial x^2{z}^2}+2\frac{{\partial }^{4}u}{\partial y^2{z}^2}.$$

$$\frac{{\partial }^{4}u}{\partial x^4}+\frac{{\partial }^{4}u}{\partial y^4}+2\frac{{\partial }^{4}u}{\partial x^2{y}^2}=f(x,y).$$

$$42u_{i,j,k}-12(u_{i-1,j,k}+u_{i+1,j,k}+u_{i,j-1,k}+u_{i,j+1,k}+u_{i,j,k-1}+u_{i,j,k+1})$$

$$+u_{i-2,j,k}+u_{i+2,j,k}+u_{i,j-2,k}+u_{i,j+2,k}+u_{i,j,k-2}+u_{i,j,k+2}+$$

$$+2(u_{i-1,j-1,k}+u_{i-1,j+1,k}+u_{i+1,j-1,k}+u_{i+1,j+1,k}+u_{i-1,j,k-1}+u_{i+1,j,k-1}+u_{i,j-1,k-1}+u_{i,j+1,k-1}$$

$$+u_{i-1,j,k+1}+u_{i+1,j,k+1}+u_{i,j-1,k+1}+u_{i,j+1,k+1})=h^4{f}_{i,j,k}.$$

$$\frac{u_{1,j,k}-u_{-1,j,k}}{2h}={\left(\frac{\partial u}{\partial x}\right)}_{0,j,k},$$

$$u_{-1,j,k}=u_{1,j,k}-2h{\left(\frac{\partial u}{\partial x}\right)}_{0,j,k}.$$

$$\frac{u_{1,j,k}-2u_{0,j,k}+u_{-1,j,k}}{h^2}={\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_{0,j,k},$$

$$u_{-1,j,k}=-u_{1,j,k}+2u_{0,j,k}+h^2{\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_{0,j,k}.$$

$$A_h{u}_h=f_h,$$

$$e^i(x_k)=A(x_k)h_i^2+O(h_i^4),$$

$${\tilde{u}}_k^1:=\frac{4u_k^1-u_k^0}{3}=u(x_k)+O(h_0^4),k=j,j+1.$$

$${\tilde{u}}_{j+1/2}^1:=u_{j+1/2}^1+\frac{1}{6}(u_j^1-u_j^0+u_{j+1}^1-u_{j+1}^0)=u(x_{j+1/2})+O(h_0^4),$$

$$A(x_k)=\frac{4}{3h_0^2}(u_k^0-u_k^1)+O(h_0^2),k=j,j+1.$$

$$A(x_{j+1/2})=\frac{1}{2}(A(x_j)+A(x_{j+1}))+O(h_0^2).$$

$$A(x_{j+1/2})=\frac{2}{3h_0^2}(u_j^0-u_j^1)+\frac{2}{3h_0^2}(u_{j+1}^0-u_{j+1}^1)+O(h_0^2).$$

$$u_{j+1/2}^1=u(x_{j+1/2})+\frac{1}{4}A(x_{j+1/2})h_0^2+O(h_0^4),$$

$$\Big{\{}x_{j},x_{j+1/4},x_{j+1/2},x_{j+3/4},x_{j+1}\Big{\}}$$

$$c{u}^1+(1-c)u^0=u^2+O(h_0^4),$$

$$w_k^2:=\frac{5u_k^1-u_k^0}{4}=u_k^2+O(h_0^4),k=j,j+1,$$

$$u_{j+1/2}^2=u_{j+1/2}^1-\frac{3}{16}A(x_{j+1/2})h_0^2+O(h_0^4),$$

$$w_{j+1/2}^2:=u_{j+1/2}^1+\frac{1}{8}(u_j^1-u_j^0+u_{j+1}^1-u_{j+1}^0)=u_{j+1/2}^2+O(h_0^4).$$

$$w_{j+1/4}^2$$

$$w_{j+3/4}^2$$

$$w_{j+1/4}^2=u_{j+1/4}^2+O(h_0^3),$$

$$w_{j+3/4}^2=u_{j+3/4}^2+O(h_0^3),$$

$$w^2(\xi ,\eta ,\zeta )=\sum _i{N}_i(\xi ,\eta ,\zeta )w_i^2,$$

$$N_i(\xi ,\eta ,\zeta )=l_i^2(\xi )l_i^2(\eta )l_i^2(\zeta )$$

Full text

1 Introduction

In this paper, we consider the following three-dimensional (3D) biharmonic equation

[TABLE]

with Dirichlet boundary conditions of first kind

[TABLE]

or Dirichlet boundary conditions of second kind

[TABLE]

The biharmonic operator $\Delta^{2}$ in three-dimensional (3D) Cartesian coordinates can be written as

[TABLE]

And the two dimensional (2D) version of Eq. (1) is

[TABLE]

The biharmonic equation is a fourth-order partial differential equation which arises in areas of continuum mechanics, including linear elasticity theory, phase-field models and Stokes flows. Due to the significance of the biharmonic equation, a large number of numerical methods for solving the biharmonic equations have been proposed [6, 7, 1, 2, 3, 4, 5, 14, 8, 9, 11, 15, 13, 12, 10, 16]. Most of these works focus on two-dimensional case. There has been very little work devoted to solving the 3D biharmonic equations. The main reason is that 3D problems require large computational power and memory storage [14, 8].

Various methods for the numerical solutions of the biharmonic equations have been considered in the literature. A popular technique is to split $\Delta^{2}u=f$ into two coupled Poisson equations for $u$ and $v$: $\Delta u=v,\Delta v=f$, each equation can be solved by using fast Poisson solvers. The coupled method has been widely used by many authors [2, 6, 7]. As it is mentioned in [2, 6, 7], the main difficulty for the coupled (splitting or mixed) method is that the boundary conditions for the newly introduced variable $v$ are undefined and needs to be approximated accurately, and the computational results strongly depends on the choice of the approximation of missing boundary values for $v$.

Another conventional approach for solving the 3D biharmonic equations is to directly discretize Eq. (1) on a uniform grid using a 25-point computational stencil with truncation error of order $h^{2}$, which is derived by Ribeiro Dos Santos [5] in 1967. This conventional 25-point difference approximation connects the value of $u$ at grid $(x_{i},y_{j},z_{k})$ in terms of 24 neighboring values in a $5\times 5\times 5$ cube. Thus, this direct method need to be modified at grid points near the boundaries. As mentioned in [1, 13, 14], there are serious computational difficulties with solution of the linear systems obtained by the 13-point discretization of the 2D biharmonic equation and the 25-point discretization of 3D biharmonic equation. Dehghan and Mohebbi [8] also pointed that this direct method can only be used for moderate values of grid width $h$ and the well-known iterative methods such Jacobi or Gauss-Seidel either converge very slowly or diverge.

The combined compact difference method is another popular method for solving the biharmonic equation [8, 14]. For example, Altas et al. [14] proposed a fourth-order, combined compact formulation, where The unknown solution and its first derivatives are carried as unkonws at grid point and computed simultaneously, for the 3D biharmonci equation with Dirichlet boundary conditions of first kind. In 2006, Dehghan et al. [8] proposed two combined compact difference schemes for solve 3D biharmonic equation with Dirichlet boundary conditions of second kind, which use the known solution and its second derivatives as unknowns. In these combined compact difference methods, there is no need to modify the difference scheme at grid points near the boundaries, and the given Dirichlet boundary conditions are exactly satisfied and no approximations need to be carried out at the boundaries, in contrary to the coupled method. However, these combined compact difference methods introduce extra amount of computation, and the classical iterations for solving the resulting linear system suffer from slow convergence. Multigrid methods give good results in [8] and [14]. However, numerical results in [8] and [14] are reported only up to $32\times 32\times 32$ and $64\times 64\times 64$ grids, respectively. To the best of our knowledge, there is no numerical results for solving the 3D biharmonic equations with large-scale discretized meshes.

In this paper, we propose an efficient extrapolation cascadic multigrid method based on the conventional 25-point approximation to solve 3D biharmonic equations with both first and second boundary conditions. In our method, the conventional 25-point difference scheme is used to approximate the 3D biharmonic equation (1). In order to overcome the serious computational difficulties with solution of the resulting linear system, by combining Richardson extrapolation and quadratic interpolation on numerical solutions on current and previous grids, we obtain quite good initial guess of the iterative solution on the next finer grid, and then adopt the bi-conjugate gradient (Bi-CG) method to solve the large linear system efficiently. Our method has been used to solve 3D biharmonic problems with more than 135 million unknowns with only several iterations.

The rest of the paper is organized as follows: Section 2 presents the 25-point difference approximation for the 3D biharmonic equation and its modification of the difference scheme at grid points near boundaries. Section 3 reviews the classical V-cycle and W-cycle multigrid methods. In Section 4, we present a new EXCMG method to solve the linear three-dimensional biharmonic equation (1). Section 5 describes the Bi-CG solver in our new EXCMG method. Section 6 provides the numerical results to demonstrate the high efficiency and accuracy of the proposed method, and conclusions are given in the final section.

2 Second-order Finite Difference Discretization

We consider a cubic domain $\Omega=[0,1]\times[0,1]\times[0,1]$. Let $N=1/h$ be the numbers of uniform intervals along all the $x$, $y$ and $z$ directions. We discretize the domain with unequal meshsizes $h=1/N$ in all $x,y$ and $z$ coordinate directions. The grid points are ($x_{i},y_{j},z_{k}$), with $x_{i}=ih,y_{j}=jh$ and $z_{k}=kh,i,j,k=0,1,\cdots,N$. The quantity $u_{i,j,k}$ represents the numerical solution at ($x_{i},y_{j},z_{k}$).

Then the value on the boundary points $u_{i,j,k}$ can be evaluated directly from the Dirichlet boundary condition. For internal grid points ($i=2,\cdots,N-2,j=2,\cdots,N-2,k=2,\cdots,N-2$), the 25-point second-order difference scheme for 3D biharmonic equation was derived [5, 14]:

[TABLE]

Note that $u_{i,j,k}$ is connected to grid points two grids away in each direction from the point $(x_{i},y_{j},z_{k})$. Thus, the above difference formulation (2) for the grid points near the domain boundary $\partial\Omega$ involves at least one value of point outside the domain, and these points outside the domain are fictitious points which need to be replaced by the internal points through the boundary condition. These could be done for both first and second kind of boundary conditions.

For the first kind of boundary condition. For example, for $i=1$, the point ($x_{i-2},y_{j},z_{k}$) lies outside the computational domain, and the value on the fictitious point ($x_{-1},y_{j},z_{k}$) can be obtained through the following central difference formula called the reflection formulas [3, 4]:

[TABLE]

where $\left(\frac{\partial u}{\partial x}\right)_{0,j,k}$ can be obtained from the boundary condition (2) and $u_{-1,j,k}$ is given by

[TABLE]

For the second kind of boundary condition. For example, for $i=1$, the point ($x_{i-2},y_{j},z_{k}$) lies outside the computational domain, and the value on the fictitious point ($x_{-1},y_{j},z_{k}$) can also be obtained through the following central difference formula called the reflection formulas:

[TABLE]

where $u_{0,j,k}$ and $\left(\frac{\partial^{2}u}{\partial x^{2}}\right)_{0,j,k}$ can be obtained from the boundary condition (2) and $u_{-1,j,k}$ is given by

[TABLE]

We use $u_{h}$ and $u_{h/2}$ to represent the finite difference solutions of equation (1) with mesh sizes $h$ and $\frac{h}{2}$ respectively. Afterward, a matrix form, which express the finite difference scheme (2) and an equation set including formulas of the grid points near the boundary, can be obtained as below:

[TABLE]

Where $A_{h}$ is not a symmetry positive definite matrix, and the right hand-side vector of (2) and an equation set including the formulas of the grid points near the boundary are expressed by $f_{h}$.

Note that the discretization equations for grid points that away from the boundary and that near the boundary are different, one must distinguish all possible cases. Although there are a little bit troublesome to treat all cases (there are totally 27 cases with 27 different equations), by moving the known boundary values into the right hand-side of the system, it is convenient to solve these equations which only involves unknown on the grid points.

3 Classical Multigrid Method

Since the 1970s’, many scholars have done researches on the classical multigrid method. Through deep researches on it for about fifty years, the classical multigrid method gradually forms its own comprehensive system. Including the interpolation, restriction and iteration, the classical multigrid method starts from the fine grid, goes to coarse grid and then returns to the fine grid. The classical multigrid methods contain V-cycle and W-cycle.

The classical multigrid method is introduced in detail with several steps. First, the specific smoother is used to smooth the current approximation on the fine grid. To obtain more oscillatory error components, we compute the residual and transfer it to the coarser grid with restriction. Next, we solve the residual equation on the coarser grid with the application of the number ($\gamma$) of cycles. From the fine grid to the coarsest grid and back to the fine grid is called a cycle. Then, we acquire the improved approximation on the fine grid by interpolating the correction back to the fine grid. Finally, we smooth the obtained approximation on the fine grid with the smoother again. If $\gamma$=1, call it V-cycle. And if $\gamma$=2, call it W-cycle. We take the four-level structures of V-cycle and W-cycle in Fig.1 for instances to illustrate that.

Remark 1

When the $\gamma$-cycle is performed on the coarsest grid, direct solver is used to solve the residual equation.

4 Extrapolation Cascadic Multigrid Methods

It is an important issue to find approaches to solving the linear equation with enormous unknowns, which is obtained by FE and FD discretizations. Therefore, many authors paid great attention on it and presented multigrid methods including the MG method, the CMG method and the EXCMG method. The MG method has had a nearly integrated system through many scholars’ hard work in the past several decades. However, its algorithm is complex. Then the CMG method proposed by Deuflhard and Bornemann in [17] only use the interpolation and iteration so that its algorithm which is easy to operate is appealing. Furthermore, in 2008, the EXCMG method was proposed by Chen et al. [26] and the cores of it are Richardson extrapolation and quadric interpolation. Compared with the CMG method, the EXCMG method provides a much better initial guess for the iteration solution on the next finer grid. In this section, we propose a new EXCMG method combined with the second-order incompact FD discretization for solving the linear three-dimensional biharmonic equation.

4.1 Description of the EXCMG Algorithm

In Algorithm 1, H is the size of the coarsest grid. L, the positive integer, denotes the total number of grids except first two embedded grids and indicates that the finest grids’ size is $\frac{H}{2^{L+1}}$. For the sizes of first two coarse grids are small, DSOLVE, a direct solver, is applied on the first two coarse grids (see line 1-2 in the Algorithm 1). In addition, procedure $\textrm{EXP}_{finite}(u_{2h},u_{4h})$ represents the third-order approximation of the FD solution $u_{h}$ which is obtained by Richardson extrapolation and quadratic interpolation from numerical solutions $u_{2h}$ and $u_{4h}$. Meanwhile, a selective step is presented in the Algorithm 1 above where $\textrm{EXP}_{true}(u_{h},u_{2h})$ refers to a higher-order solution which is extrapolated on the finest grid with the mesh size, $h$, from two second-order numerical solutions $u_{h}$ and $u_{2h}$.

The details of the procedure of extrapolation and quadratic interpolation are introduced next subsection 4.2. The difference between our new EXCMG method and existing EXCMG method are discussed below:

Instead of applying the second-order linear FE method, a second-order incompact difference scheme is used to discretize the 3D biharmonic equation in our new EXCMG method. 2. 2.

Rather than perform the fixed number of iterations used in the existing EXCMG method, we introduce a relative residual tolerance $\epsilon_{i}$ into the Bi-CG solver (see line 7 in the Algorithm 1), which enables us to avoid the difficulty of determining the number of iterations at every grid level and obtain numerical solutions with desired accuracy conveniently. 3. 3.

In our new EXCMG method, we take the Bi-CG solver as smoother instead of the CG solver (see line 8 in Algorithm 1). The Bi-CG is more suitable for positive definite matrix which is not symmetric compared with the CG solver. 4. 4.

Through $\textrm{EXP}_{true}(u_{h},u_{2h})$, a higher-order extrapolated solution $\tilde{u}_{h}$ is obtained easily, which improves the accuracy of the numerical solution $u_{h}$ (see line 11 in Algorithm 1).

4.2 Extrapolation and Quadratic Interpolation

The Richardson extrapolation is a well-known method for producing more accurate solutions of many problems in numerical analysis. Marchuk and Shaidurov [21] researched the application of the Richardson extrapolation on the FD method systematically in 1983. Since then, this technique has been well demonstrated in the frame of the FE and FD methods [18, 21, 20, 22, 23, 24, 25, 33, 34, 35, 36].

In next three subsections, we will give the explanation for how to obtain higher-oder accuracy solution on the fine grid. Moreover, how to acquire a third-order approximation of the second-order FD method on the next finer grid is illustrated as well. Meanwhile, we can regard it as another critical application of the extrapolation method which produces good initial guesses for iterative solutions.

4.2.1 Extrapolation for the True Solution

For simplicity, we first consider the three-level of embedded grids $Z_{i}$(i=0, 1, 2) with mesh sizes $h_{i}$ = $h_{0}$/$2^{i}$ in one dimension. In addition, let $e^{i}$ = $u^{i}$-u be the error of the second-order incompact FD solution $u^{i}$ with mesh size $h_{i}$. We make an assumption that the error at the node has the following form:

[TABLE]

where A(x) is a properly smooth function independent of $h_{i}$. We will verify the error expansion (12) by numerical results in Sect. 5.

Through the equation (12), the Richardson extrapolation formula at the coarse grid point is obtained

[TABLE]

Then, a midpoint extrapolation formula is obtained by linear interpolation

[TABLE]

whose accuracy is fourth-order at fine grid points.

From equation (12), it is easy to obtain

[TABLE]

Through the error estimate of the linear interpolation

[TABLE]

and substituting equation (15) to equation (16), yield

[TABLE]

Since

[TABLE]

The midpoint extrapolation formula (14) is obtained by the use of equation (17).

4.2.2 Extrapolation for the FD Solution

In this subsection, given solutions $u^{0}$ and $u^{1}$ of the second-order FD method, we will explain how to construct a third-order approximation $w^{2}$ of the FD solution $u^{2}$ by using extrapolation and interpolation methods in detail.

We divide the coarse element ($x_{j}$,$x_{j+1}$) into four uniform elements by adding one midpoint and two four equal points which are on the left side and right side of the midpoint. As a result, a set which contains five points and belongs to fine mesh $Z_{2}$ is obtained

[TABLE]

To acquire the more accurate approximation of FD solution $u^{2}$, the given solutions $u^{0}$ and $u^{1}$ are combined linearly. Therefore, here assume the existence of a constant c such that

[TABLE]

For obtaining the value of the constant c, we substitute the error expansion (12) into (19) and obtain c=5/4. Afterward, we obtain formulas of node extrapolation at points $x_{j}$ and $x_{j+1}$.

[TABLE]

Next, derivate the midpoint $x_{j+1/2}$’s extrapolation formula. First, use the error expansion (12) again and obtain the formula below

[TABLE]

Then substitute the (17) to (21) for eliminating the unknown A($x_{j+1/2}$) and the following extrapolation formula of $x_{j+1/2}$ is yielded

[TABLE]

Finally, since the values of three points $w_{j}^{2}$, $w_{j+1/2}^{0}$ and $w_{j+1}^{0}$ have been obtained, we can derivate the extrapolation formulas at the points $x_{j+1/4}^{2}$ and $x_{j+3/4}^{2}$ blow through the use of the quadratic interpolation method

[TABLE]

From the polynomial interpolation’s theory, it is easy to demonstrate that the third-order approximation of the FD solution can be presented by formulas (25) and (26). i.e.,

[TABLE]

[TABLE]

4.2.3 Application of Extrapolation and Quartic Interpolation on Three-Dimension

In this subsection, how to acquire an accurate third-order approximation $w^{2}$ of the FD solution $u^{2}$ will be explained for embedded cubic mesh as shown in Fig.3. The specific steps of the construction of the third-order approximation $u^{2}$ are illustrated below:

Step 1: Corner nodes (1, 5, 21, 25, 101, 105, 121, 125): Use the extrapolation formula (20) to obtain the approximations of the values of 8 corner nodes on interpolation cells.

Step 2: Midpoints of edges (3, 11, 15, 23, 51, 55, 71, 75, 103, 111, 115, 123): Use the midpoint extrapolation formula (22) to obtain the approximations of the values of 19 midpoints of edges on interpolation cells.

Step 3: Centers of faces (13, 53, 65, 61, 77, 113): View the center of each face as the midpoint of the two face diagonals on interpolation cells. To obtain the approximations of the values of them, use the midpoint extrapolation formula (22) obtaining two approximations, calculate the arithmetic mean of the two obtained approximations and treat it as the approximate value at the center of each face.

Step 4: Center of the hexahedral element (63): View the center of the hexahedral element as the midpoint of the four space diagonals on interpolation cells. To obtain the approximation of the value of it, use the midpoint extrapolation formula (22) obtaining four approximations, calculate the arithmetic mean of the four obtained approximations and treat it as the approximate value at the hexahedral element.

Step 5: Other 98 fine grid points: The approximations of remaining 98 ($5^{3}$ - 27) grid points can be obtained by using tri-quadratic interpolation with the known 27-node (8 corner nodes, 12 midpoints of edges, 6 centers of faces and 1 center of the hexahedral element) values.

The tri-quartic interpolation formula at natural coordinates ($\xi,\eta,\zeta$) is defined as

[TABLE]

where the shape functions $N_{i}$ can be written below:

[TABLE]

where $l_{i}^{2}$(x) (0 $\leq$ i $\leq$ 2) are the Lagrange basis polynomials of degree 2 which are defined as

[TABLE]

and ($\xi$, $\eta$, $\zeta$) is the natural coordinate of node i (1 $\leq$ i $\leq$ 27).

5 Bi-Conjugate Gradient Method

The Bi-Conjugate Gradient (Bi-CG) method is an algorithm which is focus on solving linear equation systems

[TABLE]

Compared with the Conjugate Gradient (CG) method which needs matrix A to be self-joint, the Bi-CG method does not require matrix A self-joint but require it to multiply conjugate transpose $A^{*}$. In addition, the Bi-CG method replaces the residual’s orthogonal sequence in the CG method with two sequences which are mutually orthogonal. In the Bi-CG method, the residual $r_{j}$ is orthogonal with a set of vectors $\hat{r}_{0}$, $\hat{r}_{1}$ … $\hat{r}_{j-1}$ and $\hat{r}_{j}$ is also orthogonal with $r_{0}$, $r_{1}$ … $r_{j-1}$. These relationships can be achieved by two three-term recurrence relations of vectors \{$$\hat{r}_{j}$$\} and \{$$r_{j}$$\}. Meanwhile, the Bi-CG method terminates within at most n steps when A is an n by n matrix. The preconditioned version and the unpreconditioned version of algorithms of the Bi-CG method are described as follows:

In the Bi-CG Method with the preconditioner algorithm below, $x_{k}^{*}$ is adjoint, $\bar{\alpha}$ is the complex conjugate and the calculated $r_{k}$ and $r_{k}^{*}$ satisfy the following equations respectively

[TABLE]

In this paper, we adopt the Bi-CG method with the preconditioner as the iteration solver in our new EXCMG method.

6 Numerical Experiments

Test Problem 1. The exact solution of the test problem 1 introduced in [19] can be written as

[TABLE]

Applying the biharmonic operator on the exact solution, we can obtain the forcing term f(x, y, z) as follows:

[TABLE]

Obtain first boundary data from the exact solution while obtaining second boundary data $f_{2}$ by taking partial derivative for the exact solution.

Results listed on the table 1 of numerical experiments are performed with EXCMGbi-cg, using 3.6 thousand unknowns on the coarsest grid 32$\times$32$\times$32 and more than 135 million unknowns on the finest grid 512$\times$512$\times$512. In the table 1, “Iter” denotes the number of iterations needed for the Bi-CG solver to achieve the relative residual less than the given tolerance. Additionally, the last row in the table provides the the $L^{\infty}$-error and $L^{2}$-error of the extrapolated solution $u_{h}$ on the finest grid, and the amount of computational cost of the ${EXCMG}_{bi-cg}$ method in terms of a work unit (WU) on the finest grid, which is defined as the total computation required to perform one relaxation sweep on the finest grid. We use the same notations in all tables.

From the results in the table 1, it is clear that the numerical solution $u_{h}$ reaches almost full second-order accuracy while the initial guess $w_{h}$ is third-order approximation to numerical solution $u_{h}$. In addition, the extrapolated solution $\tilde{u}_{h}$ increase the numerical solution’s accuracy greatly. What’s more, the number of iterations is reduced significantly while the grids are finer and finer and this feature is especially important while solving large linear systems. We will introduce this feature specifically in the following text.

First, we we define the error ratio $r_{h}$ as

[TABLE]

For the order of $||w_{h}-u_{h}||_{2}$ is one higher than the order of $||u_{h}-u||_{2}$, the error ratio $r_{h}$ is almost $\frac{1}{2^{n}}$ where n denotes the level of the grid. As the grid becomes finer, $w_{h}$ is much closer to $u_{h}$ especially on the finest grid. Therefore, when the grid is fine enough, the error of $||u_{h}-w_{h}||_{2}$ is smaller so that the number of iterations is reduced.

For test problem 1, on the finest grid 512$\times$512$\times$512, the error ratio $r_{h}$ is 0.028. It is obvious that the $r_{h}$ is so small on the finest grid that we only need to perform one iteration on the finest grid. The number of iteration is reduced significantly.

Test Problem 2. The exact solution of the test problem 2 can be written as

[TABLE]

Applying the biharmonic operator on the exact solution, we can obtain the forcing term f(x, y, z) as follows:

[TABLE]

Obtain first boundary data from the exact solution while obtaining second boundary data $f_{2}$ by taking partial derivative for the exact solution.

Again, results listed on the table 2 of numerical experiments are performed on five level grids with 3.6 thousand unknowns on the coarsest grid 32$\times$32$\times$32 and more than 135 million unknowns on the finest grid 512$\times$512$\times$512. Moreover, from table 2, we can see that numerical solution $u_{h}$ reaches almost full second-order accuracy, the initial guess $w_{h}$ is third-order approximation to numerical solution $u_{h}$, while the extrapolated solution $\tilde{u}_{h}$ increase the numerical solution’s accuracy significantly. If we use $w_{h}$ as numerical solution on the finest grid 512$\times$512$\times$512, the error ratio $r_{h}$ is already 0.27. Thus only six iterations are needed to perform to achieve the expected accuracy.

Test Problem 3. The exact solution of the test problem 3 can be written as

[TABLE]

Applying the biharmonic operator on the exact solution, we can obtain the forcing term f(x, y, z) as follows:

[TABLE]

Obtain first boundary data from the exact solution while obtaining second boundary data $f_{2}$ by taking partial derivative for the exact solution.

Again, five level grids are used with 3.6 thousand unknowns on the coarsest grid 32$\times$32$\times$32 and more than 135 million unknowns on the finest grid 512$\times$512$\times$512. In addition, from table 3, we can see that numerical solution $u_{h}$ reaches almost full second-order accuracy, the initial guess $w_{h}$ is third-order approximation to numerical solution $u_{h}$, while the extrapolated solution $\tilde{u}_{h}$ increase the numerical solution’s accuracy greatly. On the finest grid 512$\times$512$\times$512, if using $w_{h}$ as numerical solution, the error ratio $r_{h}$ is already equal to 0.13. Thus we only need to perform six iterations to achieve the expected accuracy.

Test Problem 4. The exact solution of the test problem 4 can be written as

[TABLE]

Applying the biharmonic operator on the exact solution, we can obtain the forcing term f(x, y, z) as follows:

[TABLE]

Obtain first boundary data from the exact solution while obtaining second boundary data $f_{2}$ by taking partial derivative for the exact solution.

Again, results listed on the table 4 of numerical experiments are performed with EXCMGbi-cg, using 3.6 thousand unknowns on the coarsest grid 32$\times$32$\times$32 and more than 135 million unknowns on the finest grid 512$\times$512$\times$512. Besides, from table 4, we can see that numerical solution $u_{h}$ reaches almost full second-order accuracy, the initial guess $w_{h}$ is third-order approximation to numerical solution $u_{h}$, while the extrapolated solution $\tilde{u}_{h}$ increase the numerical solution’s accuracy greatly. On the finest grid 512$\times$512$\times$512, using $w_{h}$ as numerical solution, the error ratio $r_{h}$ is already equal to 0.090. Thus we only need to perform iterations eight times to achieve the expected accuracy.

Test Problem 5. The exact solution of the test problem 5 can be written as

[TABLE]

Obtain first boundary data from the exact solution while obtaining second boundary data $f_{2}$ by taking partial derivative for the exact solution.

Again, we use five level grids which have 3.6 thousand unknowns on the coarsest grid 32$\times$32$\times$32 and more than 135 million unknowns on the finest grid 512$\times$512$\times$512. Additionally, from table 5, we can see that numerical solution $u_{h}$ reaches almost full second-order accuracy, the initial guess $w_{h}$ is third-order approximation to numerical solution $u_{h}$, while the extrapolated solution $\tilde{u}_{h}$ increase the numerical solution’s accuracy significantly.

7 Conclusion

In this work, we propose a new extrapolation cascadic multigrid method $EXCMG_{bi-cg}$ to solve the linear three-dimensional biharmonic equation. By applying the Richardson extrapolation and quadratic interpolation methods on numerical solutions which are on current and previous grids, much better initial guesses of iterative solutions are obtained on the next finer grid so that the iterative time for Bi-CG solver is reduced. It is the main advantage of our work. Additionally, the introduction of the relative residual tolerance into our work enables us to obtain the desired accuracy conveniently. Furthermore, reducing computational time and the number of iteration, the numerical results of tests demonstrate that the $EXCMG_{bi-cg}$ method is efficient and particularly suitable for solving large scale problems.