# Smoothing Analysis of Distributive Red-Black Jacobi Relaxation for Solving 2D Stokes Flow by Multigrid Method

Authors: Xingwen Zhu; Lixiang Zhang
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Smoothing analysis process of distributive red-black Jacobi relaxation in multigrid method for solving 2D Stokes flow is mainly investigated on the nonstaggered grid by using local Fourier analysis (LFA). For multigrid relaxation, the nonstaggered discretizing scheme of Stokes flow is generally stabilized by adding an artificial pressure term. Therefore, an important problem is how to determine the zone of parameter in adding artificial pressure term in order to make stabilization of the algorithm for multigrid relaxation. To end that, a distributive red-black Jacobi relaxation technique for the 2D Stokes flow is established. According to the 2h-harmonics invariant subspaces in LFA, the Fourier representation of the distributive red-black Jacobi relaxation for discretizing Stokes flow is given by the form of square matrix, whose eigenvalues are meanwhile analytically computed. Based on optimal one-stage relaxation, a mathematical relation of the parameter in artificial pressure term between the optimal relaxation parameter and related smoothing factor is well yielded. The analysis results show that the numerical schemes for solving 2D Stokes flow by multigrid method on the distributive red-black Jacobi relaxation have a specified convergence parameter zone of the added  artificial pressure term.


---

## Body

## 1. Introduction

Multigrid methods [1–7] are generally considered as one of the fastest numerical methods which have an optimally computational complexity for solving partial differential equations (PDEs), especially for 3D steady incompressible Newtonian flow governed by Navier-Stokes equations.In multigrid methods, smoothing relaxations play an important role. Several multigrid relaxation methods were developed for solving PDEs, which are roughly classified into two categories, collective and decoupled relaxations [8]. The collective relaxations are considered as a straightforward generalization of the scalar case [2]. The early decoupled relaxation is on a distributive Gauss-Seidel relaxation [9]. Gradually, it is generalized to an incomplete LU factorization relaxation [10]. Recently, Stokes system with distributive Gauss-Seidel relaxation based on the least squares commutator has been researched [11]. Much of the relaxations for Stokes system is seen in [12, 13].For multigrid methods, LFA is a very useful tool to design efficient algorithms and to predict convergence factors for solving PDEs with high order accuracy [1–7]. Distributive relaxation for poroelasticity equations is optimized by LFA [14]. Using LFA, textbook efficiency multigrid solver for compressible Navier-Stokes equations is designed [15]. All-at-once multigrid approach for optimality systems with LFA is discussed in detail, and an analytical expression of the convergence factors is given by using symbolic computation [16–18].The smoothing analysis of the distributive relaxations for solving 2D Stokes flow is investigated with LFA. As we know, the discretizing Stokes flow in computational domain is not stable by means of standard central differencing on nonstaggered grid. Thus, in order to overcome the stability problem, an artificial pressure term is generally added by the method in [1, 2]. The optimal one-stage relaxation parameter and related smoothing factor of the distributive relaxation with the red-black Jacobi point relaxation need to be developed. In deriving an explicit formulation of the smoothing factor for the multigrid method, the symbolic operation process is carried out by using the MATLAB and Mathematica software, especially, by the cylindrical algebraic decomposition (CAD) function in the Mathematica build-in command [19].


## 2. Discretizing Stokes Flow and LFA

### 2.1. Discrete Stokes Flow

3D steady incompressible Newtonian flow governed by Navier-Stokes equations is given as(1)-Δu⃑+∇⃑p=f⃑x,y,z∈Ω,∇⃑·u⃑=0x,y,z∈Ω,u⃑=g⃑x,y,z∈∂Ω,where u⃑=(u(x,y,z),v(x,y,z),w(x,y,z)) is the velocity field, p=p(x,y,z) represents the pressure, f⃑=(f1(x,y,z),f2(x,y,z),f3(x,y,z)) is the external force field, (x,y,z)∈Ω⊆R3, and ∂Ω is the Dirichlet boundary of the computing domain. From (1), 2D Stokes operator is written as(2)L=-Δ0∂x0-Δ∂y∂x∂y0on nonstaggered grid (3)Gh={x⃑=x,y:=k1h,k2h∣k1,k2∈Z2}.Discretizing Stokes operator (2) by means of standard central differencing is given as(4)Lh′=-Δh0∂xh0-Δh∂yh∂xh∂yh0,where h denotes the uniform mesh size and -Δh, ∂xh, and ∂yh are the second-order difference operator with the following discrete stencils:(5)-Δh=1h2-1-14-1-1h,∂xh=12h-101h,∂yh=12h10-1h.From [1], the above nonstaggered schemes (4) are not stable. Stabilization may be achieved by adding an artificial elliptic pressure term -ch2Δh to the continuity equation in (2) [1, 2, 6]. With discrete operator in (5) and parameter c, the discrete Stokes operator is changed to(6)Lh=-Δh0∂xh0-Δh∂yh∂xh∂yh-ch2Δh.


### 2.2. Elements of LFA in Multigrid

In LFA, a current approximation and its corresponding error and residual are represented by a linear combination of certain exponential functions, for example, Fourier modes, which form a unitary basis in space on a bounded infinite grid functions [1–7].From [1, 2], on nonstaggered grid (3), a unitary basis of the Fourier modes is defined by(7)φhθ⃑,x⃑:=exp⁡iθ⃑·x⃑h,in which θ⃑=(θ1,θ2)∈Θ:=(-π,π]2 is called Fourier frequency, x⃑∈Gh, and complex unit i=-1. Thus, a Fourier space is yielded as(8)F(θ⃑):=span⁡{φhθ⃑,x⃑∣θ⃑∈Θ}.From [1–7], applying (3) and (7), for 2D scalar discrete operator Dh with discrete stencil (9)Dh=lk⃑h,where lk⃑∈R and k⃑∈J⊂Z2 containing (0,0); the Fourier mode of (9) is defined by (10)D~hθ⃑:=∑k⃑∈Jlk⃑exp⁡⁡iθ⃑·k⃑,with θ⃑·k⃑=θ1k1+θ2k2, subjected to(11)Dhφhθ⃑,x⃑=D~hθ⃑φhθ⃑,x⃑.A main idea of LFA is to analyze relaxation properties in multigrid for (6) by evaluating their effects on the Fourier components. From [2, 14, 16], if standard coarsening in 2D is selected, each low frequency θ⃑=θ⃑00∈Θlow2h=(-π/2,π/2]2 is coupled with three high frequencies {θ⃑11,θ⃑10,θ⃑01}∈Θhigh2h in the transition from Gh to G2h, where Θhigh2h=Θ∖Θlow2h, and(12)θ⃑α⃑=θ⃑-α1signθ1,α2signθ2π,where α⃑∈Λ={00,11,10,01} and α⃑=(α1,α2) are denoted by (α1,α2):=α1α2. In this paper, standard coarsening is to be assumed. Then, the Fourier space (8) is subdivided into 2h-harmonics subspaces(13)F2hθ⃑≔span⁡φhθ⃑00,x⃑,φhθ⃑11,x⃑,φhθ⃑10,x⃑,φhθ⃑01,x⃑.


## 2.1. Discrete Stokes Flow

3D steady incompressible Newtonian flow governed by Navier-Stokes equations is given as(1)-Δu⃑+∇⃑p=f⃑x,y,z∈Ω,∇⃑·u⃑=0x,y,z∈Ω,u⃑=g⃑x,y,z∈∂Ω,where u⃑=(u(x,y,z),v(x,y,z),w(x,y,z)) is the velocity field, p=p(x,y,z) represents the pressure, f⃑=(f1(x,y,z),f2(x,y,z),f3(x,y,z)) is the external force field, (x,y,z)∈Ω⊆R3, and ∂Ω is the Dirichlet boundary of the computing domain. From (1), 2D Stokes operator is written as(2)L=-Δ0∂x0-Δ∂y∂x∂y0on nonstaggered grid (3)Gh={x⃑=x,y:=k1h,k2h∣k1,k2∈Z2}.Discretizing Stokes operator (2) by means of standard central differencing is given as(4)Lh′=-Δh0∂xh0-Δh∂yh∂xh∂yh0,where h denotes the uniform mesh size and -Δh, ∂xh, and ∂yh are the second-order difference operator with the following discrete stencils:(5)-Δh=1h2-1-14-1-1h,∂xh=12h-101h,∂yh=12h10-1h.From [1], the above nonstaggered schemes (4) are not stable. Stabilization may be achieved by adding an artificial elliptic pressure term -ch2Δh to the continuity equation in (2) [1, 2, 6]. With discrete operator in (5) and parameter c, the discrete Stokes operator is changed to(6)Lh=-Δh0∂xh0-Δh∂yh∂xh∂yh-ch2Δh.


## 2.2. Elements of LFA in Multigrid

In LFA, a current approximation and its corresponding error and residual are represented by a linear combination of certain exponential functions, for example, Fourier modes, which form a unitary basis in space on a bounded infinite grid functions [1–7].From [1, 2], on nonstaggered grid (3), a unitary basis of the Fourier modes is defined by(7)φhθ⃑,x⃑:=exp⁡iθ⃑·x⃑h,in which θ⃑=(θ1,θ2)∈Θ:=(-π,π]2 is called Fourier frequency, x⃑∈Gh, and complex unit i=-1. Thus, a Fourier space is yielded as(8)F(θ⃑):=span⁡{φhθ⃑,x⃑∣θ⃑∈Θ}.From [1–7], applying (3) and (7), for 2D scalar discrete operator Dh with discrete stencil (9)Dh=lk⃑h,where lk⃑∈R and k⃑∈J⊂Z2 containing (0,0); the Fourier mode of (9) is defined by (10)D~hθ⃑:=∑k⃑∈Jlk⃑exp⁡⁡iθ⃑·k⃑,with θ⃑·k⃑=θ1k1+θ2k2, subjected to(11)Dhφhθ⃑,x⃑=D~hθ⃑φhθ⃑,x⃑.A main idea of LFA is to analyze relaxation properties in multigrid for (6) by evaluating their effects on the Fourier components. From [2, 14, 16], if standard coarsening in 2D is selected, each low frequency θ⃑=θ⃑00∈Θlow2h=(-π/2,π/2]2 is coupled with three high frequencies {θ⃑11,θ⃑10,θ⃑01}∈Θhigh2h in the transition from Gh to G2h, where Θhigh2h=Θ∖Θlow2h, and(12)θ⃑α⃑=θ⃑-α1signθ1,α2signθ2π,where α⃑∈Λ={00,11,10,01} and α⃑=(α1,α2) are denoted by (α1,α2):=α1α2. In this paper, standard coarsening is to be assumed. Then, the Fourier space (8) is subdivided into 2h-harmonics subspaces(13)F2hθ⃑≔span⁡φhθ⃑00,x⃑,φhθ⃑11,x⃑,φhθ⃑10,x⃑,φhθ⃑01,x⃑.


## 3. Smoothing Process

### 3.1. Distributive Relaxation of System (6)

From [1, 2, 7], a distributive operator for the discrete system (6) is constructed as(14)Ch=Ih0-∂xh0Ih-∂yh00-Δh,where Ih is the unit operator with discrete stencil [1]h. From (14), the discrete system (6) is transformed as(15)LhCh=-Δh000-Δh0∂xh∂yhch2Δh2-Δ2h,where the discrete stencils of Δh2 and -Δ2h are (16)Δh2=1h412-821-820-812-821h,-Δ2h=14h200-10000000-1040-10000000-100h=14h2-1-14-1-12h.From (9)–(11), the Fourier modes of the scalar discrete operators of (16) are(17)Δ~h2(θ⃑)=-Δ~hθ⃑2,-Δ~2h(θ⃑)=-∂~xhθ⃑2-∂~yhθ⃑2,where(18)-Δ~hθ⃑=1h24-exp⁡⁡-iθ1-exp⁡⁡iθ1-exp⁡⁡(-iθ2)-exp⁡⁡(iθ2)=1h2(4-2cos⁡θ1-2cos⁡θ2),(19)∂~xhθ⃑=12hexp⁡iθ1-exp⁡-iθ1=1hisinθ1,(20)∂~yhθ⃑=12hexp⁡iθ2-exp⁡-iθ2=1hisinθ2.


### 3.2. Optimal One-Stage Relaxation

For the discrete scalar operator of (15), standard coarsening and an ideal coarse grid correction operator [2] are applied as(21)Qh2hF2h(θ⃑)=:Q⌢h2h=diag0,1,1,1∈C4×4,where Q⌢h2h is the Fourier representation of the operator Qh2h with subspace (13), which suppresses the low frequency error components and makes the high frequency components unchanged. Then, from [2], the smoothing factor for discrete operator (9) is defined by (22)ρ(n,Dh)=sup⁡θ⃑∈ΘlowρQ⌢h2hS⌢hθ⃑,ωn1/n.It implies that the asymptotic error reduction of the high frequency error components is given by n sweeps of the relaxation method, where S⌢h(θ⃑,ω) is the Fourier representation of the relaxation operator Sh(ω) on subspace (13) and ω is the relaxation parameter.From [2, 14], a good smoothing factor is obtained by using one-stage parameter ω in the relaxation operator Sh(ω); the optimal smoothing factor and related smoothing parameter are given by(23)ωopt=22-Smax⁡-Smin⁡,ρopt=Smax⁡-Smin⁡2-Smax⁡-Smin⁡,where Smax⁡ and Smin⁡ are the max and min eigenvalues of the product matrix Q⌢h2hS⌢h(θ⃑,1) with the relaxation parameter ω=1 for θ∈Θlow2h. From [2, 19], the smoothing factor of (6) with the distributive relaxation (14) is determined by the diagonal blocks of the transformed system (15), which is given by(24)ρn,Lh=max⁡ρn,-Δh,ρn,ch2Δh2-Δ2h.


### 3.3. Optimal Smoothing for Stokes Flow

The red-black Jacobi point relaxationShRB is applied to (15) to discuss the optimal smoothing problems for Stokes flow. From [1, 2, 14], the operator ShRB makes the 2h-harmonics subspace (13) invariant; that is,(25)ShRBF2h(θ⃑)=:S⌢hRBθ⃑∈C4×4,where S⌢hRB(θ⃑) is the Fourier representation of ShRB(ω) with relaxation parameter ω=1 and is given as (26)S⌢hRBθ⃑,1=S⌢hRBθ⃑=S⌢hBθ⃑S⌢hRθ⃑=12A00+1-A11+100-A00-1A11+10000A10+1-A01+100-A10+1A01+1·12A00+1A11-100A00-1A11+10000A10+1A01-100A10-1A01+1,in which(27)Aα⃑=1-D~h(θ⃑α⃑)D~h0θ⃑α⃑denotes the Fourier mode of the point Jacobi relaxation for the discrete operator (9) on subspace (13) and D~h0(θ⃑α⃑) is the Fourier mode of the discrete operator with the stencil [l(0,0)]h in (9). For the sake of convenient discussion in the following, two variables are introduced as(28)s1=sin2θ102=sin2θ12,s2=sin2θ202=sin2θ22.Thus, θ⃑=(θ1,θ2)∈Θlow2h=(-π/2,π/2]2 is transformed to s⃑=(s1,s2)∈Slow=[0,1/2]2.Theorem 1.
For the Poisson operator-Δh, the optimal one-stage relaxation parameter and related smoothing factor of the red-black Jacobi point relaxation are stated as(29)ωopt=1615,ρopt=15.Proof.
For the red-black Jacobi point relaxation for the Poisson operatorDh=-Δh, substituting (12), (18), and (28) into (26) and (27), and from (5), the product of (21) and (25) is written as(30)Q⌢h2hS⌢hθ⃑,1=Q⌢h2hS⌢hθ⃑=120000s1+s21-s1-s2s1+s2s1+s2-10000s1-s2s1-s2+1s2-s1s1-s2+100s1-s2s2-s1+1s2-s1s2-s1+1.Thus, eigenvalues of (30) are obtained as(31)λ1=0,λ2=s1-s22,λ3=0,λ4=(s1+s2)(s1+s2-1)2.From (31), the max and min eigenvalues of (30) are yielded as(32)Smax⁡=max⁡(s1,s2)∈0,1/22{λ1,λ2,λ3,λ4}=max⁡(s1,s2)∈0,1/22λ2=14,Smin⁡=min⁡(s1,s2)∈0,1/22{λ1,λ2,λ3,λ4}=min⁡(s1,s2)∈0,1/22λ4=-18.From (23) and (32), (29) is obtained. Theorem 1 holds.Next,ρ(n,ch2Δh2-Δ2h) for the red-black Jacobi point relaxation need to be computed. Meanwhile, the smoothing factor of distributive relaxation (15) is given as follows.Theorem 2.
For the discrete operatorch2Δh2-Δ2h with c>0, the optimal one-stage relaxation parameter and related smoothing factor of the red-black Jacobi point relaxation are given by(33)ωopt=1+20c1+16c0<c≤13221+20c21+56c+1744c2132<c≤112,ρopt=11+16c0<c≤1321+24c+1104c21+56c+1744c2132<c≤112.Proof.
For the discrete operator(34)Dh=ch2Δh2-Δ2h,from (17)–(20), the Fourier mode of (34) is given by (35)D~hθ⃑=1h24c2-cos⁡θ1-cos⁡θ22+sin2θ1+sin2θ2.Thus, when the red-black point relaxation is applied to (34), from (16), substituting (12), (28), and (35) into (26) and (27), the product of (21) and (25) is(36)Q⌢h2hS⌢hRBθ⃑,1=Q⌢h2hS⌢hRBθ⃑=14diagR11,R22,where both R11 and R22 are 2 × 2 square matrices, whose expressions are below:(37)R11=0001·A00B+1-A11B+1-A00B-1A11B+1·A00R+1A11R-1A00R-1A11R+1=41+20c2·001+4-1-36c+16cs1+s2·s1-s12+s2-s22+4cs1+s2264c-1+s1+s2s1-s12+s2-s22+4c-2+s1+s22,R22=1001·A10B+1-A01B+1-A10B+1A01B+1·A10R+1A01R-1A10R-1A01R+1=41+20c2·1+24c+80c2-48c+4s1+s2+384c2s1-s2+192c2+4s12+s22-256c2-64cs13-s23-64cs2+128cs22+32cs1s22s2-2s1-11-12c-64cs1-s24cs1-s22+8cs1-s2+4c+s1-s12+s2-s2264cs1-s24cs1-s22-8cs1-s2+4c+s1-s12+s2-s221+24c+80c2-48c+4s1+s2-384c2s1-s2+4+192c2s12+s22+256c2-64cs13-s23-64cs1+128cs12+32cs1s22s1-2s2+11-12c.Thus, the eigenvalues of matrix (36) are obtained as(38)λ1=0,λ2=64c1+20c2-1+s1+s2s1-s12+s2-s22+4c-2+s1+s22,(39)λ3,4=11+20c21+24c+80c2-(4+80c)(s1+s2)+(4+64c+192c2)(s12+s22)+(32c-384c2)s1s2±32c(s1-s2)1+80c2+24c+(-64c2+64c+4)(s12+s22)-(80c+4)(s1+s2)+(128c2+32c)s1s2.By using the MATLAB and Mathematica software with cylindrical algebraic decomposition function [19], for s⃑=(s1,s2)∈(0,1/2)2, there is no extreme value for (39); when 0<c≤1/32, one of extreme values of (38) is obtained as(40)s1*=64c2+3+40c-348c-6,s2*=64c2+3+40c-348c-6.Thus, for s⃑∈Slow=[0,1/2]2, besides (40), the possible extreme values of the eigenvalues of matrix (36) are placed on the boundary of Slow. From -1≤λk≤1 with k=1,…,4, then 0<c≤1/12. Noting that (40) exists with 0<c≤1/32. From (38)–(40), when 0<c≤1/12, for s⃑∈Slow, the max and min eigenvalues of (36) are yielded as(41)Smax⁡=λ3,4(0,0)=1+4c1+20c,Smin⁡=λ3,412,12=4c-11+20c0<c≤132λ2(0,0)=-32c20c+12132<c≤112.Substituting (41) into (23), (33) is obtained. Theorem 2 holds.From (33), 1/2≤ρopt(ch2Δh2-Δ2h)<1 holds with 0<c≤1/12. Therefore, from Theorems 1 and 2, when 0<c≤1/12, the smoothing factor of (6) with the distributive relaxation (14) is as(42)12≤ρoptLh=max⁡⁡ρopt-Δh,ρoptch2Δh2-Δ2h=ρopt(ch2Δh2-Δ2h)<1.


## 3.1. Distributive Relaxation of System (6)

From [1, 2, 7], a distributive operator for the discrete system (6) is constructed as(14)Ch=Ih0-∂xh0Ih-∂yh00-Δh,where Ih is the unit operator with discrete stencil [1]h. From (14), the discrete system (6) is transformed as(15)LhCh=-Δh000-Δh0∂xh∂yhch2Δh2-Δ2h,where the discrete stencils of Δh2 and -Δ2h are (16)Δh2=1h412-821-820-812-821h,-Δ2h=14h200-10000000-1040-10000000-100h=14h2-1-14-1-12h.From (9)–(11), the Fourier modes of the scalar discrete operators of (16) are(17)Δ~h2(θ⃑)=-Δ~hθ⃑2,-Δ~2h(θ⃑)=-∂~xhθ⃑2-∂~yhθ⃑2,where(18)-Δ~hθ⃑=1h24-exp⁡⁡-iθ1-exp⁡⁡iθ1-exp⁡⁡(-iθ2)-exp⁡⁡(iθ2)=1h2(4-2cos⁡θ1-2cos⁡θ2),(19)∂~xhθ⃑=12hexp⁡iθ1-exp⁡-iθ1=1hisinθ1,(20)∂~yhθ⃑=12hexp⁡iθ2-exp⁡-iθ2=1hisinθ2.


## 3.2. Optimal One-Stage Relaxation

For the discrete scalar operator of (15), standard coarsening and an ideal coarse grid correction operator [2] are applied as(21)Qh2hF2h(θ⃑)=:Q⌢h2h=diag0,1,1,1∈C4×4,where Q⌢h2h is the Fourier representation of the operator Qh2h with subspace (13), which suppresses the low frequency error components and makes the high frequency components unchanged. Then, from [2], the smoothing factor for discrete operator (9) is defined by (22)ρ(n,Dh)=sup⁡θ⃑∈ΘlowρQ⌢h2hS⌢hθ⃑,ωn1/n.It implies that the asymptotic error reduction of the high frequency error components is given by n sweeps of the relaxation method, where S⌢h(θ⃑,ω) is the Fourier representation of the relaxation operator Sh(ω) on subspace (13) and ω is the relaxation parameter.From [2, 14], a good smoothing factor is obtained by using one-stage parameter ω in the relaxation operator Sh(ω); the optimal smoothing factor and related smoothing parameter are given by(23)ωopt=22-Smax⁡-Smin⁡,ρopt=Smax⁡-Smin⁡2-Smax⁡-Smin⁡,where Smax⁡ and Smin⁡ are the max and min eigenvalues of the product matrix Q⌢h2hS⌢h(θ⃑,1) with the relaxation parameter ω=1 for θ∈Θlow2h. From [2, 19], the smoothing factor of (6) with the distributive relaxation (14) is determined by the diagonal blocks of the transformed system (15), which is given by(24)ρn,Lh=max⁡ρn,-Δh,ρn,ch2Δh2-Δ2h.


## 3.3. Optimal Smoothing for Stokes Flow

The red-black Jacobi point relaxationShRB is applied to (15) to discuss the optimal smoothing problems for Stokes flow. From [1, 2, 14], the operator ShRB makes the 2h-harmonics subspace (13) invariant; that is,(25)ShRBF2h(θ⃑)=:S⌢hRBθ⃑∈C4×4,where S⌢hRB(θ⃑) is the Fourier representation of ShRB(ω) with relaxation parameter ω=1 and is given as (26)S⌢hRBθ⃑,1=S⌢hRBθ⃑=S⌢hBθ⃑S⌢hRθ⃑=12A00+1-A11+100-A00-1A11+10000A10+1-A01+100-A10+1A01+1·12A00+1A11-100A00-1A11+10000A10+1A01-100A10-1A01+1,in which(27)Aα⃑=1-D~h(θ⃑α⃑)D~h0θ⃑α⃑denotes the Fourier mode of the point Jacobi relaxation for the discrete operator (9) on subspace (13) and D~h0(θ⃑α⃑) is the Fourier mode of the discrete operator with the stencil [l(0,0)]h in (9). For the sake of convenient discussion in the following, two variables are introduced as(28)s1=sin2θ102=sin2θ12,s2=sin2θ202=sin2θ22.Thus, θ⃑=(θ1,θ2)∈Θlow2h=(-π/2,π/2]2 is transformed to s⃑=(s1,s2)∈Slow=[0,1/2]2.Theorem 1.
For the Poisson operator-Δh, the optimal one-stage relaxation parameter and related smoothing factor of the red-black Jacobi point relaxation are stated as(29)ωopt=1615,ρopt=15.Proof.
For the red-black Jacobi point relaxation for the Poisson operatorDh=-Δh, substituting (12), (18), and (28) into (26) and (27), and from (5), the product of (21) and (25) is written as(30)Q⌢h2hS⌢hθ⃑,1=Q⌢h2hS⌢hθ⃑=120000s1+s21-s1-s2s1+s2s1+s2-10000s1-s2s1-s2+1s2-s1s1-s2+100s1-s2s2-s1+1s2-s1s2-s1+1.Thus, eigenvalues of (30) are obtained as(31)λ1=0,λ2=s1-s22,λ3=0,λ4=(s1+s2)(s1+s2-1)2.From (31), the max and min eigenvalues of (30) are yielded as(32)Smax⁡=max⁡(s1,s2)∈0,1/22{λ1,λ2,λ3,λ4}=max⁡(s1,s2)∈0,1/22λ2=14,Smin⁡=min⁡(s1,s2)∈0,1/22{λ1,λ2,λ3,λ4}=min⁡(s1,s2)∈0,1/22λ4=-18.From (23) and (32), (29) is obtained. Theorem 1 holds.Next,ρ(n,ch2Δh2-Δ2h) for the red-black Jacobi point relaxation need to be computed. Meanwhile, the smoothing factor of distributive relaxation (15) is given as follows.Theorem 2.
For the discrete operatorch2Δh2-Δ2h with c>0, the optimal one-stage relaxation parameter and related smoothing factor of the red-black Jacobi point relaxation are given by(33)ωopt=1+20c1+16c0<c≤13221+20c21+56c+1744c2132<c≤112,ρopt=11+16c0<c≤1321+24c+1104c21+56c+1744c2132<c≤112.Proof.
For the discrete operator(34)Dh=ch2Δh2-Δ2h,from (17)–(20), the Fourier mode of (34) is given by (35)D~hθ⃑=1h24c2-cos⁡θ1-cos⁡θ22+sin2θ1+sin2θ2.Thus, when the red-black point relaxation is applied to (34), from (16), substituting (12), (28), and (35) into (26) and (27), the product of (21) and (25) is(36)Q⌢h2hS⌢hRBθ⃑,1=Q⌢h2hS⌢hRBθ⃑=14diagR11,R22,where both R11 and R22 are 2 × 2 square matrices, whose expressions are below:(37)R11=0001·A00B+1-A11B+1-A00B-1A11B+1·A00R+1A11R-1A00R-1A11R+1=41+20c2·001+4-1-36c+16cs1+s2·s1-s12+s2-s22+4cs1+s2264c-1+s1+s2s1-s12+s2-s22+4c-2+s1+s22,R22=1001·A10B+1-A01B+1-A10B+1A01B+1·A10R+1A01R-1A10R-1A01R+1=41+20c2·1+24c+80c2-48c+4s1+s2+384c2s1-s2+192c2+4s12+s22-256c2-64cs13-s23-64cs2+128cs22+32cs1s22s2-2s1-11-12c-64cs1-s24cs1-s22+8cs1-s2+4c+s1-s12+s2-s2264cs1-s24cs1-s22-8cs1-s2+4c+s1-s12+s2-s221+24c+80c2-48c+4s1+s2-384c2s1-s2+4+192c2s12+s22+256c2-64cs13-s23-64cs1+128cs12+32cs1s22s1-2s2+11-12c.Thus, the eigenvalues of matrix (36) are obtained as(38)λ1=0,λ2=64c1+20c2-1+s1+s2s1-s12+s2-s22+4c-2+s1+s22,(39)λ3,4=11+20c21+24c+80c2-(4+80c)(s1+s2)+(4+64c+192c2)(s12+s22)+(32c-384c2)s1s2±32c(s1-s2)1+80c2+24c+(-64c2+64c+4)(s12+s22)-(80c+4)(s1+s2)+(128c2+32c)s1s2.By using the MATLAB and Mathematica software with cylindrical algebraic decomposition function [19], for s⃑=(s1,s2)∈(0,1/2)2, there is no extreme value for (39); when 0<c≤1/32, one of extreme values of (38) is obtained as(40)s1*=64c2+3+40c-348c-6,s2*=64c2+3+40c-348c-6.Thus, for s⃑∈Slow=[0,1/2]2, besides (40), the possible extreme values of the eigenvalues of matrix (36) are placed on the boundary of Slow. From -1≤λk≤1 with k=1,…,4, then 0<c≤1/12. Noting that (40) exists with 0<c≤1/32. From (38)–(40), when 0<c≤1/12, for s⃑∈Slow, the max and min eigenvalues of (36) are yielded as(41)Smax⁡=λ3,4(0,0)=1+4c1+20c,Smin⁡=λ3,412,12=4c-11+20c0<c≤132λ2(0,0)=-32c20c+12132<c≤112.Substituting (41) into (23), (33) is obtained. Theorem 2 holds.From (33), 1/2≤ρopt(ch2Δh2-Δ2h)<1 holds with 0<c≤1/12. Therefore, from Theorems 1 and 2, when 0<c≤1/12, the smoothing factor of (6) with the distributive relaxation (14) is as(42)12≤ρoptLh=max⁡⁡ρopt-Δh,ρoptch2Δh2-Δ2h=ρopt(ch2Δh2-Δ2h)<1.


## 4. Conclusions

The smoothing analysis process of the distributive red-black Jacobi point relaxation for solving 2D Stokes flow is analytically presented. Applying (28), the Fourier modes with the trigonometric functions for the discrete operator and relaxation are mapped to rational functions. So, it is possible to apply the cylindrical algebraic decomposition function in the Mathematica software to realize complex smoothing analysis, and the computation process is simplified. The analytical expressions of the smoothing factor for the distributive red-black Jacobi point relaxation are obtained, which is an upper bound for the smoothing rates and is independent of the mesh size with the parameter c. Obviously, it is valuable to understand numerical experiments in multigrid method.



---

*Source: 572198-2015-03-15.xml*