Next Article in Journal
Mapping Topology of Skyrmions and Fractional Quantum Hall Droplets to Nuclear EFT for Ultra-Dense Baryonic Matter
Previous Article in Journal
A Subgradient-Type Extrapolation Cyclic Method for Solving an Equilibrium Problem over the Common Fixed-Point Sets
Previous Article in Special Issue
Analytical Investigation of Fractional-Order Korteweg–De-Vries-Type Equations under Atangana–Baleanu–Caputo Operator: Modeling Nonlinear Waves in a Plasma and Fluid
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Efficient Red–Black Skewed Modified Accelerated Arithmetic Mean Iterative Method for Solving Two-Dimensional Poisson Equation

by
Azali Saudi
1,* and
A’qilah Ahmad Dahalan
2,3
1
Faculty of Computing and Informatics, Universiti Malaysia Sabah, Kota Kinabalu 88450, Malaysia
2
Department of Mathematics, Universiti Pertahanan Nasional Malaysia, Kuala Lumpur 57000, Malaysia
3
CONFIRM Centre for SMART Manufacturing, University of Limerick, V94 C928 Limerick, Ireland
*
Author to whom correspondence should be addressed.
Symmetry 2022, 14(5), 993; https://doi.org/10.3390/sym14050993
Submission received: 15 March 2022 / Revised: 3 April 2022 / Accepted: 10 April 2022 / Published: 12 May 2022
(This article belongs to the Special Issue Advance in Partial Differential Equations of Applied Mathematics)

Abstract

:
This paper presents the extended variants to the established two-stage Arithmetic Mean (AM) method known as the Modified Accelerated Arithmetic Mean (MAAM) and Skewed Modified Accelerated Arithmetic Mean (SkMAAM) methods to solve the two-dimensional elliptic problem. The existing two-stage AM and its skewed variants apply one weighted parameter for the computation of nodes in Levels 1 and 2. The suggested MAAM and SkMAAM methods employ red–black ordering with two different weighted parameters and an additional two distinct accelerated parameters for red and black nodes, respectively. By carefully choosing optimum parameter values, the proposed MAAM and SkMAAM improve the computational execution of the algorithm. With red–black ordering, the computational molecules of red and black nodes are symmetrical, in which the computation of red nodes applies the updated values of their four neighbouring black nodes and vice versa. These symmetrical computational molecules of red and black nodes can be seen for the modified variants MAM and MAAM, and their corresponding skewed variants SkMAM and SkMAAM. The proposed MAAM and SkMAAM methods are compared to the existing AM and Modified AM (MAM) and their corresponding skewed variants, namely the Skewed AM (SkAM) and Skewed MAM (SkMAM) methods. The performance of the newly proposed MAAM and SkMAAM methods is compared against the existing methods in terms of computational complexity and actual execution time. It is shown in the simulation results that the skewed variants are superior to their corresponding regular variants, in which the SkMAAM method gives the best performance.

1. Introduction

The two-stage iterative methods had been studied extensively in the past to solve linear system problems [1]. In [2], a two-stage method known as the Arithmetic Mean (AM) method was introduced. The AM method was applied recently to solve a linear Fredholm integral equation [3]. In [4], a variant of a two-stage approach based on the Alternating Decomposition Explicit method was explored to solve the diffusion equation. An application of the two-stage arithmetic method to solve the glioma growth model was presented in [5]. A combination of AM and Newton methods was reported in [3]. In [6], Sulaiman et al. introduced the Half-Sweep Arithmetic Mean method, in which the half-sweep approach has been successfully applied on the skewed grid to reduce the computational complexity drastically. The existing variants of the AM method employ a weighted parameter to speed up the convergence rate. This sole weighted parameter method, also widely known as the Successive Overrelaxation (SOR) method, was extensively studied by Young [7] to solve the partial differential equation of elliptic type and its impact in the 20th century was presented in [8].
In recent studies, the SOR and its variants were applied to solve problems in image processing [9] and agent navigation [10]. An extension to the SOR known as the Accelerated Overrelaxation (AOR) method that employs two acceleration parameters was later introduced by Hadjidimos [11]. As pointed out in [11], if fully exploited, the presence of two parameters in AOR and its variant [12,13] will provide a faster convergence rate than any other methods of the same type. Moreover, the red–black ordering strategy was applied in a previous work where two different weighted parameters were used for the calculation of red and black nodes to improve the execution time [14,15]. It is interesting to note that the numerical solutions to the fractional version of the partial differential equation were actively studied in recent works. Ali and Abdullah [16] proposed an implicit finite difference method to solve a one-dimensional wave equation. In [17], an efficient computational approach for a local fractional Poisson equation in fractal media was suggested. In future work, we aim to consider employing the two-stage iterative methods for solving fractional partial differential equations.
In this study, we introduced two variants of AM methods, namely the MAAM and SkMAAM methods, for solving elliptic partial differential equations in two-dimensional space. We applied the AOR scheme to the existing AM method to construct a new variant of the AM method known as the Modified Accelerated Arithmetic Method (MAAM) method, which is implemented on a regular grid using a red–black ordering strategy. Apart from MAAM, its variant, namely the Skewed Modified Accelerated Arithmetic Method (SkMAAM), was implemented on the skewed grid. It is widely known that the half-sweep approach can be applied to the skewed grid, which is capable of reducing the amount of computation by half.
In our previous work, the half-sweep variants of the AM and MAM methods, known as Half-Sweep AM (HSAM) [18] and Half-Sweep MAM (HSMAM) [19], were applied to solve the two-dimensional elliptic equation. The performance of the proposed MAAM and SkMAAM methods was compared with the AM and MAM methods, and the current state-of-the-art HSAM method (also known as the SkAM method) and SkMAM method. The encouraging results given by the two-stage iterative methods motivate us to develop an extension to the previous work. All skewed variants considered in this study apply the rotated grid to reduce the amount of computation. Thus, the main contribution of this work is that the newly proposed MAAM and SkMAAM methods introduce the use of accelerated parameters that are shown to provide faster convergence when optimum values are chosen.
The description of the finite difference methods of SOR and AOR is provided in Section 2. In Section 3, we present the description of the existing AM method and its skewed variant in two-dimensional space implementation. The MAM and SkMAM methods are also presented in Section 3. In Section 4, we describe the proposed MAAM and SkMAAM methods and explain the implementation of both methods for computing the solution of the elliptic equation. The novelty of the proposed methods by introducing additional different accelerated parameters for red and black nodes is explained. The obtained numerical results are provided and summarised in Section 5. The computational complexities of the considered iterative methods are discussed in Section 5. Finally, the conclusions are drawn in Section 6.

2. Finite Difference Approximation Equations

By considering the following elliptic boundary value problem
2 u = 2 u x 2 + 2 u y 2 = f ( x , y ) , ( x , y ) Ω ,
with Dirichlet boundary conditions
u ( x , y ) = g ( x , y ) , ( x , y ) Ω ,
where Ω is a bounded region in 2 . The finite difference method is used to discretise the Poisson Equation (1), which is often applied to model the phenomena of fluid dynamics and heat conduction processes. Equation (1) can be approximated using the standard five-point difference approximation formula given as
u i 1 , j + u i + 1 , j + u i , j 1 + u i , j + 1 4 u i , j = h 2 f i , j ,
where h denotes the spacing between node points in i and j directions of a rectangular grid.
If the axis of the rectangular grid is rotated by 45 [20], we can obtain the skewed five-point difference approximation formula given as
u i 1 , j 1 + u i + 1 , j 1 + u i 1 , j + 1 + u i + 1 , j + 1 4 u i , j = 2 h 2 f i , j .
Equations (3) and (4) are then used for the iterative methods implemented on the regular or skewed grids. Accordingly, the iterative schemes for the standard five-point difference formula on the regular and skewed grids can be written as
u i , j ( k + 1 ) = 1 4 u i 1 , j ( k + 1 ) + u i + 1 , j ( k ) + u i , j 1 ( k + 1 ) + u i , j + 1 ( k ) h 2 f i , j
and
u i , j ( k + 1 ) = 1 4 u i 1 , j 1 ( k + 1 ) + u i + 1 , j 1 ( k + 1 ) + u i 1 , j + 1 ( k ) + u i + 1 , j + 1 ( k ) 2 h 2 f i , j ,
respectively. Consequently, the SOR iterative schemes for the corresponding regular and skewed five-point difference formula can be expressed as
u i , j ( k + 1 ) = ω 4 u i 1 , j ( k + 1 ) + u i + 1 , j ( k ) + u i , j 1 ( k + 1 ) + u i , j + 1 ( k ) h 2 f i , j + ( 1 ω ) u i , j ( k ) ,
u i , j ( k + 1 ) = ω 4 u i 1 , j 1 ( k + 1 ) + u i + 1 , j 1 ( k + 1 ) + u i 1 , j + 1 ( k ) + u i + 1 , j + 1 ( k ) 2 h 2 f i , j + ( 1 ω ) u i , j ( k ) .
As defined by Hadjidimos [11], to obtain the AOR iterative scheme for the standard five-point approximation, we replace u i 1 , j ( k + 1 ) and u i , j 1 ( k + 1 ) with u i 1 , j ( k ) and u i , j 1 ( k ) , respectively, and add the terms r 4 ( u i 1 , j ( k + 1 ) u i 1 , j ( k ) ) and r 4 ( u i , j 1 ( k + 1 ) u i , j 1 ( k ) ) . Hence, the AOR iterative scheme for the standard five-point formula on the regular grid can be written as
u i , j ( k + 1 ) = ω 4 u i 1 , j ( k ) + u i + 1 , j ( k ) + u i , j 1 ( k ) + u i , j + 1 ( k ) h 2 f i , j + ( 1 ω ) u i , j ( k ) + r 4 u i 1 , j ( k + 1 ) u i 1 , j ( k ) + u i , j 1 ( k + 1 ) u i , j 1 ( k ) .
Similarly, the AOR iterative scheme for the five-point difference formula on the skewed grid can be constructed and is given as
u i , j ( k + 1 ) = ω 4 u i 1 , j 1 ( k ) + u i + 1 , j 1 ( k ) + u i 1 , j + 1 ( k ) + u i + 1 , j + 1 ( k ) 2 h 2 f i , j + ( 1 ω ) u i , j ( k ) r 4 u i 1 , j 1 ( k + 1 ) u i 1 , j 1 ( k ) + u i + 1 , j 1 ( k + 1 ) u i + 1 , j 1 ( k ) .
The application of these finite difference Formulae (5)–(10) to approximate problem (1) will generate a large and sparse linear system that can be stated in matrix form as
A u = f
where A and f are known, and u is unknown. Matrix A can be decomposed into
A = D L T
where D, L and T are diagonal, strictly lower and strictly upper triangular matrices, respectively. The general iterative scheme of the generated linear system can be written as
D u ( k + 1 ) = L u ( k + 1 ) + T u ( k ) + D u ( k ) + f = L u ( k + 1 ) + ( T + D ) u ( k ) + f .
Accordingly, the SOR iterative scheme can be written as
D u ( k + 1 ) = ω L u ( k + 1 ) + ω T u ( k ) + ( 1 ω ) D u ( k ) + ω f = ω L u ( k + 1 ) + ( ω T + ( 1 ω ) D ) u ( k ) + ω f .
The AOR formula may be obtained from (14) by replacing ω L u ( k + 1 ) with ω L u ( k ) and adding the term r L ( u ( k + 1 ) u ( k ) ) as shown below:
D u ( k + 1 ) = ω L u ( k ) + ω T u ( k ) + ( 1 ω ) D u ( k ) + ω f + r L ( u ( k + 1 ) u ( k ) ) = ( ω L + ω T + ( 1 ω ) D ) u ( k ) + ω f + r L ( u ( k + 1 ) u ( k ) ) .
From now on, we shall call ω the weighted parameter and r the acceleration parameter. Furthermore, it can be noted that the AOR can be treated as an extrapolation of the SOR method with overrelaxation parameter r and extrapolation parameter s = ω / r , as described in [11]. From Equations (13)–(15), u ( k + 1 ) denotes the unknown vector at iteration ( k + 1 ) , which we attempt to compute efficiently.
The classical Equation (13) does not provide a convergence optimisation mechanism. In Equation (14), the SOR scheme provides one overrelaxation parameter only to speed up the rate of convergence. Using Equation (15), two parameters can be adjusted to obtain the optimal convergence rate for solving the problem (1). Different values of ω and r can be chosen within the range between 1 and 2, in which good parameter values will give an acceptable optimum number of iterations.

3. Arithmetic Mean (AM) Method

The iterative process of the two-stage AM iterative method involves solving two independent systems such as u ^ ( 1 ) and u ^ ( 2 ) , which are calculated at Levels 1 and 2, respectively. The average sum of the two matrices is computed at Level 3. By applying Equation (14), the general iterative scheme for the AM method can be defined as
D u ^ ( 1 ) = ω ( L u ^ ( 1 ) + T u ( k ) ) + ( 1 ω ) D u ( k ) + ω f , D u ^ ( 2 ) = ω ( L u ( k ) + T u ^ ( 2 ) ) + ( 1 ω ) D u ( k ) + ω f , u ( k + 1 ) = 1 2 u ^ ( 1 ) + u ^ ( 2 ) ,
where the optimal value of the weighted parameter is in the range 1 < ω < 2 . Note that if ω = 1 , the standard Gauss–Seidel method is obtained. From Equation (16), it can be seen that the formulae applied in Levels 1 and 2 are symmetrical, as illustrated in the portion of the computational grid of the regular and skewed cases shown in Figure 1 and Figure 2, respectively.
The iterative method (16) is known to have well-defined structures that are suitable for parallel implementation [21]. The lower triangular system and the upper triangular system in this method can be executed independently on two different processors. The convergence proof of the iterative method is described in the following theorems [21].
Theorem 1.
Let A = ( a i j ) be an N × N irreducibly diagonally dominant matrix with a i j 0 , for all i j and a i i > 0 , i = 1 , 2 , , n . Then, the iterative method (16) is convergent for all 0 < ω < 2 .
Proof of Theorem 1. 
A detailed proof is described in [21].    □
With the AM method, all nodes are involved in the computation, as shown in Figure 3a. The weighted parameter ω is utilised so that an optimal convergence rate can be obtained by carefully choosing a good value in the range between 1 and 2. Note that in Level 1, the iteration continues from lower node to upper node, while in Level 2, the iteration is carried out in reverse order. Level 3 computes the average sum of the updated values of u ( 1 ) and u ( 2 ) obtained from Levels 1 and 2, respectively. The updating process is repeated until the error rate ϵ satisfies the specified maximum tolerance value.
The SkAM method utilises a skewed grid where computational nodes are equally divided by half into black and white colors; see Figure 3b. The following condition is applied in the computation process so that only black nodes are involved during the iteration procedure:
if (i+j) mod 2 is equal to 0 then
    Update u at point ( i , j ) for iteration k
else
    Continue
end if
By using the above condition, all white nodes are ignored in the loop of the SkAM algorithm.

3.1. Red–Black Modified AM (MAM) Method

The MAM method employs two different weighted parameters, α r and α b , for the respective red and black nodes, and is defined as
D u ^ ( 1 ) = α c ( L u ^ ( 1 ) + T u ^ ( 1 ) ) + ( 1 α c ) D u ( k ) + α c f , D u ^ ( 2 ) = α c ( L u ^ ( 2 ) + T u ^ ( 2 ) ) + ( 1 α c ) D u ( k ) + α c f , u ( k + 1 ) = 1 2 u ^ ( 1 ) + u ^ ( 2 ) ,
where c = { r , b } , and α r and α b are the weighted parameters, for the red and black nodes, respectively. The optimal values of the relaxation parameters are in the range 1 < α r , α b < 2 . It is worth noting that if α r = α b = 1 , we will obtain the standard Gauss–Seidel scheme.
All nodes are computed for the MAM method, where red and black nodes are computed in alternating order, as further explained in Section 3.2. The weighted parameters α 1 and α 2 are used for red and black iterations, respectively, to achieve the best convergence rate by carefully selecting good values in the range of 1 to 2. It is worth noting that in Level 1, iteration proceeds from lower to upper node, whereas in Level 2, iteration proceeds in the other direction. Level 3 computes the average of the total of the updated u ^ ( 1 ) and u ^ ( 2 ) values obtained from Levels 1 and 2. The procedure of updating is continued until the error rate ϵ meets the maximum tolerance value.
The SkMAM method is implemented on the skewed grid, where only red and black nodes are involved in the computations. Since the red and black nodes represent 50% of the total nodes in the grid, the computation is reduced approximately by half. All white nodes are ignored during the iteration process.

3.2. Red–Black Ordering

The MAM and SkMAM methods employ red–black ordering, as depicted in Figure 4, where the computational nodes are divided into red and black nodes. For the regular grid applied in the MAM method, half of the nodes are in red color and the other half are in black color. In the case of a skewed grid used for the SkMAM method, half of the nodes are in white color, 25% are red nodes and the remaining 25% are black nodes.
With the red–black strategy, two different weighted parameters, α 1 and α 2 , are applied to substitute parameter ω for the respective red and black nodes. Additionally, the accelerated parameter r is replaced with two distinct accelerated parameters, β 1 and β 2 . By providing a wider choice of values in the range between 1 and 2 and with appropriate parameter values, this procedure will result in an acceptable optimum number of iterations, and thus speed up the computational time. It can be observed clearly in Figure 5 that the computational molecules for red and black nodes are symmetrical, where the black node applies the updated values of its four neighbouring red nodes and vice versa.
By utilising the red–black ordering approach, red and black nodes are computed alternately. One iteration will involve red nodes only, while in the next iteration, only black nodes are computed. For a skewed grid, all white nodes are ignored during the iterations. These white nodes are only computed once after the convergence is achieved.

4. Red–Black Modified Accelerated Arithmetic Mean (MAAM) Method

The formulation of the proposed two-stage MAAM method also employs a red–black strategy, so an independent equation is applied for the red and black nodes. Unlike the MAM method, which employs weighted parameters only, additional accelerated parameters are utilised for this newly proposed method. It can be constructed based on the AOR iterative scheme (15) to solve the given problem and is written as
D u ^ ( 1 ) = α c ( L + T ) u ( k ) + ( 1 α c ) D u ( k ) + α c f + β d ( L ( u ^ ( 1 ) u ( k ) ) + T ( u ^ ( 1 ) u ( k ) ) ) , D u ^ ( 2 ) = α c ( L + T ) u ( k ) + ( 1 α c ) D u ( k ) + α c f + β d ( L ( u ^ ( 2 ) u ( k ) ) + T ( u ^ ( 2 ) u ( k ) ) ) , u ( k + 1 ) = 1 2 u ^ ( 1 ) + u ^ ( 2 ) ,
where c = { r , b } , d = { 1 , 2 } , and α r and α b are the weighted parameters, for the respective red and black nodes, and the constant values β 1 and β 2 are the accelerated parameters at their corresponding Levels 1 and 2. The optimal values for all parameters are in the range 1 < α r , α b , β 1 , β 2 < 2 . If α r = α b = β 1 = β 2 , the original AM method is obtained. Moreover, if α r = α b = β 1 = β 2 = 1 , we obtain the standard Gauss–Seidel method. A proof identical to that given above for Theorem 1 guarantees that the iterative method (18) is convergent for 0 < α r , α b , β 1 , β 2 < 2 .
The iterative schemes for the MAAM are based on Equation (9), and are written as
u i , j ( d ) = α r 4 u i 1 , j ( k ) + u i + 1 , j ( k ) + u i , j 1 ( k ) + u i , j + 1 ( k ) h 2 f i , j + ( 1 α r ) u i , j ( k )
and
u i , j ( d ) = α b 4 u i 1 , j ( k ) + u i + 1 , j ( k ) + u i , j 1 ( k ) + u i , j + 1 ( k ) h 2 f i , j + ( 1 α b ) u i , j ( k ) + β d 4 u i 1 , j ( d ) u i 1 , j ( k ) + u i + 1 , j ( d ) u i + 1 , j ( k ) + u i , j 1 ( d ) u i , j 1 ( k ) + u i , j + 1 ( d ) u i , j + 1 ( k ) ,
for the respective red and black nodes, where d denotes Level 1 or 2. Note that Equation (19) is simply a weighted Jacobi scheme, since it uses the previous values of black nodes. In contrast, Equation (20) utilises the updated values obtained from the prior procedure that computes all red nodes.
Likewise, based on Equation (10), the iterative schemes for the SkMAAM method can be obtained and are written as
u i , j ( d ) = α r 4 u i 1 , j 1 ( k ) + u i + 1 , j 1 ( k ) + u i 1 , j + 1 ( k ) + u i + 1 , j + 1 ( k ) 2 h 2 f i , j + ( 1 α r ) u i , j ( k )
and
u i , j ( d ) = α b 4 u i 1 , j 1 ( k ) + u i + 1 , j 1 ( k ) + u i 1 , j + 1 ( k ) + u i + 1 , j + 1 ( k ) 2 h 2 f i , j + ( 1 α b ) u i , j ( k ) + β d 4 u i 1 , j 1 ( d ) u i 1 , j 1 ( k ) + u i + 1 , j 1 ( d ) u d , j 1 ( k ) + u i 1 , j + 1 ( d ) u i 1 , j + 1 ( k ) + u i + 1 , j + 1 ( d ) u i + 1 , j + 1 ( k ) ,
for the corresponding red and black nodes, where d is Level 1 or 2.
By using the iterative schemes given in Equations (19)–(22), the corresponding implementations of the MAAM and SkMAAM methods are described in Algorithms 1 and 2, respectively.
Algorithm 1: Modified Accelerated Arithmetic Mean (MAAM)
1: k = 0
2: repeat
3:    Level 1:
4:        Compute all red nodes u i , j ( 1 ) using Equation (19) with c = r , d = 1
5:        Compute all black nodes u i , j ( 1 ) using Equation (20) with c = b , d = 1
6:    Level 2:
7:        Compute all red nodes u i , j ( 2 ) using Equation (19) with c = r , d = 2
8:        Compute all black nodes u i , j ( 2 ) using Equation (20) with c = b , d = 2
9:    Level 3:
10:        Compute all nodes
11:          u i , j ( k + 1 ) = 1 2 u i , j ( 1 ) + u i , j ( 2 )
12:     k = k + 1
13: until | | u i , j ( k + 1 ) u i , j ( k ) | | < ϵ
As described in the above MAAM algorithm, all nodes are involved during the iteration procedure, in which parameters α r , α b , β 1 and β 2 are used to optimise the convergence rate, and ϵ is the tolerance error rate. In practice, the optimal values of these parameters are often very close to each other, in which when β 2 β 1 α b α r , the convergence is usually found slightly faster.
Algorithm 2: Skewed Modified Accelerated Arithmetic Mean (SkMAAM)
1: k = 0
2: repeat
3:    Level 1:
4:        Compute all red nodes u i , j ( 1 ) using Equation (21) with c = r , d = 1
5:        Compute all black nodes u i , j ( 1 ) using Equation (22) with c = b , d = 1
6:    Level 2:
7:        Compute all red nodes u i , j ( 2 ) using Equation (21) with c = r , d = 2
8:        Compute all black nodes u i , j ( 2 ) using Equation (22) with c = b , d = 2
9:    Level 3:
10:        Compute all red and black nodes
11:          u i , j ( k + 1 ) = 1 2 u i , j ( 1 ) + u i , j ( 2 )
12:     k = k + 1
13: until | | u i , j ( k + 1 ) u i , j ( k ) | | < ϵ
14: Compute all white nodes using Equation (5)
Similar to the SkMAM algorithm, only black nodes are computed in the above SkMAAM algorithm. As described above, parameters α r , α b , β 1 and β 2 are applied to obtain the optimal convergence rate. By checking if ( i + j ) obtains an even value, all white nodes are ignored during the iteration procedure. Once the maximum error ϵ converges to the specified tolerance value, the remaining white nodes are computed using the standard Formula (5).

5. Numerical Results

In this section, we present some numerical examples in order to compare the considered iterative methods that were described in the previous sections. For simplicity, the iterative methods were applied to Laplace’s equation,
2 u = 2 u x 2 + 2 u y 2 = 0 , ( x , y ) Ω
and the Dirichlet boundary conditions
u ( x , 0 ) = 0.1 , u ( x , m ) = 1.0 , u ( 0 , y ) = 0.1 , ( m , y ) = 1.0 ,
where 0 x m , 0 y m , and m = N 1 . These conditions are applied to an N × N square rectangle computational grid. The experiments were conducted on a Linux machine equipped with an Intel i5-3570K CPU at 3.40 GHz and 8 GB memory.

5.1. Experiment 1

Initially, small grid sizes were chosen to examine the behaviour of the AM method, where N = 16 , 32 , 48 , 64 , 80 and 96. Different weighted parameter ω values were tested in the range of 1 < ω < 2 . The convergence rate to stop the iteration procedure was 1 . 0 2 ϵ 1 . 0 6 . The number of iterations k and maximum error E were recorded. The boundary conditions (24) were applied and initial inner nodes were set to u ( x * , y * ) = 0.1 , where 1 x * N 2 and 1 y * N 2 .
From Table 1, it can be observed that if the required accuracy is low ϵ = 1 . 0 2 , the iterations count k does not increase even when a larger grid size is used. However, when higher accuracy ϵ = 1 . 0 6 is required, the number of iterations k increases very quickly. It should also be noted that the optimal convergence rate is obtained when the tolerance value is in the range of 1.70 ω 1.80 since, when ω is very close to 2, it will result in a very high iterations count. Therefore, these parameter values were then applied in Experiment 2 so that an optimal convergence rate could be obtained with acceptable accuracy of the solutions.
Figure 6 illustrates that for grid size N = 16 , it requires at least ϵ = 0.01 so that the solutions can be obtained with acceptable accuracy. When larger grid size N = 96 is used, higher accuracy ϵ = 0.0001 is required, as shown in Figure 7.

Stability Analysis

We define that the proposed methods are stable if the obtained maximum absolute error becomes smaller as the computation proceeds. In order to observe this behaviour, the AM and MAM methods were tested with several different parameters. The resulting maximum absolute errors obtained from both methods were observed against the number of iterations. It was assumed that the same behaviour would be achieved for the proposed methods.
Figure 8 demonstrates the stability of the AM and MAM methods with several different parameters. It can be seen clearly from the graphs that as the number of iterations increases, the maximum error becomes smaller to achieve the convergence criterion. This stable condition was true for different grid sizes ( N = { 64 , 128 , 256 , 512 } ), several weighted parameter values ( 1.4 ω 1.7 , 1.2 α r , α b 1.8 ) and various tolerance errors ( 1 . 0 6 ϵ 1 . 0 4 ).

5.2. Experiment 2

In the final experiments, the convergence rate was ϵ = 1.0 × 10 7 , and the chosen values of ω and r were in the range of between 1.70 and 1.90 . Sufficient different grid sizes of N = 128 , 256 , 384 , 512 , 640 and 768 were chosen to record the iteration counts and execution time of the considered methods.
Table 2 shows that the proposed MAAM outperforms the existing AM and MAM methods in terms of iterations and execution time. Moreover, the skewed variants clearly outperform their conventional counterparts, with the SkMAAM method providing the greatest results. The MAAM and SkMAAM methods are also more flexible than the previous AM, MAM, SkAM and SkMAM methods, because they include additional accelerated parameters, allowing for a greater range of parameter values. As shown in the maximum error value, all methods achieve an acceptable tolerance rate in terms of accuracy.
Table 3 shows the reduction percentages in terms of iterations and time. The MAM method performs better than the AM method. In comparison to the AM and MAM methods, the suggested MAAM method decreases iterations and execution time by roughly 33% and 24%, respectively. The skewed methods reduce iterations and execution time by 34% to 44% and 52% to 61%, respectively, when compared to their standard variants. Meanwhile, the improved SkMAAM method has a reduced iteration count and execution time compared to the SkMAAM method by around 14% and 11%, respectively. The SkMAAM is also clearly better than the SkAM method by providing much lower iteration counts and execution time.
Figure 9 shows the comparison between the considered methods in terms of the number of iterations, where the skewed variants (shown in solid lines) require a lower number of iterations compared to their corresponding standard variants (dashed lines). Meanwhile, Figure 10 illustrates the execution time of the considered iterative methods, in which the SkMAAM method has consistently outperformed the other methods. It is also shown that the execution time increases rapidly compared to the size of the grid. The graph shows that the iterations of all methods grow almost linearly with the increase in grid size.

6. Discussion

Only one weighted value is used as a relaxation parameter to speed up the computation for the existing AM and its skewed SkAM variant. For the modified variants of the MAM and MAAM methods and their corresponding skewed methods, two different weighted values are applied for red and black nodes, respectively. Meanwhile, the proposed MAAM and its skewed variant also apply an additional two distinct accelerated parameters. As a result, there is a greater variety of parameter tuning options available. When good parameter values are chosen, as shown in the numerical experiment findings, the MAAM variants execute quicker than the existing AM and MAM variants.
All skewed variants, SkAM, SkMAM and SkMAAM, outperform their corresponding regular variants since they only consider half of the total nodes. Only red and black nodes are computed in the iteration process when using a skewed grid. After the convergence requirement is met, the remaining white nodes are computed only once. The skewed variants lower the number of iterations and execution time by half in terms of computational complexity. In terms of accuracy, all methods achieve an acceptable maximum error rate.
Computation at Level 1 begins by computing all red nodes. After that, all black nodes are computed, in which each black node utilises the updated red nodes. Similarly, at Level 2, the process begins by computing all red nodes and is followed by the computation of all black nodes. Finally, the average of the updated values obtained from Levels 1 and 2 is computed at Level 3. As a result, since these methods require two independent computations that may be performed in parallel, their parallel versions should outperform their sequential implementations. The application of the red–black strategy allows each red node to be computed concurrently and vice versa. Therefore, the computation of red and black nodes can be carried out in alternating order, where the parallel algorithms can be implemented easily. Furthermore, the availability of a fast and affordable Graphical Processing Unit (GPU) makes such parallel implementation very attractive. The suggested improved accelerated two-stage iterative methods with red–black ordering can also be applied to solve other types of partial differential equations.
Although additional parameters provide a wider choice of options, they pose a challenge in finding good optimum values. Since there is no formula currently available to immediately give the optimum value, several trial runs have to be conducted to study the response of the AM method to the chosen parameters, as explained in Experiment 1. It is assumed that the good optimum values found for the AM method are also valid for the other methods since their iterative schemes are all derived from the same finite difference approximation equations.

7. Conclusions

All of the methods investigated in this study have two levels of computations, each of which is totally independent, making them ideal for parallel implementation. As a result, such a parallel solution should greatly reduce the execution time. In addition to the existing relaxation parameter in the AM and SkAM, the use of two weighted parameters provides different relaxation values for red and black nodes applied in the MAM, SkMAM, MAAM and SkMAAM methods. The introduction of two accelerated parameters to the respective red and black nodes further improves the overall performance of the newly proposed MAAM and SkMAAM methods.
The computational complexity of the iterative method implemented on the skewed grid can be decreased by around half by using the half-sweep iteration procedure. Because only half of the total nodes in the computational grid are considered throughout the iteration process, the skewed variant implemented using the half-sweep technique outperforms the regular variant constructed using the full-sweep technique on a regular grid. We plan to look into parallel implementations of the proposed methods in future work. More advanced approaches such as a quarter-sweep technique and implementation on a block grid would be examined. The potential application of the proposed methods for solving fractional-order partial differential equations in the extended work is seriously considered.

Author Contributions

Conceptualisation, A.S.; methodology, A.S.; software, A.S.; validation, A.S. and A.A.D.; formal analysis, A.S.; investigation, A.A.D.; resources, A.A.D.; data curation, A.S.; writing—original draft preparation, A.S. and A.A.D.; writing—review and editing, A.S. and A.A.D.; visualisation, A.S.; supervision, A.S.; project administration, A.S.; funding acquisition, A.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Universiti Malaysia Sabah grant number GTA00240-2019.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
SORSuccessive Overrelaxation
AORAccelerated Overrelaxation
AMArithmetic Mean
MAMModified Arithmetic Mean
MAAMModified Accelerated Arithmetic Mean
SkAMSkewed Arithmetic Mean
SkMAMSkewed Modified Arithmetic Mean
SkMAAMSkewed Modified Accelerated Arithmetic Mean

References

  1. Evans, D.J.; Sahimi, M.S. The alternating group explicit (age) iterative method for solving parabolic equations i: 2-dimensional problems. Int. J. Comput. Math. 1988, 24, 311–341. [Google Scholar] [CrossRef]
  2. Galligani, I.; Ruggiero, V. The arithmetic mean method for solving essentially positive systems on a vector computer. Int. J. Comput. Math. 1990, 32, 113–121. [Google Scholar] [CrossRef]
  3. Muthuvalu, M.S.; Galligani, E.; Ali, M.K.M.; Sulaiman, J.; Lebelo, R.S. Performance analysis of Arithmetic Mean method for solving composite 6-point closed Newton–Cotes quadrature algebraic equation. Asian-Eur. J. Math. 2019, 12, 1950061. [Google Scholar] [CrossRef]
  4. Sahimi, M.S.; Khatim, M. The Reduced Iterative Alternating Decomposition Explicit (RIADE) method for the diffusion equation. Pertanika J. Sci. Tech. 2001, 9, 13–20. [Google Scholar] [CrossRef]
  5. Hussain, A.; Muthuvalu, M.S.; Faye, I.; Ali, M.K.M.; Lebelo, R.S. Numerical Study of Glioma Growth Model with Treatment Using the Two-Stage Gauss-Seidel Method. J. Phys. Conf. Ser. 2018, 1123, 012040. [Google Scholar] [CrossRef]
  6. Sulaiman, J.; Othman, M.; Hasan, M.K. A new Half-Sweep Arithmetic Mean (HSAM) algorithm for two-point boundary value problems. In Proceedings of the International Conference on Statistics and Mathematics and Its Application in the Development of Science and Technology, Bandung, Indonesia, 4–6 October 2004; pp. 169–173. [Google Scholar]
  7. Young, D.M. Iterative Solution of Large Linear Systems; Academic Press: Cambridge, MA, USA, 1971. [Google Scholar]
  8. Saad, Y.; van der Vorst, H.A. Iterative solution of linear systems in the 20th century. J. Comput. Appl. Math. 2000, 123, 1–33. [Google Scholar] [CrossRef] [Green Version]
  9. Hong, E.J.; Saudi, A.; Sulaiman, J. Numerical Evaluation of Quarter-Sweep SOR Iteration for Solving Poisson Image Blending Problem. In Proceedings of the 2018 IEEE International Conference on Artificial Intelligence in Engineering and Technology (IICAIET), Kota Kinabalu, Malaysia, 8 November 2018; pp. 1–5. [Google Scholar] [CrossRef]
  10. Musli, F.A.; Saudi, A. Agent Navigation via Harmonic Potentials with Half-Sweep Kaudd Successive over Relaxation (HSKSOR) Method. In Proceedings of the 2019 IEEE 9th Symposium on Computer Applications Industrial Electronics (ISCAIE), Kota Kinabalu, Malaysia, 27–28 April 2019; pp. 335–339. [Google Scholar] [CrossRef]
  11. Hadjidimos, A. Accelerated overrelaxation method. Math. Comput. 1978, 32, 149–157. [Google Scholar] [CrossRef]
  12. Sunarto, A.; Agarwal, P.; Sulaiman, J.; Chew, J.V.L.; Aruchunan, E. Iterative method for solving one-dimensional fractional mathematical physics model via quarter-sweep and PAOR. Adv. Differ. Equ. 2021, 2021, 147. [Google Scholar] [CrossRef]
  13. Huang, Z.G.; Cui, J.J. Accelerated Relaxation Modulus-Based Matrix Splitting Iteration Method for Linear Complementarity Problems. Bull. Malays. Math. Sci. Soc. 2021, 44, 2175–2213. [Google Scholar] [CrossRef]
  14. Dahalan, A.; Saudi, A.; Sulaiman, J. Autonomous navigation on modified AOR iterative method in static indoor environment. J. Phys. Conf. Ser. 2019, 1366, 012020. [Google Scholar] [CrossRef]
  15. Li, R.; Gong, L.; Xu, M. A heterogeneous parallel Red-Black SOR technique and the numerical study on SIMPLE. J. Supercomput. 2020, 76, 9585–9608. [Google Scholar] [CrossRef]
  16. Ali, U.; Abdullah, F.A. Modified implicit difference method for one-dimensional fractional wave equation. AIP Conf. Proc. 2019, 2184, 060021. [Google Scholar] [CrossRef]
  17. Singh, J.; Ahmadian, A.; Rathore, S.; Kumar, D.; Baleanu, D.; Salimi, M.; Salahshour, S. An efficient computational approach for local fractional Poisson equation in fractal media. Numer. Methods Partial Differ. Equ. 2021, 37, 1439–1448. [Google Scholar] [CrossRef]
  18. Saudi, A.; Sulaiman, J. Half-Sweep Arithmetic Mean Method for Solving 2D Elliptic Equation. Int. J. Appl. Math. Stat. 2019, 58, 52–59. [Google Scholar]
  19. Saudi, A.; Sulaiman, J. An Efficient Two-Stage Half-Sweep Modified Arithmetic Mean (HSMAM) Method for the Solution of 2D Elliptic Equation. Adv. Sci. Lett. 2018, 24, 1417–1921. [Google Scholar] [CrossRef]
  20. Abdullah, A.R. The four point Explicit Decoupled Group (EDG) method: A fast Poisson solver. Int. J. Comput. Math. 1991, 38, 61–70. [Google Scholar] [CrossRef]
  21. Ruggiero, V.; Galligani, E. An iterative method for large sparse linear systems on a vector computer. Comput. Math. Appl. 1990, 20, 25–28. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Portion of regular computational grid about point ( i , j ) for AM method at (a) Level 1 and (b) Level 2.
Figure 1. Portion of regular computational grid about point ( i , j ) for AM method at (a) Level 1 and (b) Level 2.
Symmetry 14 00993 g001
Figure 2. Portion of skewed computational grid about point ( i , j ) for SkAM method at (a) Level 1 and (b) Level 2.
Figure 2. Portion of skewed computational grid about point ( i , j ) for SkAM method at (a) Level 1 and (b) Level 2.
Symmetry 14 00993 g002
Figure 3. The standard computational mesh of (a) regular and (b) skewed cases for the AM and SkAM methods, respectively.
Figure 3. The standard computational mesh of (a) regular and (b) skewed cases for the AM and SkAM methods, respectively.
Symmetry 14 00993 g003
Figure 4. The red–black computational mesh for (a) regular and (b) skewed cases.
Figure 4. The red–black computational mesh for (a) regular and (b) skewed cases.
Symmetry 14 00993 g004
Figure 5. The computational molecules for finite difference approximation of regular (left a,b) and skewed (right c,d) cases. The molecules of red and black nodes are symmetrical for both cases, in which the black node applies the updated values of its four neighbouring red nodes and vice versa.
Figure 5. The computational molecules for finite difference approximation of regular (left a,b) and skewed (right c,d) cases. The molecules of red and black nodes are symmetrical for both cases, in which the black node applies the updated values of its four neighbouring red nodes and vice versa.
Symmetry 14 00993 g005
Figure 6. Solution of Laplace’s equation using different maximum error rate and grid size N = 16 . Graphs (a,b) give lower accuracy, whereas (c,d) illustrate the solutions that require ϵ = 0.01 toobtain acceptable accuracy.
Figure 6. Solution of Laplace’s equation using different maximum error rate and grid size N = 16 . Graphs (a,b) give lower accuracy, whereas (c,d) illustrate the solutions that require ϵ = 0.01 toobtain acceptable accuracy.
Symmetry 14 00993 g006
Figure 7. Solution of Laplace’s equation for grid size N = 96 with different maximum error rate. Acceptable accuracy is obtained when ϵ = 1 . 0 4 . (a) ϵ = 0.01 ; (b) ϵ = 0.001 ; (c) ϵ = 0.0001 ; (d) ϵ = 0.00001 .
Figure 7. Solution of Laplace’s equation for grid size N = 96 with different maximum error rate. Acceptable accuracy is obtained when ϵ = 1 . 0 4 . (a) ϵ = 0.01 ; (b) ϵ = 0.001 ; (c) ϵ = 0.0001 ; (d) ϵ = 0.00001 .
Symmetry 14 00993 g007
Figure 8. As the computation proceeds, the maximum errors become smaller. (a) N = 64 , ω = 1.4 , ϵ = 1 . 0 4 ; (b) N = 64 , α r = 1.2 , α b = 1.4 , ϵ = 1 . 0 5 ; (c) N = 128 , ω = 1.5 , ϵ = 1 . 0 5 ; (d) N = 128 , α 5 = 1.5 , α b = 1.6 , ϵ = 1 . 0 5 ; (e) N = 256 , ω = 1.5 , ϵ = 1 . 0 5 ; (f) N = 256 , α 5 = 1.5 , α b = 1.6 , ϵ = 1 . 0 5 ; (g) N = 512 , ω = 1.7 , ϵ = 1 . 0 6 ; (h) N = 512 , α r = 1.7 , α b = 1.8 , ϵ = 1 . 0 6 .
Figure 8. As the computation proceeds, the maximum errors become smaller. (a) N = 64 , ω = 1.4 , ϵ = 1 . 0 4 ; (b) N = 64 , α r = 1.2 , α b = 1.4 , ϵ = 1 . 0 5 ; (c) N = 128 , ω = 1.5 , ϵ = 1 . 0 5 ; (d) N = 128 , α 5 = 1.5 , α b = 1.6 , ϵ = 1 . 0 5 ; (e) N = 256 , ω = 1.5 , ϵ = 1 . 0 5 ; (f) N = 256 , α 5 = 1.5 , α b = 1.6 , ϵ = 1 . 0 5 ; (g) N = 512 , ω = 1.7 , ϵ = 1 . 0 6 ; (h) N = 512 , α r = 1.7 , α b = 1.8 , ϵ = 1 . 0 6 .
Symmetry 14 00993 g008
Figure 9. The number of iterations for the considered iterative methods.
Figure 9. The number of iterations for the considered iterative methods.
Symmetry 14 00993 g009
Figure 10. The execution time of the considered iterative methods.
Figure 10. The execution time of the considered iterative methods.
Symmetry 14 00993 g010
Table 1. Performance of AM method in terms of iterations count (k) and maximum error (E). N is the grid size and ω denotes the weighted parameter.
Table 1. Performance of AM method in terms of iterations count (k) and maximum error (E). N is the grid size and ω denotes the weighted parameter.
ω 1.201.701.801.99
ϵ 1.0 2 1.0 4 1.0 6 1.0 6
N k E k E k E k E
16 24 9.21 × 10 4 34 8.69 × 10 6 55 8.84 × 10 8 340 9.64 × 10 8
32 34 9.87 × 10 4 89 9.77 × 10 6 111 9.42 × 10 8 583 9.77 × 10 8
48 35 9.78 × 10 4 164 9.99 × 10 6 210 9.70 × 10 8 679 9.92 × 10 8
64 35 9.78 × 10 4 256 9.84 × 10 6 336 9.89 × 10 8 712 9.93 × 10 8
80 35 9.78 × 10 4 359 9.98 × 10 6 489 9.95 × 10 8 769 9.66 × 10 8
96 35 9.78 × 10 4 473 9.93 × 10 6 668 9.90 × 10 8 768 9.67 × 10 8
Table 2. Number of iterations (k) and execution time (t) in seconds of the considered iterative methods. N is the size of grid, α 1 and α 2 are the weighted parameters, β 1 and β 2 are the accelerated parameters, and E is the maximum error.
Table 2. Number of iterations (k) and execution time (t) in seconds of the considered iterative methods. N is the size of grid, α 1 and α 2 are the weighted parameters, β 1 and β 2 are the accelerated parameters, and E is the maximum error.
MethodsN α 1 α 2 β 1 β 2 ktE
REGULAR MESHAM1281.701.701.701.7016410.83 9.9836 × 10 8
2561.711.711.711.71531211.05 9.9880 × 10 8
3841.721.721.721.7210,22846.55 9.9965 × 10 8
5121.731.731.731.7315,961135.95 9.9987 × 10 8
6401.741.741.741.7422,214292.58 9.9986 × 10 8
7681.751.751.751.7528,763570.83 9.9999 × 10 8
MAM1281.721.701.721.7214530.73 9.9659 × 10 8
2561.731.711.731.7347019.41 9.9999 × 10 8
3841.741.721.741.74903740.99 9.9947 × 10 8
5121.751.731.751.7514,067122.14 9.9988 × 10 8
6401.761.741.761.7619,519255.97 9.9994 × 10 8
7681.771.751.771.7725,186492.20 9.9993 × 10 8
MAAM1281.721.701.741.8011160.55 9.9433 × 10 8
2561.731.711.751.8136047.56 9.9814 × 10 8
3841.741.721.761.82687731.47 9.9943 × 10 8
5121.751.731.771.8310,60691.70 9.9967 × 10 8
6401.761.741.781.8414,554188.85 9.9993 × 10 8
7681.771.751.791.8518,539363.39 9.9977 × 10 8
SKEWED MESHSkAM1281.701.701.701.708850.32 9.9756 × 10 8
2561.711.711.711.7129134.21 9.9991 × 10 8
3841.721.721.721.72567418.54 9.9892 × 10 8
5121.731.731.731.73893653.10 9.9978 × 10 8
6401.741.741.741.7412,538117.61 9.9961 × 10 8
7681.751.751.751.7516,352222.95 9.9992 × 10 8
SkMAM1281.721.701.721.728000.31 9.9445 × 10 8
2561.731.711.731.7326343.81 9.9707 × 10 8
3841.741.721.741.74512016.68 9.9929 × 10 8
5121.751.731.751.75804648.31 9.9943 × 10 8
6401.761.741.761.7611,258106.15 9.9955 × 10 8
7681.771.751.771.7714,637210.31 9.9968 × 10 8
SkMAAM1281.721.701.741.806910.27 9.9883 × 10 8
2561.731.711.751.8122763.56 9.9751 × 10 8
3841.741.721.761.82440915.60 9.9919 × 10 8
5121.751.731.771.83689642.77 9.9993 × 10 8
6401.761.741.781.84959795.11 9.9937 × 10 8
7681.771.751.791.8512,398168.90 9.9991 × 10 8
Table 3. Average reduction percentage in terms of iterations count and execution time.
Table 3. Average reduction percentage in terms of iterations count and execution time.
METHODSIterations ReductionTime Reduction
MAM vs. AM11.84%12.55%
MAAM vs. AM33.41%33.68%
MAAM vs. MAM24.48%24.14%
SkAM vs. AM44.41%60.87%
SkMAM vs. MAM43.21%58.77%
SkMAAM vs. MAAM35.50%51.79%
SkMAM vs. SkAM9.93%7.85%
SkMAAM vs. SkAM22.76%18.29%
SkMAAM vs. SkMAM14.24%11.25%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Saudi, A.; Dahalan, A.A. An Efficient Red–Black Skewed Modified Accelerated Arithmetic Mean Iterative Method for Solving Two-Dimensional Poisson Equation. Symmetry 2022, 14, 993. https://doi.org/10.3390/sym14050993

AMA Style

Saudi A, Dahalan AA. An Efficient Red–Black Skewed Modified Accelerated Arithmetic Mean Iterative Method for Solving Two-Dimensional Poisson Equation. Symmetry. 2022; 14(5):993. https://doi.org/10.3390/sym14050993

Chicago/Turabian Style

Saudi, Azali, and A’qilah Ahmad Dahalan. 2022. "An Efficient Red–Black Skewed Modified Accelerated Arithmetic Mean Iterative Method for Solving Two-Dimensional Poisson Equation" Symmetry 14, no. 5: 993. https://doi.org/10.3390/sym14050993

APA Style

Saudi, A., & Dahalan, A. A. (2022). An Efficient Red–Black Skewed Modified Accelerated Arithmetic Mean Iterative Method for Solving Two-Dimensional Poisson Equation. Symmetry, 14(5), 993. https://doi.org/10.3390/sym14050993

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop