Algorithms for solving the parametric self-adjoint 2D elliptic boundary-value problem using high-accuracy finite element method Текст научной статьи по специальности «Физика»

UDC 517.958, 530.145.6, 519.632.4 DOI: 10.22363/2312-9735-2017-25-1-36-55

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. Gozdz^

* 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-Sklodowska, Lublin, Poland

We consider the calculation schemes for solving elliptic boundary-value problems (BVPs) within the framework of the Kantorovich method that provides the reduction of an elliptic BVP to a system of coupled second-order ordinary differential equations (ODEs). The surface basis functions of the expansion depend on the independent variable of the ODEs parametrically. Here we use the basis functions calculated by means of the finite element method(FEM), as well as the probe parametric surface basis functions calculated in the analytical form.

We propose new calculation schemes and algorithms for solving the parametric self-adjoint elliptic boundary-value problem (BVP) in a 2D finite domain, using high-accuracy finite element method (FEM) with rectangular and triangular elements. The algorithm and the programs calculate with the given accuracy the eigenvalues, the surface eigenfunctions and their first derivatives with respect to the parameter of the BVP for parametric self-adjoint elliptic differential equation with the Dirichlet and/or Neumann type boundary conditions on the 2D finite domain, and the potential matrix elements, expressed as integrals of the products of surface eigenfunctions and/or their first derivatives with respect to the parameter. The parametric eigenvalues (potential curves) and the potential matrix elements computed by the program can be used for solving bound-state and multi-channel scattering problems for systems of coupled second-order ODEs by means of the Kantorovich method.

We demonstrate the efficiency of the proposed calculation schemes and algorithms in benchmark calculations of 2D elliptic BVPs describing quadrupole vibrations of a collective nuclear model.

Key words and phrases: parametric elliptic boundary-value problem, finite element method, Kantorovich method, systems of second-order ordinary differential equations

1. 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-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

Received 9th January, 2017.

The authors thank Artur Dobrowolski for collaboration.

This work was supported by the Polish-French COPIN collaboration of the project 04-113, the Bogoliubov-Infeld JINR program and the Russian Foundation for Basic Research (grants 17-01-00298, 17-01-00785).

The reported study was funded within the Agreement N 02.a03.21.0008 dated 24.04.2016 between the Ministry of Education and Science of the Russian Federation and RUDN University.

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-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.

2. Kantorovich method with the etalon potential

Let us consider the BVP in the domain Q(xf ,xs) C Rn_ 1 x R1:

(-¿lb é è+flfr+F ■x-)-E) *<" ■I-)=0- (1)

where D(xf; xs) is a self-adjoint elliptic differential operator in the finite region fi(xf ; xs) C Rn_ 1, E is the spectral parameter, corresponding to the energy of the quantum system, fsi(xs) > 0, dXsfsi(xs) and V(xf,xs) dXsV(xf,xs) are real-valued continuous bounded functions in fi(xf ,xs), and tf(xf ,xs) satisfies the Dirichlet and/or Neumann boundary condition (BC) at the boundary dfi = dfi(xf ,xs) of the domain fi(xf ,xs) and the orthonormalization conditions

{^il^j ) = fsi(xs)^i(xf ,Xs)Vj (xf ,Xs)dxnf1 dxs=5ij. (2) jQ(xf ,xs)

The solution tf(xf , £s)eW|(fi) of the BVP (1) is sought in the form of the Kan-torovich expansion [3]

Jmax

tfi(Xf ,Xs)=YI (xf ; Xs)Xji(Xs), (3)

3=1

using the set of parametric eigenfunctions $ j(xf; xs) G T(xs) ~ W22(Q(xf; xs)) of the parametric BVP in the domain Q(xf; xs) C R"-1

(D(Xf; Xs)-£j(xs)) $(xf; ®s)=0. (4)

For example, D(xf; xs) at n — 1 = 2 is determined in a conventional form in Section 3. To avoid cumbersome notations, below we will use the simplest definition of D(xf; xs) at n — 1 = 1

&2

D(xf; Xs) = — +V0(xf; Xs,g(xs)), (5)

where V0(xf; xs,g(xs)) is the etalon potential defined in the interval Xf G (xmin(xs), xmax(xs))=QXf (xs) and depending on the variable xs G QXs as a parameter. We assume that these functions obey the BCs

$0(xfn(xs); Xs)=0, $0(xJax(xs); xs )=0 (6)

at the boundary points {^min(^s), xmax(xs)}=dQXf (xs), of the interval QXf (xs). The eigenfunctions satisfy the orthonormality condition

-f v-s/

<$i = J $i(Xf; Xs)$ j(Xf; Xs)dXf(7)

X^in (Xs )

Here e(xs) : e1(xs) < ... < £jmax (xs) < ... is the desired set of real eigenvalues.

During the simulation the etalon potential V0(xf; xs,g(xs)) in Eq. (5) will be chosen as V0(xf; xs,g(xs)) = V(xf ,xs) or calculated separately for different values of xs G QXs from the conditions

Xm (x s)

min / (V(xf ,Xs)—Vo(xf;Xs,g(xs))2 dxf. (8)

g(xs) J

xmin(xs)

If this parametric eigenvalue problem has no analytical solution, then it is solved numerically by the FEM: at n — 1 = 1 using the program ODPEVP [6] and at n — 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 x(t)(xs,E) =

(—1 zsibdï ^ ^ 1+

+ fs2(Xs)dQ(Xs) + 1 dfs2 x(i),g )=0. (9)

fsl(Xs) d^s fsl(Xs) dxs J s '

Here I, W(xs) and Q(xs) are the jmax x jmax matrices

W^(xs Hij(xs)+Vij(xs), Iij = % (10)

fs3(Xs) J fsl(Xs) J J

max / \ X f (X s)

0Ф

0xs dxs

TT / \ TT i \ I d$i(xf ;xs) дФу (xf ;xs) , , Нгз (xs) = Hji(xs)= -jrf--^-f-dxf , (11)

xffin (x3)

xf (xs)

Ягз(xs)=-Qji(xs)=- j Ф^;xs)дФdx]x's) dxf. (12)

xffin(x3)

The effective potentials Vij(xs)=Vji(xs) are calculated by integrating the difference 5V(xf ,xs) = V(xf, xs)-V0(xf ;xs, g(xs)) with the basis functions:

max Xf (xs)

vij (xs)

Vij (xs)= J Ф^/ ;xs)5V (xf ,Xs)Фj (xf ;xs)dxf. (13)

afin(as)

If the etalon potential is equal to the original one, V0(xf ;xs, g(xs))=V(xf ,xs) then Vij(xs)=0. If the original potential V(xf ,xs) is presented in the tabular form, e.g., as in Ref. [4], then the etalon potential V0(xf ;xs, g(xs)) can be considered as a certain approximation or interpolation of the original one, and the approximation error 5V(xf ,xs) = 0. In this case the etalon potential Va(xf ;xs, g(xs)) can be determined minimizing 5V(xf ,xs) in accordance with the condition (8) of the least squares method. Then the integrals Vij (xs) of the approximation error 5V(xf ,xs) 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 E : E\ < E2 < ... < Ev < ... obey the BCs at the points xts={x'min,x'max}=dQXB, bounding the interval qXb and satisfy the orthonormality conditions

X^(xs)=0,xs=xmin,xrx, / fsi(xs)(x{l)(xs))TX(j)(xs)dxs=öij. (14)

x

3. FEM algorithm for solving the parametric 2D BVP

Let us consider a boundary value problem for the parametric self-adjoint 2D PDE in the domain üx,y = (xmin,xmax) x (ymin, ymax)

(D(x, y; z)-Ei(z)) фi(x, y; z~)=0, 1 д д 1 д д (15)

D=D(x,У;z)=-—öxh(x, у)ox-мщу) д~уh(x,у)д~у+u(x,y; z),

with the Dirichlet and/or Neumann boundary conditions

lim f 2(x, у)0V;Z) = 0 or Ф^г, y;z) = 0, y G [ymin, ?/max], x^xt dx

,дфi(x, y;z)

lim f5(x, y)--- = 0 or Ф^, yt ;z) = 0, x G [xmin,xmax],

y^yt oy

where t = min, max. Here z e Qz = [zmin, zmax] is a parameter, the functions f1(x) > 0, h(x,y) > 0, f3(x) > ° f4(y) > 0, fh(x,y) > ° and дxf2(x,y), dy f5(x,y), U(x,y; z), 9U(a'z 'z) and are continuous and bounded for (x,y) e Q,xyy. Also assume that

the BVP (15), (16) has only the discrete spectrum, so that e(xs) : e1(xs) < ... < ejmax(xs) < ... is the desired set of real eigenvalues. The eigenfunctions satisfy the orthonormality conditions

x------y

il$j) = J J h(x)h(y)$i(x,y; z)$j(x,y; z)dx dy=5ij. (17)

min ...mm

Xmin y

The FEM calculation scheme is derived from the variational functional

Xmax ymax

, /mt1 i , i \ d § i(x, y; z) d $ .j(x, y; z)

ID—e(z)l$j)=Jdxj dy(U(y)f2(x,y) ^ j qx ' +

xmin ymin

+ h(x) f, y) d$i(x,y; z) d$j(x,y; z) + f3(x) 5 ' dy dy

+ $i(x, y; z)fi(x)f4(y)(U(x, y; z) — e(z))$j(x, y; z)j . (18)

1. The domain A = [^min,^max] x [ymin,ymax] is covered by the system of n x m subdomains A^- = [xi-1,xi] x [yj-1,yj] in such a way that A = (J™=1 UAij. In each subdomain A^ the nodes {xP,r}P=0 and {yjj,r}P

j,r> r=0

hX

hy

xPr=Xi-1 + — r, hX=Xi—Xi-1, yP r=yj-1+--r, hy =yj—yj-1,

and the Lagrange elements {$p,r(x)}p=0 and {^p,r(y)}p=0

r=0 ,P

<r (x)= n TP _ % , (y)= n

P P _ ^

s=0,s=r 'l,r 2,,s

P

y — yP,s

P P

s=0,s=r yj,r — yj,s

are determined. By means of the Lagrange elements (x) and (y), we define the set of piecewise polynomial functions Np(x) and Mp(y) as follows:

iP

4>Pt0(x), x e A1j,

0 ' x e A1j,

$tr(x), X e Aij ,

0, x e Aij,

1=0,

l=r+p(i—1), r=1,p—1,

NP (x) =

€,p(x), x e Aij,

C1

n^^ x e Ai+1j ,

0, xe AijU Ai+1j,

$n,P(x), X e Anj ,

0,

xi A

nj

l=ip, i=1, n—1, l=np,

for any and

M[(y) =

tf,o(v), 0,

(y),

0,

€P(y),

^+1,0 (y), 0,

0,

y G Ai1,

y G An, y G Aa, y G An, y G Aij,

y G Kj+U

y G Ai^ Aij+1,

y G Aim,

y G Aгm,

1=0,

l=r+p(j-1), r=1,p-1,

1=jl^ j=1,m-1,

l=mp,

for any .

The functions {Nf (x)}™f0 and {Mlf(y)}]=p0 form a basis in the space of polynomials of the p-th order. Now, the function $(x, y;z) G T^hh ~ 'H1(fihx,hy) is approximated by a finite sum of piecewise polynomial functions Nf(x) and Mf (y)

mp np

\x, y;z) = ^ ^ £Jv (*)K (X)K (y)-

(19)

lu=0 lx = 0

2. The domain fi(x, y) = U^=1 Aq, specified as a polygon in the plane (x = z1,y = z2) G 'R2, is covered with finite elements, the triangles Aq with the vertices (Z11, Z21), (Z12, Z22), (-^13, ^23) (here Zik = Zik;q, i = 1, 2, k = 1, 3, q = 1,Q). On each of the triangles Aq (the boundary is considered to belong to the triangle) the shape functions (pf(z1, z2) are introduced. For this purpose we divide the sides of the triangle into p 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 p, 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 p.

Three straight lines from different families intersect in one point Ai G Aq, which will be numbered by the triplet (n1,n2,n3), ni > 0, n1 + n2 + n3 = p, where n1, n2 and n3 are the numbers of the straight lines passing parallel to the side of the triangle that does not contain the vertex (z11, z21), (z12, z22) and (z13, z23), respectively. The coordinates of this point zl = (z1l, z2l) are determined by the expression (z1l, z2l) = (Z11, Z2\)n1/p + (Z12, Z22 )n2 /p + (Z13, Z23)n3/p.

As shape functions we use the Lagrange triangular polynomials <pf (z1, z2) of the order p that satisfy the condition (pf(z1 y, z2l') = 5ly, i.e., equal 1 in one of the points Al and zero in the other points.

In this method the piecewise polynomial functions Nf(z1, z2) in the domain fi are constructed by joining the shape functions (pf(z1, z2) in the triangle Aq:

Nf (Z1, Z2) = ^( Z1, Z2 ),Al G Aq ;0,AlG AqJ

and possess the following properties: functions Nf(z1, z2) are continuous in the domain fi; the functions Nf(z1, z2) equal 1 in one of the points Al and zero in the rest points; Nf(z1y, z2y) = 5ll' in the entire domain fi. Here I takes the values I = 1,N.

The functions Np(x, y) form a basis in the space of polynomials of the p-th order. Now, the function $(x,y; z) G Jj ~ H1(Qhx,hv) is approximated by a finite sum of piecewise basis functions Nf (x,y)

N

<-l

$h(X,y; z) = Y, il(z)^P (x,y). (20)

1=1

The vector function ^h = {¡¿?lv(z)}™= o"p=o or £h = (z)}{=1 has a generalized first-order partial derivative and belongs to the Sobolev space H1(Qhx,hv) [10]. After substituting the expansion (19) or (20) into the variational functional (18) and minimizing it [10,11], we obtain the generalized eigenvalue problem

AP£h = eh BP£h. (21)

Here AP is the stiffness matrix; BP is the positive definite mass matrix; £h is the vector approximating the solution on the finite-element grid; and eh is the corresponding eigenvalue. The matrices AP and BP have the following form:

■m n m n

AP = EEapi, Bp = ££bp, or Ap={ap,}N=1, Bp={6p,}N=1,

=1 =1 =1 =1

where the local matrices ap^ and bp^ are calculated for rectangular elements as

+1 +1

+ f3(x) h(X,y)^,q (X)^,r (X)(hy )2 dVy dVy +

1 h'X hy-

+ h(x)h(y)U(x, y; z)4>pq(x)^p,r(xWp^(y)^p„(y)f^Ty d^x dVy,

+1 +1 IV

/( h'x h ■

J h(x)f4(y)$q(x)4Ptr(xWhlAyWp,v(f)yydVxdVy, -1 -1

x = Xi-1 +0.5hX(1 + rjx), y = yj-1 +0.5h| (1 + r]y), y,,v,q,r = 0,p, or the matrix elements ap,, and bp,, are calculated for triangular elements as

(22)

apu, = / ¡1(x) ¡4(y)^p(z1, Z2)pp (Z1, Z2)U (Z1, Z2 dz2+

A

, iff \f f \d^i(z u Z2)d^p, (ZU Z2)

+ f 2(x, y)U(y)—^--^-dZ1 dZ2+

a,

, f h(x) t , Z2)d^p(Z1, Z2) , ,

+ y)—^--—-- dz 1 d Z2 =

J f3(x) d Z2 d Z2

A„

= IJI J h(x) h(y)vP(z[, (4, 4)U (Z1, z2;z)dz[ d4+ A

+IJI-1 /№.v)M {^n-^f-^^+(^-«1)^^)

x ( fe^) ^f^ + dz1 d^ +

+ I J- / ^^ " T + ^)

A

X (( ,11-,13)^^+(, 12-,11)^f^) d,1 d4, (23)

tfv = J h (x) f 4 (y)VPl( Z1, (Z1, Z2)dZ 1 d ,2 =

A,

= Ulj h (x) f4(yM(zf1, z2)vf (z[, z2)dz1 d zf2, (24)

A

where Aq is a triangular element with the vertices (z11, z21)(z12, z22)(z13, z23), A is the triangle with vertices (0,0)(1, 0)(0,1), and

J =

- 1 - 1

dz[ dz'2

dz2 dz2

- 1 - 2

-Z22Z13 + Z22Z11 + Z21Z13 + Z23Z12 - Z23ZU - Z21Z12

is the Jacobian of the transformation from the global frame (x = z1,y = z2) to the local one (z\, z'2); d^dz2 = Jdz\dz'2 and (z1, z2) are expressed via (z1, z'2) by the relations:

,1 = ,U + ( ,12-,11), 1 + ( Z13 - ,n) Z'2, Z2 = ,21 + ( ,22-,21),1 + ( ,23 - ,21) (25)

In this case we have explicit expression for shape functions (pf(z\, ,2):

«1-1 / / / I «2-1 / /I V.3-1 / I,

VP(^ Z2) = TT 1 - - Z2 -n1/p TT Z1 n2 /p TT ,2 - n3/p

l ^ 2 , n n1 /p -n'1/p n2 /p -n'2/p *-\n3/p -n'3/p'

«1=0 n2 = 0 n3=0

Remark. We start from the initial (global) coordinate frame (x = z1,y = 2 ) in which the coordinates of vertices of the triangle Aq are equal to (Z", z21), (z12, z22), (z13, z23), and the local coordinate frame (x = z[,y = ,2) in which the coordinates of the same vertices of the triangle A are equal to ( 11, 21) = (0, 0), ( , 12, ,22) = (1, 0), ( , 13, ,23) = (0, 1).

We seek the relation between the global and the local coordinates of the triangle vertices in the form

,1 = C10 + 011^1 + C12Z2, Z2 = C20 + 021^1 + C22Z2. (26)

X

X

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 c**:

Z11 = C10, Z21 = C20, Z12 = C10 + C11, Z22 = C20 + C21, ZVi = CW + C12, Z23 = C20 + C22.

Substituting the calculated coefficients into Eq.(26), we arrive at the formula (25) for the transformation of coordinates (z1, z2) ^ ( z[, z'2), from which we express the inverse transformation of coordinates ( z[, z'2) ^ (z1, z2)

Z[ = — J-1( Z23Z11— Z21Z13) + J -1( Z23 Z21) Z1— J-1( Z13 — Z11) Z2,

Z*2 = J-1( Z22Z11— Z21Z12) + J-1( — ¿22 + Z21) Z1+J-1( Z12- ^11) Z2,

J =(-Z22Z13 + ^22^11+ Z21Z13+Z23Z12 -Z23Z11-Z21Z12),

where J-1 is the Jacobian of the inverse transformation from the local frame (z'1, z'2) to the global one (x = z1,y = z2): d z[dz2 = J-1dz1 dz2.

The integrals (22) and (23)-(24) are evaluated using the Gaussian quadrature of the order p +1. Note that in the present approach the maximum half-band width of the matrices Ap and Bp is small compared to their dimension and is not greater than (p+1)(p+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 banded-matrix 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 Ap in Eq. (21) is not positively defined, the problem (21) is replaced with the following problem:

Ap = £h Bp Ap = Ap - aBp. (27)

The number a (the shift of the energy spectrum) is chosen such that the matrix Ap is positive defined. The eigenvector of the problem (27) is the same, and eh = eh + a.

4. 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 z, we find that dz$i(x, y; z) is a solution of the following boundary problem

(D(x, y;z)—ez))

dÇ^x^ y;z) dz

d

— (U(x, y;z)—£i(z))

$iy;z),

(28)

д2$í(x, y;z) d®^^ y;z)

lim M^ y)-^- = 0 or -^- = ^ y e [Уmin, y max],

x yx-t OXOZ OZ

,д2$i(x, y;z) d$i(x, yt;z)

lim J 5 (x, y)-^T- = 0 or --- = 0, X e [Xmin,Xmax],

y^yt oyoz OZ

where t = min, max. The parametric BVP (28), (29) has a unique solution, if and only if it satisfies the conditions

#max Umax

f j dx dyf1(x)f4(y)(Mx, y;z)) °U ^(x, y;z), (29)

^min Umin ^max Umax

J J dx dyfr(x) f4(y)(*i(x, y; z)) d^(xdzy ]z) =0. (30)

^min Umin

Below we present an efficient numerical method that allows the calculation of dz^¿(x, y; z) 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

Hij(z)=Hji(z)= dx f1 (x) dy f4(y)

д$i(x, y;z) ^^x^ y;z) I 'xyjiyy' 3~z

^min

W ^ ~> I j ; i » -z --Z

(31)

#max ^max

Ql3(z)=-Q3l(*)=-/ dx h(x) j dy f4(y)^(x, y,z) ^^V]Z).

^min ^min

The boundary problem (28)-(30) is reduced to the linear system of inhomogeneous algebraic equations with respect to the unknown -$,h/-z:

4 - (AP = ^ - - (-A - £bp) * <32)

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

«f = '• (f) B'i" = 0. th = ^<h <33)

Then the potential matrix elements H"(z) and Qhj(z) (31) corresponding to Eqs. (11)-(12) can be calculated using the formulas

T

HhW= (f) BPf• Qh« = - (<h)TBPf

iS- • QW = - Un B"^-. (34)

Since eh 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

Lv = b, Lw = d, (35)

with the non-degenerate matrix L and the right-hand sides b and d

Lss,, (s -S)(s'-S) = 0, 8Sa,, ( s -S)(s' - S) = 0,

{

h i ba, S = S, f LsS, S =

hs [0, s = S, ds [0, s = S,

where S is the number of the element of the vector Bp^h having the greatest absolute value.

Step k2. Evaluate the coefficient 7

7 = , 71 = vTBpÇh, 72 = wT BpÇh, Ds = (Bp^h)s. (36)

(Ds " 72)

Step k3. Evaluate the vector dz$,h

d^h ( Vs " 7Ws, s = S, dz = 1 7, s = S.

(37)

From the above consideration it is evident, that the computed derivative has the same accuracy as the calculated eigenfunction.

Let D(x, y; z) in Eq. (15) be a continuous and bounded positively defined operator on the space V1 with the energy norm, £i(z), ^¿(x, y;z) £ V2 being the exact solutions of Eqs. (15)-(17), and z), y; z) £ V1 being the corresponding numerical solutions.

Then the following estimates are valid [10,12]

^¿(z) - e'l(z)\^ cih2p, ||^(x, y;z) - $£(x, y;z)\\Q < c2hp+1, (38)

where

^max ?/max

■ •¿(x^lS = / /A(x, »)dxA(x,,)d»*(x, ^¿(x,*,),

xmin ymin

h is the largest distance between any two points in Aq (see [12], p. 161), p is the order of the finite elements, i is the number of the corresponding solutions, and the constants c1 and c2 are independent of the step h.

The following theorem can be formulated.

Theorem 1. Let D(x, y;z) in Eq. (15) be a continuous and bounded positively defined operator on the space V1 with the energy norm. Also let dzU(x, y; z) be continuous and bounded for each value of the parameter z. Then for the exact values of the solutions dze^(z), dz^¿(x, y;z) £ V2, H¿j(z), Q¿j(z) from (28)-(31) and the corresponding numerical values dzell(z) estimates are valid:

dz$h(x, y;z) eH1

Hhj(z), Q%(z) from (32)-(34) the following

de i(z) de h(z)

d

d

< c3h2p,

d$i(x, y;z) d$h(x, y;z)

d

d

< c4hp+1,

\ Qij(z) " Qhi(z)\ < c5h2p, \Hij(z) " Hih,(z)\ < c6h2p,

(39)

where h is the largest distance between any two points of the finite element Aq, p is the order of finite elements, i, j are the numbers of the corresponding solutions, and the constants c3, c4, c5 and c6 are independent of the step h.

0

The proof is straightforward following the proof schemes in accordance with [10,13].

5.1.

5. Benchmark calculations 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 (x, y) e 0(x, y)

d2 d2

-m - d? - V ^ = 0

(40)

with the Dirichlet or Neumann conditions at the boundary d0(x, y) of the region 0(x, y). In both cases the eigenvalues si are integer [14-16].

The BVPs were solved on the uniform finite-element grid composed of n2 equilateral triangles with the side equal to 4^/(3n). In Table 1 we present the estimations of the finite-element scheme order p depending on the size of elements.

Table 1

Discrepancies Ssi(n) = (n) — £i | of the first e\ =3 and the fourth e4 = 3 eigenvalues of the BVP(40) with Dirichlet and Neumann boundary conditions, respectively, and the Runge coefficients Ru« = log2((Ssi(n) — ¿£i(2n))/(Ô£i(2n) — Ô£i(4n))) for the schemes of the first p =1 and the second p = 2 orders; n is the number of elements

p n Ô£i (n) Ô£!(2n) Ô£i (4n) Rui

1 6 0.06914677126132 0.01717343333794 0.00428595766335 2.011

2 3 0.00470310603030 0.00030790240009 0.00001932528605 3.928

p n Ô £ 4 (n) Ô£4 (2n) Ô£4 (4n) Ru4

1 6 0.06914677126135 0.01717343333808 0.00428595766371 2.011

2 3 0.00470310603029 0.00030790240005 0.00001932528593 3.928

They are seen to correspond rigorously to the theoretical estimations (38) of the order 0(h2p) 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 1. The first eight eigenfunctions of the BVPs (40) in the triangle region with the Neumann boundary conditions. The eigenvalues s% = 0,1,1, 3, 4,4, 7, 7 at

i = 1,..., 8.

Figure 2. The first eight eigenfunctions of the BVP (40) in the triangle region with the Dirichlet boundary conditions. The eigenvalues £i = 3, 7, 7,12,13,13,19,19

at i = 1,..., 8.

5.2. Solution of the parametric 2D BVP for the oscillator potential

We consider the parametric 2D BVP (15)-(17) at fi (x) = f2 (x, y) = f3 (x) = f4 (y) = f5 (x,y) = 1, with the potential function

U (x, y ; * ) = U (Z! ,Z2 ; ¿3 ) = (Z! - Z3 )2 + (3/2)2 ^

(41)

which has the known spectrum Sj={j1 ¿2} = (2(ji — 1) + 1) + 3/2(2(j2 — 1) + 1) and the eigenfunctions (zi ,z2; z3) = ^ -i (zi; z3)$j2-i (z2;0), where ^ (zi; z3) and

-i (z2; 0) are determined in (45) at ui (xs) = 1, z0i (xs = z3) = z3 and (xs) = 3/2, ¿02(^5) = 0, respectively. The matrix elements Qi={il,¿2},j={ji,32} and #¿={¿1,¿2},j={ji,32} are calculated in the analytical form

Qij (xs )=^2 32 sign(ii -H )

V2n

-¿11,1

Hij (xs )=ô.

(2n+1)

12 32

2

ôji -

h ,0

\sJn(n 1) 2

-ii |,2^

(42)

where n= max(ii ,ji) —1.

Solving the BVP, we use 294 finite elements, i.e. equilateral triangles with the side equal to 1 that form a regular hexagon with the vertices {(7 cos nn/3, 7 sin ^n/3)}^=0. In each element the interpolation polynomials of the six order were applied. Stiffness and mass matrices with dimensions of 5167 x 5167 were used. The calculated eigenvalues e1} are: ^ = (2.500000001, 4.500000002, 5.500000005, 6.500000014, 7.500000025), and the exact ones ^ are: ^ = (2.5, 4.5, 5.5, 6.5, 7.5).

The calculated matrices Q^ and H^, i,j = 1,..., 5 from Eq. (34) of the parametric 2D BVP with the potential function (41) are

Qij =

/-0.00000000 0.70710678 -0.00000000 -0.00000000 0.00000000\

-0.70710678 -0.00000000 -0.00000000 1.00000000 0.00000000

-0.00000000 0.00000000 0.00000000 -0.00000000 0.70710678

0.00000000 -1.00000000 0.00000000 -0.00000000 0.00000000

-0.00000000 0.00000000 -0.70710678 0.00000000 0.00000000

Tjh

0.50000000 -0.00000000

-0.00000000 1.49999999

0.00000000 -0.00000000

-0.70710678 -0.00000000

0.00000000 -0.00000000

0.00000000 -0.70710678 0.00000000

-0.00000000 -0.00000000 -0.00000000

0.49999999 -0.00000000 0.00000000

-0.00000000 2.50000001 -0.00000000

0.00000000 -0.00000000 1.49999999

From the comparison of the calculated and the exact values of si and Q^, H^ from (42), one can see that the achieved accuracy of results is of the order of 8-9 significant digits.

The parametric surface eigenfunctions and their derivatives with respect to the parameter under the dipole splitting are shown in Fig. 3.

2* * 1 • a « 1 * *

-4 -22-■■■ - n 9 v 2 4 * • « • * *

■ ■ ■ -4 -2- 4

2-■ "

■ ■ ■ ■

-22-

1

\v2 4 -4

2

■ M

: «1 4-4

2

4

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.

1

4

5.3. The solution of 2D BVP for the CSv oscillator potential

Let us consider the 2D BVP (15)-(17) for x=a22, y=a20 with the potential function U(x, y) = 2mV(x, y) and the spectral parameter £ = 2mE, and fi (x) = f2 (x, y) = fs (x) = f4 (y) = f5 (x, y) = 1. The quadrupole potential energy is approximated by the quartic potential:

V (a22 ,a20 )=ci (a222+alo)+C2 (a^o-a^/3)+cs (a^ +y 2 )2+Co (43)

We use the set of parameters ^=-120, c2 =240, cs=1200, c0=65/16 that provide a crude approximation for the shape of 156 Gdg2, which has been fitted in the following points1: the minima at (a22,a20)=(0,1/4), V(0,1/4)=0; the maxima at (a22,a20)=(0, 0), V(0, 0)=65/16; and the saddle point at (a22, a2o)=(0,-1/5), V(0,-1/5)=729/400 (see Fig. 4).

We choose the mass parameter to be m=B2=124. Thus, there are ground and doubly degenerate excited states, localized in three wells.

The BVP was solved using three algorithms (and the forth in Section 5.4): 1. The solution of the BVP was calculated using the FEM scheme from Section 3 on the rectangular grid [-0.4, -0.3,..., 0.4] x [-0.4, -0.3,..., 0.4] with Lagrange interpolation polynomials of the order p = 12. The first 18 eigenvalues were calculated with 9 significant digits and are presented in Table 2.

1 The fitting points in our parametrization are related to those of Ref. [4] as a22 = V2a22, &20 = «20 •

Figure 4. The potential energy of the C^v oscillator having the quadrupole shape, the eigenfunction (a22,a2o) of the ground state (irr. A1) and the degenerate eigenfunctions (a22 ,a20) (irr. E1) and (a22 ,a20) (irr. E2)

Table 2

The first energy levels of the 2D BVP using triangular e^ and rectangular finite elements and the Kantorovich method using jmax=28 parametric basis functions, classified by irrs (Class.) of the C3v point group, Ei (MeV)

i ci £h Eh=eh/2m Class.

1 381.754344 381.754355 381.754351 1.53933206 A1

2 387.240633 387.240644 387.240641 1.56145419 E1

3 387.240633 387.240646 387.240641 1.56145419 E2

4 617.024951 617.024967 617.024963 2.48800388 E1

5 617.024951 617.024989 617.024963 2.48800388 E2

6 667.104970 667.105020 667.104992 2.68993948 A2

7 695.166557 695.166590 695.166575 2.80309103 A1

8 785.680037 785.680100 785.680078 3.16806483 E1

9 785.680037 785.680136 785.680078 3.16806483 E2

10 898.045395 898.045497 898.045434 3.62115094 A1

11 915.823095 915.823200 915.823167 3.69283535 E1

12 915.823095 915.823309 915.823167 3.69283535 E2

13 993.158636 993.158784 993.158708 4.00467221 E2

14 993.158636 993.158872 993.158708 4.00467221 E1

15 1063.73690 1063.73709 1063.73692 4.28926178 A1

16 1119.21670 1119.21668 1119.21649 4.51296973 A2

17 1174.72840 1174.71177 1174.71166 4.73674057 E1

18 1174.75531 1174.71183 1174.71166 4.73674057 E2

2. In the solution of the BVP we used 54 finite elements in the form of equilateral triangles with the side equal to 1/6, forming a regular hexagon with the vertices {(0.5cosnn/3, 0.5sin™/3)}^=0. In each finite element the interpolation polynomials of the fifth order were applied. The stiffness and the mass 721 x 721 matrices were used. The calculated eigenvalues are presented in Table 2.

3. The problem (1)-(43) was also solved using the Kantorovich method implemented in the FEM program KANTBP 2 [5], with jmax=28 parametric basis functions calculated using the FEM program ODPEVP [6]. The solution was calculated in the domain -\/a20+a22< 1/2 with Dirichlet BCs at the boundary y/'a|0+®12=1 /2 using the scheme presented in Section 2. The calculations were performed in the case of (xf, xs)= (a22, a20) as well as (xf ,xs)=(a20, a22) on the FE grid {-1/2, -1/3, -1/6, 0,1/6,1/3,1/2} with Lagrange interpolation polynomials of the order p = 12.

The point symmetry group C3v 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.

5.4. Parametric surface functions for KM in the analytical form

Let us consider the BVP for Eq. (5) with the etalon potential V0(xf ;xs, g(xs)):

( C)2 \

I -Cx1+Vo(xf;xs,g)-£i(xs) I $i(xf;xs)=0,

Vo(xf ,xs, g)=Vo(xs)+u2(xs)(xf - Zo(xs))2,

(44)

where g(xs) is the set of parameters, g(xs) = {V0(xs), u2(xs), z0(xs)}. In the considered case the parametric eigenvalue problem (5)-(7) has an exact solution, i.e., the parametric eigenfunctions (xf; xs) and potential curves £i (xs) are expressed in the analytical form

£i (xs) =Vo(xs)+u(xs)(2(i-1)+1), (xf; xs) = W ^¡ls) exp(-w(xs)(xf-zo(xs))2/2)), (45)

k

(Xf Zo((^)1 _ (Xf ;Xs) _2 (Xf ]Xs).

vi—l Vl—1

The integration in the effective potentials (11) with the basis functions (45) is carried out analytically, which yields the expressions

^ , x . v2n^(xs) dz0(xs) , v/n(n_I) dw(Xs) J. A Q,J(xs) = Slgn(J_) I-2--d^5I_I'1--4--d^Ôl_,2j ,

n2+n+1 -\/n(n_1)(n_2)(n_3) \ /dw(xs)x 2

h*j(xs)= l s^X) Si-i '0--16^X7)-^H-dX-J +

u(Xs)(2n+1) u(Xs)\/n(n_1) \ f dz0(Xs)\2

+ I-2-s*-i'0--2-_

I n\/2n i/2n(n_1)(n_2l \ dzo(xs) dw(xs)

^ï^xsTu+ "13 ) dir:

where n=max(i,j)—1. The effective potentials are calculated by integrating the difference V(xf, xs)-V0(xf ; xs) with the basis functions:

X f (xb)

Vij (xs)= J $i(xf ; xs)(V (Xf ,xs)-Vo(xf ; xs))$j (xf ; xs)dxf.

xmin(xs)

During the simulation the adiabatic parameters V0(xs), u(xs), z0(xs) of the etalon potential (44) were found from the conditions

xf (xB)

J (V (Xf Xs)-Vo(xf ; ^s))2 dxf. (46)

. \ \ , - / ~ s

Vo(x3),u2(x3),zo(x3)

min

), zo(xs

xfin(xB)

For the potential (43) from the condition (46) we have /960

u(xs)=V'2'm^J~'9^+240xs+2400x2s, z0(xs)=0, at xs=a2o, Xf=a22,

/ ^ ^ /960 9 / ^ 7(20®2-1) u(Xs)=v2m^ —+2400x9s,z0(xs)=-§0(35^9+2)' at xs=^99, Xf=090.

We performed the calculations for these parameters in the case when (xf, xs)=(a99, a90), as well as (xf ,xs)=(a90,a92). 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.

6. 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.

We proposed the construction of parametric surface functions in the analytical form as eigenfunctions of the etalon equation (44), which provides the solution of the 2D BPV with given accuracy and reduces the expenditures of computer resources compared to the conventional basis, numerically calculated using FEM. One can construct parametric functions using different types of etalon potentials, e.g., that of the two-center problem with harmonic oscillator potentials [19]. This approach can be generalized for the BVP in a multidimensional domain using, e.g., the multistep Kantorovich method [3].

The proposed algorithms and codes can be adapted and applied to the analysis of quantum transparency effect, to the study of resonance three-body scattering problems, the quantum diffusion of molecules, the penetration of micro-clusters through surfaces, the fragmentation mechanism in producing very neutron-rich light nuclei, and heavy-ion collisions, as well as the microscopic study of tetrahedral-symmetric nuclei.

References

1. S. I. Vinitskii, L. I. Ponomarev, Adiabatic representation in three-body problem with coulomb interaction, Soviet Journal of Particles and Nuclei 13 (1982) 557-587, in Russian.

2. G. A. Parker, R. T. Pack, Quantum Reactive Scattering in Three Dimensions Using Hyperspherical (APH) Coordinates., VI. Analytic Basis Method for Surface Functions, Journal of Chemical Physics 98 (1993) 6883-6896.

3. A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, A. G. Abrashkevich, V. L. Der-bov, Numerical Solution of Elliptic Boundary-Value Problems for Schrödinger-Type Equations Using the Kantorovich Method, Mathematical Modelling and Geometry 2 (2014) 54-80.

4. A. Dobrowolski, K. Mazurek, A. Gozdz, Consistent Quadrupole-Octupole Collective Model, Physical Review C 94 (2017) 054322.

5. O. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, A. G. Abrashkevich, KANTBP 2.0: New Version of a Program for Computing Energy Levels, Reaction Matrix and Radial Wave Functions in the Coupled-Channel Hyperspherical Adiabatic Approach, Computer Physics Communications 179 (2008) 685-693.

6. O. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, A. G. Abrashkevich, ODPEVP: A Program for Computing Eigenvalues and Eigenfunctions and Their First Derivatives with Respect to the Parameter of the Parametric Self-Adjoined Sturm-Liouville Problem, Computer Physics Communications 180 (2009) 1358-1375.

7. A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, A. G. Abrashkevich, POTHEA: A Program for Computing Eigenvalues and Eigenfunctions and Their First Derivatives with Respect to the Parameter of the Parametric Self-Adjoined 2D Elliptic Partial Differential Equation, Computer Physics Communications 185 (2014) 2636-2654.

8. K. Kumar, M. Baranger, Complete Numerical Solution of Bohr's Collective Hamil-tonian, Nuclear Physics A 392 (1967) 608-652.

9. A. Dobrowolski, A. GoZdZ, K. Mazurek, J. Dudek, Tetrahedral symmetry in nuclei: New predictions based on the collective model, International Journal of Modern Physics E 20 (2011) 500-506.

10. G. Strang, G. J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs, New York, 1973.

11. K. J. Bathe, Finite Element Procedures in Engineering Analysis, Englewood Cliffs, Prentice Hall, New York, 1982.

12. E. B. Becker, G. F. Carey, J. Tinsley Oden, Finite Elements. An Introduction. Vol. I, Englewood Cliffs, Prentice Hall, New Jersey, 1981.

13. O. Chuluunbaatar, The scientific doctoral thesis (2010).

14. F. Pockels, Uber die partielle differential-gleichung Au+k2u = 0 und deren auftreten in der mathematischen physik, B. G. Teubner, Leipzig, 1891.

15. M. V. Berry, M. Wilkinson, Diabolical Points in the Spectra of Triangles, Proceedings of the Royal Society A: Mathematical, Physical & Engineering Sciences 392 (1984) 15-43.

16. A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, V. L. Derbov, Solution of the Boundary-Value Problem for a Systems of ODEs of Large Dimension: Benchmark Calculations in the Framework of Kantorovich Method, Bulletin of PFUR. Series: Mathematics. Information Sciences. Physics 3 (2016) 31-37.

17. J. F. Cornwell, Group Theory in Physics, Academic Press, New York, 1984.

18. I. N. Belyaeva, N. A. Chekanov, A. A. Gusev, V. A. Rostovtsev, Yu. A. Ukolov, Y. Uwano, S. I. Vinitsky, A MAPLE Symbolic-Numeric Program for Solving the 2D-Eigenvalue Problem by a Self-Consistent Basis Method, Lecture Notes in Computer Science 3718 (2005) 32-39.

19. J. M. Eisenberg, W. Greiner, Nuclear Theory. Vol. 1, North-Holland, Amsterdam, 1970.

УДК 517.958, 530.145.6, 519.632.4 Б01: 10.22363/2312-9735-2017-25-1-36-55

Алгоритмы для решения параметрической самосопряжённой эллиптической краевой задачи в двумерной области методом конечных элементов высокого порядка точности

А. А. Гусев*, О. Чулуунбаатар*1", С. И. Виницкий**, В. Л. Дербов§,

А. Гуждж^

* Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, г. Дубна, Московская область, Россия, 141980 ^ Институт математики Монгольский национальный университет, Улан-Батор, Монголия * Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, Россия, 117198 § Саратовский государственный университет, г. Саратов ^ Институт физики, Университет им. М. Кюри-Склодовска, Люблин, Польша

Рассмотрены вычислительные схемы решения краевых эллиптических задач в рамках метода Канторовича — редукции эллиптической краевой задачи к системе обыкновенных дифференциальных уравнений второго порядка с использованием поверхностных функций, зависящих от независимой переменной системы обыкновенных дифференциальных уравнений как параметра, вычисленные как методом конечных элементов, так и пробных параметрических поверхностных базисных функций, вычисленных в аналитической виде.

Предложены новые вычислительные схемы и алгоритмы для решения параметрической самосопряжённой эллиптической краевой задачи в двумерной области, используя метод конечных элементов высокого порядка точности с прямоугольными и треугольными элементами.

Комплексы программ, реализующие алгоритмы, вычисляют с заданной точностью собственные значения, собственные функции и их первые производные по параметру, связанные с параметрической самосопряжённой краевой задачей для эллиптических дифференциальных уравнений с условиями Дирихле или Неймана на границе в конечной двумерной области, а также потенциальные матричные элементы — интегралы от произведения собственных функций и их первых производных по параметру. Параметрические собственные значения (так называемые потенциальные кривые) и матричные элементы, вычисленные с помощью комплекса программ, можно применять для решения задачи на связанные состояния и многоканальной задачи рассеяния для системы обыкновенных дифференциальных уравнений второго порядка с помощью метода Канторовича.

Эффективность предложенных схем расчёта и алгоритмов демонстрируется решением двумерных эллиптических краевых задач, описывающих квадрупольные колебания в коллективной модели атомного ядра.

Ключевые слова: параметрические эллиптические краевые задачи, метод конечных элементов, метод Канторовича, системы ОДУ второго порядка

Литература

1. Виницкий С. И., Пономарев Л. И. Адиабатическое представление в задаче трех тел с кулоновским взаимодействием // Физика элементарных частиц и атомного ядра. — 1982. — Т. 13. — С. 1336-1418.

2. Parker G. A., Pack R. T. Quantum Reactive Scattering in Three Dimensions Using Hyperspherical (APH) Coordinates., VI. Analytic Basis Method for Surface Functions // Journal of Chemical Physics. — 1993. — Vol. 98. — Pp. 6883-6896.

3. Numerical Solution of Elliptic Boundary-Value Problems for Schrodinger-Type Equations Using the Kantorovich Method / A. A. Gusev, O. Chuluunbaatar, S. I. Vinit-sky, A. G. Abrashkevich, V. L. Derbov // Mathematical Modelling and Geometry. — 2014. — Vol. 2. — Pp. 54-80.

4. Dobrowolski A., Mazurek K., Gozdz A. Consistent Quadrupole-Octupole Collective Model // Physical Review C. — 2017. — Vol. 94. — P. 054322.

5. KANTBP 2.0: New Version of a Program for Computing Energy Levels, Reaction Matrix and Radial Wave Functions in the Coupled-Channel Hyperspherical Adiabatic Approach / O. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, A. G. Abrashkevich // Computer Physics Communications. — 2008. — Vol. 179. — Pp. 685-693.

6. ODPEVP: A Program for Computing Eigenvalues and Eigenfunctions and Their First Derivatives with Respect to the Parameter of the Parametric Self-Adjoined Sturm-Liouville Problem / O. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, A. G. Abrashkevich // Computer Physics Communications. — 2009. — Vol. 180. — Pp. 1358-1375.

7. POTHEA: A Program for Computing Eigenvalues and Eigenfunctions and Their First Derivatives with Respect to the Parameter of the Parametric Self-Adjoined 2D Elliptic Partial Differential Equation / A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, A. G. Abrashkevich // Computer Physics Communications. — 2014. — Vol. 185. — Pp. 2636-2654.

8. Kumar K., Baranger M. Complete Numerical Solution of Bohr's Collective Hamil-tonian // Nuclear Physics A. — 1967. — Vol. 392. — Pp. 608-652.

9. Tetrahedral Symmetry in Nuclei: New Predictions Based on the Collective Model / A. Dobrowolski, A. Gozdz, K. Mazurek, J. Dudek // International Journal of Modern Physics E. — 2011. — Vol. 20. — Pp. 500-506.

10. Стренг Г., Фикс Г. Теория метода конечных элементов. — Москва: Мир, 1977.

11. Бате К., Вилсон Е. Численные методы анализа и метод конечных элементов. — Москва: Стройиздат, 1982.

12. Becker E. B., Carey G. F., Tinsley Oden J. Finite Elements. An Introduction. Vol. I. — New Jersey: Englewood Cliffs, Prentice Hall, 1981.

13. Chuluunbaatar O. The Scientific Doctoral Thesis. — 2010.

14. Pockels F. Uber die partielle differential-gleichung Au + k2u = 0 und deren auftreten in der mathematischen physik. — Leipzig: B. G. Teubner, 1891.

15. Berry M. V., Wilkinson M. Diabolical Points in the Spectra of Triangles // Proceedings of the Royal Society A: Mathematical, Physical & Engineering Sciences. — 1984. — Vol. 392. — Pp. 15-43.

16. Решение краевых задач для систем ОДУ большой размерности: эталонные расчеты в рамках метода Канторовича / А. А. Гусев, О. Чулуунбаатар, С. И. Ви-ницкий, В. Л. Дербов // Вестник РУДН. Серия: Математика. Информатика. Физика. — 2016. — Т. 3. — С. 31-37.

17. Cornwell J. F. Group Theory in Physics. — New York: Academic Press, 1984.

18. A MAPLE Symbolic-Numeric Program for Solving the 2D-Eigenvalue Problem by a Self-Consistent Basis Method / I. N. Belyaeva, N. A. Chekanov, A. A. Gusev, V. A. Rostovtsev, Yu. A. Ukolov, Y. Uwano, S. I. Vinitsky // Lecture Notes in Computer Science. — 2005. — Vol. 3718. — Pp. 32-39.

19. Eisenberg J. M., Greiner W. Nuclear Theory. Vol. 1. — Amsterdam: North-Holland, 1970.

© Gusev A. A., Chuluunbaatar O., Vinitsky S.I., Derbov V.L., Gozdz A., 2017