Investigations of Linear System Obtained by Element-Free Galerkin
Method – Influence of Imposing Method for Boundary Condition –Soichiro Ikuno1 , Yusuke Oishi1 , Taku Itoh2 , and Atsushi Kamitani3
1

Schoolof Computer Science, Tokyo University of Technology, Japan
2
Faculty of Science and Technology, Seikei University, Japan3
Department of Informatics, Faculty of Engineering, Yamagata University, Japan

1Introduction

A mesh generation procedure must be needed before
using general finite element method (FEM) for solution of
the partial diﬀerential equations. However, it costs a lotof time to divide the region into a set of finite elements.
Thus, the time required for discretizing the partial diﬀerential equation or solving the resulting linear system is much
shorter than that for the element generation. On the other
hand, meshless approaches do not require finite elements ofa geometrical structure. The necessary information is only
locations of nodes which are scattered in the region andon the boundary. For these reasons, various meshless approaches have been developed, such as the diﬀuse element
method, the Element-Free Galerkin (EFG) method [1] andthe Meshless Local Petrov-Galerkin (MLPG) method [2].
The purpose of the present study is to implement thevariable preconditioned Minimal Residual method withgeneralized Jacobi and to solve the linear system obtainedby EFG with two types of imposing method for the boundary condition.
2Element-Free Galerkin Method

In this section, we derive the weak form of the Poissonproblem in the analytic two-dimensional region Ω and discretize it by means of the Element-Free Galerkin (EFG)
method. For the simplicity, we assume that the regionΩ is bounded by essential boundary ΓD . The governing
equation is expressed as−∆u = f,(1)where u(x, y) denotes unknown function and f (x, y) denotes the given functions. Furthermore, the boundary conditions on ΓD is assumed as u = ū. Here, ū denotes thegiven function.
EFG is based on the global Galerkin formulation, andthe boundary conditions are satisfied by the Lagrange multiplier. That is to say, the governing equation and boundary conditions can be proved equivalent to the following
weak form:ZZpδu dΩ∇u · ∇δu dΩ −ΩΩZZδλ(u − ū) d` = 0.
(2)
λδu d` ++ΓΓThe function δu and δλ are arbitrary functions and wecall these functions test functions. Besides, u is the trial
function and λ denotes the Lagrange multiplier. Here,when the trial and test functions are taken from the samefunctional space such as u, δu ∈ span(φ1 , φ2 , · · · , φM ),λ, δλ ∈ span(N1 , N2 , · · · , Nn ), the weak form is discretizedas»–» – » –W B uf=.
(3)
BT O λgwhere φj denotes the shape function of EFG for j-th nodeand Nj denotes the shape function for j-th boundary node.
In the present study, the delta function and the first order
B-spline is adopted for the functions Nj . Furthermore, thematrices A, B and vectors f , g are defined as follows.
ZZ(W )ij =∇φi · ∇φj dΩ, (B)ij =φi Nj d`,ΩΓZZ(f )i =φi p dΩ, (g)j =Nj ū d`.
ΩΓHere ( )ij represents the (i, j)-th matrix element and ( )irepresents the i-th vector component.
We can see from (3), the left lower block matrix of thecoeﬃcient matrix becomes zero matrix. From this reasons,the stationary iterative method, such as Jacobi method
and SOR method can not be applied for the solver.

3VariableMethod

PreconditionedKrylov

Subspace

It is well known that a preconditioning strategy can improve the performance for solving a linear system Ax = b
using the Krylov subspace method and various preconditioning strategies have been developed and numericallyinvestigated. Here, A, x and b denote a coeﬃcient matrix, an unknown vector, and a known vector, respectively.
Generally, a preconditioned matrix M is determined byincomplete LU decomposition, and a vector M −1 rk is calculated at k-th iteration by using backward substitution
or incomplete Cholesky factorization. Here, rk denotes
residual vector at k-th iteration. The calculation time ofsolving linear system is relatively large for the preconditioned part.
K. Abe et al. developed new preconditioning strategy
which is called the variable preconditioning method [3].
Variable Preconditioned (VP) GCR has two nested iterations for GCR and variable preconditioning for GCR arecalled as outer-loop and inner-loop, respectively. In thepreconditioned procedure, the residual equation is solved
to determine the preconditioner for the outer-loop. The algorithm of VPGCR can be extended for symmetric matrix
solver by using other Krylov subspace method. The algorithm of variable preconditioned Minimal Residual (VP-Let x0 be an initial guess.
c1 = −1, c0 = 1, s1 = 0, s0 = 0
q0 = 0, v0 = 0, v−1 = 0
r0 = b − Ax0
Roughly√ solve Az0 = r0
β0 = (r0 , z0 ), η1 = β0
q1 = r0 /β0 , z1 = z0 /β0
for k = 1, 2, · · · , until |ηk+1 |/β0 ≤ ε dork = Azk − βk−1 qk−1
αk = (rk , zk )rk+1 = rk − αk qk
Roughly√ solve Azk+1 = rk+1
βk = (rk+1 , zk+1 )γ̂k = −c√ k αk − ck−1 sk βk−1
γk = γ̂ 2 + βk2
δk = sk αk − ck ck−1 βk−1
k = sk−1 βk−1 , ck+1 = γ̂k /γk , sk+1 = βk /γkvk+1 = (zk − k vk−1 − δk vk )/γkxk+1 = xk + ck+1 ηk vk+1
qk+1 = rk+1 /βk , zk+1 = zk+1 /βkηk+1 = sk+1 ηk
end forFig. 2. The residual history of VPMinRes. The boundary conditions are implemented by using the first order
B-spline. A: case with (5), B: case with (6), C: case with(7)Fig. 1. The algorithm of variable preconditioned Minimal
Residual (VPMinRes) method.

MinRes) is shown in Fig. 1.
In the inner-loop procedure of VP Krylov subspacemethod, the residual equation can be solved roughly byusing some iterative method with only a few iteration, anda stationary iterative method such as Jacobi method andSOR method is adopted for the solver. As we mentioned
above, however, the stationary iterative method can not
be applied for the linear system (3). From this reason, theGeneralized Jacobi method is adopted for the solver for
inner-loop [4].
To adopt Generalized Jacobi method for the inner-loop
solver, we transform the coeﬃcient matrix A of the linearsystem (3) asA = D + L + U = (D + T ) + L + (U − T ),(4)

where D, L and U denote a diagonal, lower and upper
matrix of A. In addition, the elements of the diagonal
matrix T = diag(t1 , t2 , . . . , tn ) are selected as follows.
ti = sgn(aii ) × max |aij |(5)ti = sgn(aii ) × 1.0
X
ti = sgn(aii ) ×|aij |(6)

jFig. 3. The residual history of VPMinRes. The boundary
conditions are implemented by using the Delta function.
A: case with (5), B: case with (6), C: case with (7)We can see from these figure that if the first order Bspline function is adopted for imposing method, it takes
much iteration to solve the system than that of Delta function. This result is derived from the condition number ofthe coeﬃcient matrix. Actually, the condition number ofthe coeﬃcient matrix A with first order B-spline and theDelta function are 1.9684 × 1016 and 352.10. The supposable cause is an accuracy of numerical integration of theB-spline.
References(7)[1] T. Belytschko, Y. Y. Lu and L. Gu, “Element-freeGalerkin methods,” Int. J. Numer. Methods Eng., vol.
37, pp. 229-256, 1994.

Here, aij denotes a (i, j)-th matrix element of the coeﬃcient matrix A.

[2] S. N. Atluri and T. Zhu, “A new meshless local PetrovGalerkin (MLPG) approach in computational mechanics,”
Comput. Mech., vol. 22, pp. 117-127, 1998.

4[3] K. Abe, S. L. Zhang, “A variable preconditioning using
the SOR method for GCR-like methods”, Int. J. Numer.
Anal. Model. 2, no. 2, (2005) 147-161.

i6=jNumerical Results and Conclusion

As we mentioned above, two types of imposing method
for boundary condition are adopted. One is the delta function and the other is first order B-spline for the shape functions Nj . Moreover, VPMinRes with generalized Jacobi isadopted for the solver of the linear system (3), and theresidual histories of the system are shown in Fig. 2 andFig. 3.

[4] A. Hadjidimos, “On the Optimization of the Classical Iterative Schemes for the Solution of Complex Singular LinearSystem”, SIAM. J. on Algebraic and Discrete Methods 6,(1983) 555-566.