Algorithms for Solving the Parametric Self-Adjoint 2 D Elliptic Boundary-Value Problem Using High-Accuracy Finite Element Method

Algorithms for Solving the Parametric Self-Adjoint 2D Elliptic Boundary-Value Problem Using High-Accuracy Finite Element Method A. A. Gusev*, O. Chuluunbaatar*†, S. I. Vinitsky*‡, V. L. DerbovS, A. Góźdź¶ * Joint Institute for Nuclear Research 6, Joliot-Curie, Dubna, Moscow region, Russia, 141980 † Institute of Mathematics, National University of Mongolia, Ulaanbaatar, Mongolia ‡ RUDN University (Peoples’ Friendship University of Russia) 6, Miklukho-Maklaya str., Moscow, Russia, 117198 S Saratov State University, Saratov, Russia ¶ Institute of Physics, University of M. Curie-Sk lodowska, Lublin, Poland


Introduction
The adiabatic representation is widely applied for solving multichannel scattering and bound-state problems for systems of several quantum particles in molecular, atomic and nuclear physics [1][2][3][4].
Such problems are described by elliptic boundary value problems (BVPs) in a multidimensional domain of the configuration space, solved using the Kantorovich method, i.e., the reduction to a system of self-adjoint ordinary differential equations(SODEs) using the basis of surface functions of an auxiliary BVP depending on the independent variable of the SODEs parametrically.The elements of matrices of variable coefficients of these SODEs including the matrix of the first derivatives are determined by the integrals of products of surface eigenfunctions and/or their first derivatives with respect to the parameter.Thus, the key problem of such method is to develop effective algorithms and programs for calculating with given accuracy the surface eigenfunctions and the corresponding eigenvalues of the auxiliary BVP, together with their derivatives with respect to the parameter, and the corresponding integrals that present the matrix elements of the effective potentials in the SODEs [5][6][7].
In this paper we propose new calculation schemes and algorithms for the solution of the parametric 2D elliptic boundary-value problem using high-accuracy finite element method (FEM) with rectangular and triangular elements.The algorithms were implemented in a package of programs that calculate with the given accuracy eigenvalues, eigenfunctions and their first derivatives with respect to the parameter of the parametric self-adjoint elliptic differential equations with the boundary conditions of the Dirichlet and/or Neumann type in the 2D finite domain and the integrals of products of the surface eigenfunctions and their first derivatives with respect to the parameter that express the matrix elements of the effective potentials in the SODEs.
We also propose a method of constructing the etalon potential in the auxiliary parametric BVP that allows the calculation of the parametric surface basis functions in the analytical form.These functions can be then used for the reduction of the original 2D BVP to the SODEs containing the additional potential matrix, representing the discrepancy between the original potential and the etalon one, averaged with the basis functions.The efficiency of the calculation schemes and algorithms is demonstrated by benchmark calculation of the 2D BVPs describing quadrupole vibrations in the collective nuclear model [4,8,9].
The structure of the paper is the following.In Section 2 the Kantorovich method for solving the 2D and 3D BVPs is considered.In Sections 3 and 4 the 2D FEM schemes and algorithms for solving the parametric 2D BVP and calculating derivatives with respect to the parameter together with the corresponding matrix elements are presented.In Section 5 the benchmark calculations of 2D FEM algorithms and programs are analyzed.In Conclusion we discuss the results and perspectives.
During the simulation the etalon potential  0 (  ;   , (  )) in Eq. ( 5) will be chosen as  0 (  ;   , (  )) =  (  ,   ) or calculated separately for different values of   ∈ Ω   from the conditions min If this parametric eigenvalue problem has no analytical solution, then it is solved numerically by the FEM: at  − 1 = 1 using the program ODPEVP [6] and at  − 1 = 2 using the program, implementing the algorithm presented in Section 3. Substituting the expansion (3) into Eq.( 1) with Eqs. ( 6) and ( 7) taken into account, we arrive at the set of self-adjoint ODEs for the unknown vector functions  () Here I, W(  ) and Q(  ) are the  max ×  max matrices The effective potentials   (  )=  (  ) are calculated by integrating the difference  (  ,   ) =  (  ,   )−  (  ;   , (  )) with the basis functions: If the etalon potential is equal to the original one,   (  ;   , (  ))= (  ,   ) then   (  )=0.If the original potential  (  ,   ) is presented in the tabular form, e.g., as in Ref. [4], then the etalon potential   (  ;   , (  )) can be considered as a certain approximation or interpolation of the original one, and the approximation error  (  ,   ) ̸ = 0.In this case the etalon potential   (  ;   , (  )) can be determined minimizing  (  ,   ) in accordance with the condition (8) of the least squares method.Then the integrals   (  ) of the approximation error  (  ,   ) averaged with the parametric basis functions are included into the final system of ODEs.It means that we take the approximation error of the original potential into account and require the solutions of the 2D BVP (1) to have the given accuracy, see, e.g., Section 5.4.
The solutions of the discrete spectrum  : for any  and for any .
The functions {   ()}  =0 and {   ()}  =0 form a basis in the space of polynomials of the -th order.Now, the function Φ(, ; ) ∈ ℱ ℎ  ∼ ℋ 1 (Ω ℎ  ,ℎ  ) is approximated by a finite sum of piecewise polynomial functions    () and    () 2. The domain Ω(, ) = ⋃︀  =1 ∆  , specified as a polygon in the plane ( ≡  1 ,  ≡  2 ) ∈ ℛ 2 , is covered with finite elements, the triangles ∆  with the vertices ( 11 ,  21 ), ( 12 ,  22 ), ( 13 ,  23 ) (here   ≡  ; ,  = 1, 2,  = 1, 3,  = 1, ).On each of the triangles ∆  (the boundary is considered to belong to the triangle) the shape functions    ( 1 ,  2 ) are introduced.For this purpose we divide the sides of the triangle into  equal parts and draw three families of parallel straight lines through the partition points.The straight lines of each family are numbered from 0 to , so that the line passing through the side of the triangle has the number 0, and the line passing through the opposite vertex of the triangle has the number .
Three straight lines from different families intersect in one point   ∈ ∆  , which will be numbered by the triplet ( 1 ,  2 ,  3 ),   ≥ 0,  1 +  2 +  3 = , where  1 ,  2 and  3 are the numbers of the straight lines passing parallel to the side of the triangle that does not contain the vertex ( 11 ,  21 ), ( 12 ,  22 ) and ( 13 ,  23 ), respectively.The coordinates of this point   = ( 1 ,  2 ) are determined by the expression As shape functions we use the Lagrange triangular polynomials    ( 1 ,  2 ) of the order  that satisfy the condition , equal 1 in one of the points   and zero in the other points.
In this method the piecewise polynomial functions    ( 1 ,  2 ) in the domain Ω are constructed by joining the shape functions    ( 1 ,  2 ) in the triangle ∆  : The vector function has a generalized first-order partial derivative and belongs to the Sobolev space ℋ 1 (Ω ℎ  ,ℎ  ) [10].After substituting the expansion (19) or (20) into the variational functional (18) and minimizing it [10,11], we obtain the generalized eigenvalue problem Here A  is the stiffness matrix; B  is the positive definite mass matrix;  ℎ is the vector approximating the solution on the finite-element grid; and  ℎ is the corresponding eigenvalue.The matrices A  and B  have the following form: where the local matrices a   and b   are calculated for rectangular elements as , , , ,  = 0, , or the matrix elements    ′ and    ′ are calculated for triangular elements as where ∆  is a triangular element with the vertices ( 11 ,  21 )( 12 ,  22 )( 13 ,  23 ), ∆ is the triangle with vertices (0, 0)(1, 0)(0, 1), and In this case we have explicit expression for shape functions    ( ′ 1 ,  ′ 2 ): We seek the relation between the global and the local coordinates of the triangle vertices in the form Substituting the coordinates of three vertices in the global and the local frame into Eq.(26), we obtain the system of algebraic equations for calculating the coefficients  ** :  23)-( 24) are evaluated using the Gaussian quadrature of the order  + 1.Note that in the present approach the maximum half-band width of the matrices A  and B  is small compared to their dimension and is not greater than ( + 1)( + 1).
In order to solve the generalized eigenvalue problem (21), the subspace iteration method [10,11] elaborated by Bathe [11] for the solution of large symmetric bandedmatrix eigenvalue problems has been chosen.This method uses the skyline storage mode which stores the components of the matrix column vectors within the banded region of the matrix, and is ideally suited for banded finite-element matrices.The procedure chooses a vector subspace of the full solution space and iterates upon the successive solutions in the subspace (for details, see [11]).The iterations continue until the desired set of solutions in the iteration subspace converges to within the specified tolerance on the Rayleigh quotients for the eigenpairs.If the matrix A  in Eq. ( 21) is not positively defined, the problem (21) is replaced with the following problem: The number  (the shift of the energy spectrum) is chosen such that the matrix Ã is positive defined.The eigenvector of the problem (27) is the same, and  ℎ = εℎ + .

The algorithm for calculating the parametric derivatives of eigenfunctions and the matrices of effective potentials
Taking a derivative of the boundary problem ( 15)-( 17) with respect to the parameter , we find that   Φ  (, ; ) is a solution of the following boundary problem Below we present an efficient numerical method that allows the calculation of   Φ  (, ; ) with the same accuracy as achieved for the eigenfunctions of the BVP ( 15)- (17) and the use of it for computing the matrices of the effective potentials defined as The boundary problem (28)-( 30) is reduced to the linear system of inhomogeneous algebraic equations with respect to the unknown  ℎ /: The normalization condition (17), the condition of orthogonality between the function and its parametric derivative (30), and the additional conditions (29) for the solution of (32) read as Then the potential matrix elements  ℎ  () and  ℎ  () (31) corresponding to Eqs. ( 11)-( 12) can be calculated using the formulas Since  ℎ is an eigenvalue of (21), the matrix L in Eq. ( 32) is degenerate.In this case the algorithm for solving Eq. ( 32) can be written in three steps as follows: Step k1.Calculate the solutions v and w of the auxiliary inhomogeneous systems of algebraic equations with the non-degenerate matrix L and the right-hand sides b and d where  is the number of the element of the vector B   ℎ having the greatest absolute value.
Step k2.Evaluate the coefficient Step k3.Evaluate the vector From the above consideration it is evident, that the computed derivative has the same accuracy as the calculated eigenfunction.
The following theorem can be formulated.
Theorem 1.Let (, ; ) in Eq. ( 15) be a continuous and bounded positively defined operator on the space ℋ 1 with the energy norm.Also let    (, ; ) be continuous and bounded for each value of the parameter .Then for the exact values of the solutions     (),   Φ  (, ; ) ∈ ℋ 2 ,   (),   () from ( 28)-( 31) and the corresponding numerical values 32)-(34) the following estimates are valid: where ℎ is the largest distance between any two points of the finite element ∆  ,  is the order of finite elements, ,  are the numbers of the corresponding solutions, and the constants  3 ,  4 ,  5 and  6 are independent of the step ℎ.
The proof is straightforward following the proof schemes in accordance with [10,13].
The solution of 2D BVP in the triangular domain As a benchmark example we consider the exactly solvable BVP for a membrane in the form of equilateral triangle with side equal to 4/3, in the conventional variables (, ) ∈ Ω(, ) with the Dirichlet or Neumann conditions at the boundary Ω(, ) of the region Ω(, ).
The BVPs were solved on the uniform finite-element grid composed of  2 equilateral triangles with the side equal to 4/(3).In Table 1 we present the estimations of the finite-element scheme order  depending on the size of elements.
The parametric surface eigenfunctions and their derivatives with respect to the parameter under the dipole splitting are shown in Fig. 3.
We choose the mass parameter to be = 2 =124.Thus, there are ground and doubly degenerate excited states, localized in three wells.
The point symmetry group C 3 of the problem (1), (43) has four irreducible representations (irrs.)A1, A2, E1 and E2 to classify the solutions, the E-type states being doubly degenerate [17,18].The calculated eigenvalues are presented in Table 2.The first eigenfunctions for each of the irrs.A1, E1, E2 are shown in Fig. 4.

Parametric surface functions for KM in the analytical form
Let us consider the BVP for Eq. ( 5) with the etalon potential   (  ;   , (  )): where (  ) is the set of parameters, (  ) = { 0 (  ),  2 (  ),  0 (  )}.In the considered case the parametric eigenvalue problem ( 5)-( 7) has an exact solution, i.e., the parametric eigenfunctions Φ  (  ;   ) and potential curves   (  ) are expressed in the analytical form The integration in the effective potentials (11) with the basis functions (45) is carried out analytically, which yields the expressions ( For the potential (43) from the condition (46) we have We performed the calculations for these parameters in the case when (  ,   )=( 22 ,  20 ), as well as (  ,   )=( 20 ,  22 ).The results coincide with those of the calculations of  ℎ  performed in the previous Section 5.3 and presented in Table 2 with 9 significant digits.

Conclusion
We elaborated new calculation schemes and algorithms for solving the parametric 2D elliptic BVP using the high-accuracy FEM with rectangular (22) and triangular (23)-(24) elements.The algorithm and the programs calculate with the given accuracy the eigenvalues, the eigenfunctions and their first derivatives with respect to the parameter of the parametric self-adjoint elliptic differential equations with boundary conditions of the Dirichlet and/or Neumann type in the finite 2D domain ( 15)-( 17), (18) and the corresponding inhomogeneous boundary-value problem (28)-(33), obtained by taking a parametric derivative of the original 2D BVP.The program also calculates the potential matrix elements, the integrals of the eigenfunctions multiplied by their first derivatives with respect to the parameter (34).The parametric eigenvalues (potential curves) and the matrix elements computed by the program can be used for solving the bound-state and multi-channel scattering problems for a system of the coupled second-order ODES with using the Kantorovich method.We demonstrated the efficiency of the proposed calculation schemes, algorithms and codes by the example of solving the boundary-value problem of a quadruple vibration collective nuclear model.
Table 1 Discrepancies  () = | ℎ  () −   | of the first  1 = 3 and the fourth  4 = 3 eigenvalues of the BVP(40) with Dirichlet and Neumann boundary conditions, respectively, and the Runge coefficients Ru  = log 2 ((  () −   (2))/(  (2) −   (4))) for the schemes of the first  = 1 and the second  = 2 orders;  is the number of elements to correspond rigorously to the theoretical estimations (38) of the order (ℎ 2 ) of the eigenvalues approximation.The first eight eigenfunctions of the BVPs (40) with the Neumann and Dirichlet boundary conditions are shown in Figs. 1 and 2.

Figure 3 .
Figure 3.The first four eigenfunctions of the parametric 2D BVP with the potential function (41) and their derivatives with respect to the parameter at  = 0.