The algorithms of the numerical solution to the parametric two-dimensional Boundary-Value problem and calculation derivative of solution with respect to the parameter and matrix elements by the finite-element method Текст научной статьи по специальности «Физика»

UDC 517.958:530.145.6 The Algorithms of the Numerical Solution to the Parametric Two-Dimensional Boundary-Value Problem and Calculation Derivative of Solution with Respect to the Parameter and Matrix Elements by the Finite-Element Method

A. A. Gusev

Laboratory of Information Technologies Joint Institute for Nuclear Research 6, Joliot-Curie str., Dubna, Moscow region, 141980, Russia

The effective and stable algorithms for numerical solution with the given accuracy of the parametrical two-dimensional (2D) boundary value problem (BVP)are presented. This BVP formulated for self-adjoined elliptic differential equations with the Dirichlet and/or Neumann type boundary conditions on a finite region of two variables. The original problem is reduced to the parametric homogeneous 1D BVP for a set of ordinary second order differential equations (ODEs). This reduction is implemented by using expansion of the required solution over an appropriate set of orthogonal eigenfunctions of an auxiliary Sturm-Liouville problem by one of the variables. Derivatives with respect to the parameter of eigenvalues and the corresponding vector-eigenfunctions of the reduced problem are determined as solutions of the parametric inhomogeneous 1D BVP. It is obtained by taking a derivative of the reduced problem with respect to the parameter. These problems are solved by the finite-element method with automatical shift of the spectrum. The presented algorithm implemented in Fortran 77 as the POTHEA program calculates with a given accuracy a set ~ 50 of eigenvalues (potential curves), eigenfunctions and their first derivatives with respect to the parameter, and matrix elements that are integrals of the products of eigenfunctions and/or the derivatives of the eigenfunctions with respect to the parameter. The calculated potential curves and matrix elements can be used for forming the variable coefficients matrixes of a system of ODEs which arises in the reduction of the 3D BVP (d = 3) in the framework of a coupled-channel adiabatic approach or the Kantorovich method. The efficiency and stability of the algorithm are demonstrated by numerical analysis of eigensolutions 2D BVP and evaluated matrix elements which apply to solve the 3D BVP for the Schrodinger equation in hyperspherical coordinates describing a Helium atom with zero angular momentum with help of KANTBP program.

Key words and phrases: parametrical two-dimensional boundary value problem, elliptical equation in partial derivatives, finite element method, Kantorovich method, hyper-spherical coordinates, Helium atom.

1. Introduction

Mathematical models of few-body systems in molecular, atomic and nuclear physics, as well as physics of semiconductor nanostructures are described by boundary value problems (BVPs) for the multidimensional equation of Schrodinger type in configuration space Rd, for example: spectral and optical characteristics of excited states of a Helium-like atom [1-16], photoionization and recombination of opposite charged particles (positrons, antiprotons) in the magnet-optical trap [17,18], optical absorption in quantum wells, quantum wires [19], and quantum dots [20], channeling of likely charged particles in thin doped films [21] and resonance tunneling of composite systems through repulsive barriers [22,23].

Efficient and stable algorithms for the numerical solution of such a class of BVPs [24] are based on its reduction to a system of ordinary second-order differential equations (ODEs) with respect to one of the variables by the Kantorovich method (KM) [25].

Received 5th September, 2013. The author thanks Drs. A.G. Abrashkevich and O. Chuluunbaatar, Profs. V.L. Derbov, L.A. Sevastyanov and S.I. Vinitsky for collaboration. The work was supported partially by grants 13-602-02 JINR, 11-01-00523 and 13-01-00668 RFBR..

It is quite natural to use the parametric eigenfunctions of an auxiliary parametric BVP in an appropriate region of configuration space Rd-1 as the basis for the expansion of the unknown (required) solution. It was shown that an efficient application of the KM requires the development of a set of symbolic-numerical algorithms for computing the following quantities to a prescribed accuracy:

1) parametric eigenvalues and eigenfunctions of the auxiliary parametric BVP;

2) derivatives with respect to the parameter of the eigenvalues and eigenfunctions, and matrix elements (integrals of the eigenfunctions multiplied by their derivatives with respect to the parameter);

3) asymptotics in the parameter of the eigenfunctions and of the matrix elements that appear as variable matrix coefficients in the system of ODEs by the variable;

4) asymptotics of the solutions to the system of the ODEs for small and large values of the variable;

5) solutions of BVP for the system of ODEs and numerical analysis of convergence of solutions with an increased number of equations.

The problems 1), 2), and 3) for a parametric 1D BVP were considered in [26-28] while 4) and 5) ones in [29-31]. However, the problems 1), 2), and 3) even for a parametric 2D BVP with weakly singular and long-range potentials of the Coulomb type have been discussed for a long time, for example, in [1-11], but have not been exhaustively solved till nowadays [12-14].

In this paper, we present the effective and stable algorithms of solving the problems 1) and 2) for the parametric 2D BVP for elliptic partial differential equations with boundary conditions of Dirichlet and/or Neumann type in a finite two-dimensional region which arise in the reduction of the 3D BVP (d = 3) to a system of ODEs using the KM. We seek the solution of the parametric 2DBVP in the form of expansion in the basis eigenfunctions of the auxiliary Sturm-Liouville problem for the ODE with respect to one of the variables. The coefficients of the expansion are the parametric vector-functions which are eigensolutions of the parametric homogeneous 1D BVP for a system of ODEs obtained by averaging the original problem over the basis eigen-functions. The required parametric derivatives of eigenvalues and the corresponding vector-eigenfunctions are calculated as a solution of the parametric inhomogeneous 1D BVP which is obtained by taking a derivative of the above parametric homogeneous 1D BVP with respect to the parameter.

The main goal of this paper is the following. The finite element method (FEM) [32, 33] is used to discretize both the homogenous and nonhomogeneous parametric problems and to construct numerical schemes for calculation of eigenfunctions and their derivatives with respect to the parameter within accuracy of the same order of 0(hp+1), and also the eigenvalues and the matrix elements within accuracy of the same order of 0(h2p) in step h of a finite element grid. Choice of the order of approximation p depends on the smoothness of the desired solutions. An auxiliary algorithm to calculate the lower bound of the smallest eigenvalue of the parametric generalized algebraic eigenvalue problem is presented and tested. This algorithm allows one to determine automatical shift of the spectrum and to save iterations in the calculation of a set of eigensolutions by the iterations in a subspace using the routine SSPACE [33].

The efficiency of the algorithm is demonstrated by numerical analysis of solutions of the parametric 2D BVP including evaluation of matrix elements applied to reduce the 3D BVP for the Schrodinger equation of a Helium atom with zero total angular momentum in the body-fixed hyperspherical coordinates to the 1D BVP for the ODEs with respect to the hyperradius by the Kantorovich method [10]. The corresponding benchmark calculations with a given accuracy of the ground state energy and first exited state energy of a Helium atom and their convergence versus both the number of the basis vector-eigenfunctions and the number of their components are presented.

The structure of the paper is the following. In Section 2, we present the statement of the problem and reduction of the parametric 2D BVP to the parametric 1D BVPs. In Section 3, the formulation of the parametric algebraic eigenvalue problems by the FEM is given. In Section 4, the benchmark calculations of matrix elements and the

eigenenergies of a Helium atom are presented. In Section 5, we discuss the results and perspectives.

2. Statement of the Problem

Let us consider the parametric 2D BVP for a self-adjoined second order elliptic differential equation in a two-dimensional region Qx,y = (^min,£max) x (ymin,ymax)

/1 d d 11 d d \

1 My)£ -T^T^ih(x)lL + U(^ Z) - £i(z)j Bi(x,y; z) = 0 (1)

(2)

fi(y) dy dy fs(y) f4(x) dx dx with the Dirichlet and/or Neumann type boundary conditions

lim f2(y) dBt(^,y' z) =0 or Bi(x, ymin', z) = 0,X G (^min,^max), y^ymin ay

lim f2(y) dBt(^,y' z) =0 or Bi(x, ^max; Z) = 0, X G (^min,^max),

J^^max Oy

lim /5 (x) z) =0 or Bi(Xmin, y; z) = 0,y G [ymin, ymax],

% ^^min 0X

lim ¡5(x) dBt(^,V; z) =0 or Bi(Xmax, y; z) = 0, y G [ymin, ^max].

Here z G Qz = [^min,^max] is a real-value parameter; the functions fi(y) > 0, f2(y) > 0 /a(y) > ° fi(x) > 0, ¡5(x) > ° and dy/2(y), дxf5(x), U(x,y; z), dzU(x,y; z) are continuous on the (x,y) G Qx,y. Also assume that the parametric BVP (1), (2) has only a discrete spectrum.

The main goal is realized by the following algorithm.

In Step 1 algorithm calculates a set of jmax smallest eigenvalues e1(z) < e2(z) <

IN

... <en(z), and £i(z) > a(z), and the corresponding eigenfunctions {Bj(x,y; z)}j=1 G ,y)

^max

■H],

Fz ~ L2(na,y), satisfying the orthogonality and normalization conditions

x ^max

dVfi (y) J dxf4(x)Bi(x,y; z)Bj(x,y; z) = ^, (3)

where Sij is the Kronecker symbol and a(z) > —to is the lower bound of the smallest eigenvalue of £1(z).

In Step 2 algorithm computes a set of partial derivatives of eigenvalue dej (z)/dz and partial derivatives of eigenfunctions dBj (x, y; z)/dz with an accuracy of the same orders achieved for eigenvalues and eigenfunctions of the BVP (1)-(3), respectively.

In Step 3 algorithm computes matrix elements defined by the integrals

dBi(x,y; z) dBj(x,y; z)

nj w = ±±ji w = 1 dyjiyy) 1 --Q~z---Q~z-,

Hij(z) = Hji(z) = I dyfi(y) j dxf4(x)~ ^ ^

Kmin, (4)

^min

J/max ^max

Ql3 (z) = —Q3l(z) = —J dyfi(y) J dxf4(x)Bt(x,y; z)8Bj ^ Z),

ymin ^min

with an accuracy of the same order achieved for the corresponding eigenvalues of the BVP (1)-(3).

y

X

The calculated eigenvalues £i(z) and matrix elements Hij(z), Qij(z) can be used for solving the bound state and multichannel scattering problems for a system of coupled ODEs with respect to the variable z with the help of the KANTBP programs [29-31].

2.1. Reduction of the Parametric 2D BVP to the Parametric 1D BVPs

The partial wave function Bi(x, y; z) is expanded over the orthonormal basis functions {^j (jmax ^ to) in the conventional (C) form

J max

B%(Xly; z) = £ ^j (x)éÎ\y; z), (5)

J=1

or {^j (x; y, z)}j=a1x (jmax ^ to) in the Kantorovich (K) form

max

Bi(x^ y; z)=Y^/ ^J(x; z)éj )(y; z). (6)

=1

T

In Eq. (5), the vector-functions £(i )(y; z) = (^f ) (y; z),..., ¡¿¡1^(y; z^J = &(y; z) =

(£h(y; z),... ,îjmaxi(y; z)) are unknown. The functions (x) are determined as solutions of the following eigenvalue problem in the C form:

( ■ wiè/5(a°è + Uo(x)) *j(x) = aj*j(x),

lim M®)^^ =0 or 4J (®min) = 0, (7)

■ ^ xmin dX

X-S? f5 (,)^ = 0 or jj Ow)^ where U0(x) is a known function and

dx f4(x)^i(x)^j(x) = 5ij. (8)

The functions (x; y,z) are determined as solutions of the following eigenvalue problem in the K form:

^ 1 d+ U(x; y,z^j ^j(x; y,z) = Xj(y; z)^j(x; y,z),

f4(x)dx dx

lim f5(x)d^J fo y,z) =0 or 4J (^min; y,z) = 0, (9)

x ^xmin dX

lim f5(x)d^ (XA; V,Z) =0 or ^ (w; y,z) = 0,

% ^^max dX

where U(x,y; z) is a known function of the original problem (1) and

dxfi(x)ipi(x; y,z)^j(x; y,z) = % (10)

X

Note, these problems can be numerically solved with a given accuracy by mens of the ODPEVP program [28].

After minimizing the Rayleigh—Ritz variational functional, and using the expansion (5), the parametric BVP (1)—(3) is reduced to a finite set of jmax ODEs the C or K forms

1 9 9

(D(y;z) - ei(z) I) t(i)(y;z) = 0, D(y;z) = -j-^.+ W(y;z), (11) lim /2 (y)9 ^\(y]Z) =0 or £(i)( ymin;Z) = 0,

V^Vmin By

lim f2(y)9^() (y]Z) =0 or Ç(i)(ymax;z) = 0,

V ^Vmax By

ymax

hj = % = J dy h (y) [n(i)(y;z))T t(j)(y;z),

ymin

(12)

(13)

Here I and W(y; z) are the symmetric matrices of dimension jmax x jmax in the C form

7ij(y;z) = 2/+( y) bij + / dx Î4(x)^i(x){^J ^^ y;z) - jg^ I ^J

ipj (x)

or W(y; z) = W(y, 9y ;z) are the self-adjoint matrix differential operator of dimension Jmax X jmax in the K form

Wij(y;z)

\i(y;z) + \j(y;z) + hMiw /ww 2fs(y) diJ + fi(y)Hij(V;Z)+ 1 / 9

+fy<^)]J+mq"I- (14)

Hij (y;z) = Hij (y;z)= /4 (x)

B^i(x; y, z) B^j(x; y, z)

9y

By

Qij(y;z) = -Qji (y; z) = - Î4(x) ^i(x; y, z)

(x; y, z) By '

Taking a derivative of the boundary problem (11)-(13) with respect to parameter

z, we get that Bz£(i\y;z) can be obtained as a solution of the following parametric inhomogeneous BVP:

(D(y;z) - ei ( z) I)

B^(i) (y;z) Bz

B

— (W(y;z) - £i(z) I)

£(i \y;z), (15)

lim M ) b2£(i)(y;z) = 0 or BÇ(i)(ymin;^) =0

V ^Vmin BzBy Bz '

.. ) 92t(i)(y;z) 0

lim /2 ( y)-jr^,-= 0 or

y ^Vmax BzBy

9£(i) (ymax; z) B

X

X

X

The parametric BVP (15), (16) has a unique solution if and only if the conditions are fulfilled

ymax /.x

/dtfife) («<%;,■))" ^(^=0, (17)

y

min

ymax

dz j dyh{y) (¿«(y^f (18)

y

min

In this case the required matrix elements (4) are represented by the integrals

ymax / / .x \ T / -N

Hij(z) = Hji (z) — J dy f1(y) I -9^- I

9

ymin

ymax

t 9fiu)

(19)

9t(J)(y;z)

Qij(z)- -Qji(z)- - dyfr(y) ^(i)(yz)) ^

y

2.2. Continuity Conditions for the Eigenfunction Bi(x, yz)

Since the problems (1)-(3) and (11)-(13) are homogeneous, it is necessary to use an additional condition to support the continuity of the vector-functions £(i)(y z) and matrix elements (19) with respect to the parameter z on the interval Qz — [zmin, zmaxj. We used the following additional condition:

1) At the first point z — z\ G Qz, find a value y — y0, in which the eigenfunction Bi(x0, yo', zi) reached absolute maximum value and fix the sign of the value of the eigenfunction Bi(x0, y0, z{). Here x0 G [xm;n, xmax] is a fixed point and at least one of the functions (x0) in expansion (5) is not equal zero.

2) At the next points z G Qz compute the value of the eigenfunction Bi(x0, y0, z) and compare its sign with the sign of the previous one. If they are different, change the sign of Bi(x0, y0, z) and again find a new value y — y0, in which the eigenfunction Bi(x0, y0, z\) reached absolute maximum value, and fix the sign of the value of the eigenfunction Bi(x0, y0',z).

Note that if the grids of Qz are dense, the above algorithm works well.

3. Formulation of Parametric Algebraic Eigenvalue Problems by the FEM

Let us consider a numerical algorithm for the calculation of the vector-eigenfunctions )(y;z) of the parametric boundary problem (11)-(13) and their derivative with respect to the parameter z. Computational schemes of the high order of accuracy are derived from the Rayleigh-Ritz variational functional

ymax I Jmax / , , \ 2 Jmax \

m m:/2(y)(+ ma h(y)^(yz)w^(y,*)&

£) = ymi—-y-^-'—, (20)

ymax Jmax

I ma A(y)e( y;z)dy

ymin ^=1

on the basis of the FEM. The general idea of the FEM in one-dimensional space is to divide interval [ymin, J/max] into many small domains called elements. The size of

elements can be defined very freely so that physical properties or quality behavior of solutions can be taken into account.

The interval A = [ymin, ymax] is covered by a system of n subintervals Aj =

[yj-1, yj] in such a way that A = (J"=1 Aj. In each subinterval A j the nodes

yPj,r = Vj-1 + "f ^ hj = yj — yj-u r = 0,p,

(21)

and the Lagrange elements {<ßp r(y)}P=t

r=0

tir (y)= n ^

i=0, i=r j

( yp,r yj,i)

(22)

are determined. By means of the Lagrange elements ^(y), we define a set of local functions Ni(y) as follows:

<Pp,o(У), 0,

ti.r (y),

0,

y G Al,

y G Ai, y G Aj,

y G Aj,

= 0,

I = r + p(j — 1), r = 1,p — 1,

Nj(y) =

^j,j(y), y G Aj,

(Pp+l,o(У), yG Aj+1,

0, y G A j U Aj+1,

Pip (y), y G

0, y G An,

(23)

l = jP, j = 1,n — 1, I = np.

The functions {Nf(y)}= o, L = np, form a basis in the space of polynomials of the p-th order. Now, each component of the vector-functions £(y; z) e approximated by a finite sum of local functions Nf(y)

6 ( y; z) = Ê ^ ( y;z)NP(y),

(24)

1=0

i.e. the vector-function £(y; z) has a generalized first-order partial derivative and belongs to the Sobolev space H1(Qhy) [32].

After substituting expansion (24) into the variational functional (20) and minimizing it [32,33] we obtain the generalized eigenvalue problem

AP^h = £ hBP£h. (25)

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

n n

AP = £ BP = £ ^ (26)

3=1 3=1

where the local matrices ap- and b^ are calculated as

=/{vh( y)4 d1;(V) d ^ + ( y, ( y)^r(y)}^dry,

-1

+1

(b^ = S^J /i( y)LPPj,q(y)^lr(y)ydr?, (27)

-i

y = yj-i + 0.5hj(1 + ^), q,r = 0,p. Integrals (27) are evaluated using the Gaussian quadrature formulae

Kr =£{ V h( y9) h2 +/i( y9 )W,V ( y,, z)pj>q ( y, )p*p ( y9 )} hj^p,

p h ■

(bP)qr = <V £ /l(yp)p£q(yp)p£r(yg), (28)

9=0

where yg = yj-1 + 0.5hj(1 + %), r]g and wg, g = 0,p are the Gaussian nodes and weights related to the orthogonal polynomial of order p + 1. Note in this approach maximum value of a half-band of matrices Ap and Bp is small compared to their dimension and is not greater than jmax x (p + 1).

Let D( y, z) be a continuous and bounded positively defined operator on the space Vi with energy norm, £i(z), y\z) £ V2 are the exact solutions of (11)-(13), and £h(z), £h(y;z) £ V1 are the corresponding numerical solutions. Then the following estimates are valid [32]:

Iej( y) - eH < cih2^, ||£«( y,z) - £«h||o < c2h^1, ci > 0, c2 > 0, (29)

where \\£«( y,z)\\0 = f^ dy fi(y)(^l)( y,z))T^(y,z), h is the grid step, p is the order of finite elements, j is the number of the corresponding eigensolution, and the constants i and 2 do not depend on step h. It is necessary to mention that the second estimate of Eq. (29) is valid also for solution d£,(%\y,z)/dz of problem (15)-(17). This fact guarantees the same accuracy for eigenfunctions and their derivatives within the present method.

In order to solve the generalized eigenvalue problem (25), the subspace iteration method [32,33] elaborated by Bathe [33] for the solution of large symmetric banded matrix eigenvalue problems has been chosen. This method uses a skyline storage mode which stores 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 [33]). 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 matrix Ap in Eq. (25) is not positively defined, problem (25) is replaced by the following problem:

Ap <ph = eh Bp £h, Ap = Ap - aBp. (30)

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

3.1. Calculations of Parametric Derivative of the Eigenfunctions and

Matrix Elements

The boundary problem (15)-(17) is reduced to the linear system of inhomogeneous algebraic equations

9£h u 9 £h id A" deh \

Lir^(A"-£hB"ir = * »=-(^- *BP)

- ^rBPI ih. (31)

The normalization (13) , orthogonalization (17) and additional conditions (18) read as

Bpi" = 1. (f)TBpih = 0, |h = (£")T ^ (32)

From here, the potential matrix elements H^ (z) and Qij (z) have the form

T

(z;=i dèl BP^. Q'i(z) = - (ihf BP%. (33)

H"(z) = (§) B"f, Q"(z) = - («h):

Since eh is an eigenvalue of (25), matrix L in Eq. (31) is degenerate. In this case the algorithm for solving Eq. (31) can be written in three steps as follows:

Step k1. Calculate solutions v and w of the auxiliary inhomogeneous systems of algebraic equations

Lv = b, Lw = d, (34)

with non-degenerate matrix L and right-hand sides b and d r = J LSS', (s -S)(s' ) = 0,

Lss' = I Saa', (s -S)(s' -S) = 0, (35)

T = J bs, s = S, = f LsS, s = S, bs = \ 0, s = S, ds = \ 0, s = S,

where S is the number of the greatest absolute value element of vector Bp^ Step k2. Evaluate coefficient 7

h

7 = -m 71 ,, 71 = vT B"Çh, 72 = wT B"Çh, Ds = (B"t)s. (36) (Ds - 72)

Step k3. Evaluate vector 9Çh/9z

9îî" ( vs - 7ws, s = S,

i vs

7,

9z 1 7, s = S.

(37)

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

Theorem 1. Let D( y, z) be a continuous and bounded positively defined operator on the .space H1 with energy norm. Also dzW(y; z) is continuous and bounded in each value of the parameter z. Then for exact values of solutions e¿(z), dz£i (z), Hij (z), Qij(z); £i(y;z) G H2, dz(y;z) G H2 and the corresponding numerical values eh(z), dzeh(z), Hh(z), Qh(z); £,i(y;z) G H1, dz£i(y;z) G H1 the following estimates are

valid:

de i(z) de ¿(z)

dz

dz

< czh2p,

dii (y;z) dg(y,z)

d

d

< ahp+\

Qn(z) - Qïj(z)\ < c5h2p, \Hij(z) - Hi>Az)\ < c6h2p,

(38)

(39)

where h is the maximal step of the finite-element grid, p is the order of finite-elements, i, j are the number of the corresponding solutions, and constants c3, c4, c5 and c6 do not depend on step h.

Proof is straightforward following the proof scheme, in accordance with [32,34].

0

3.2. Finding the Lower Bound for the Lowest Eigenvalue of the Generalized Eigenvalue Problem

In the general case, the boundary value problem (11)-(13) has both negative and positive eigenvalues because we did not require the positive definiteness of matrix W( y, z) providing non-negative eigenvalues of 0 < ei( z) < e 2( z) < ... < e N (z). Because of this fact, the stiffness matrix Ap = Ap(z) in the general case is not positive defined.

Therefore, it is difficult to get a close enough lower estimate a = a(z) < eh(z) of the smallest eigenvalue ei(z) for the problem (25), for example, using the known theoretical estimates of Gershgorin disks [35], for each real value of the parameter z from a given interval Qz. Based on the known properties of the decomposition of F = LDLT for symmetric matrices F: for symmetric positive defined matrix F all elements of the diagonal of the matrix D are positive defined [36], we proposed and tested the following algorithm:

Step 1. Calculate LDLT factorization of Ap — aBp.

Step 2. If some elements of the diagonal matrix D are less than zero then put a = a - 1 and go to Step 3, else go to Step 5.

Step 3. Calculate LDLT factorization of Ap - aBp.

Step 4. If some elements of the diagonal matrix D are less than zero then put a = a - 1 and go to Step 3, else put a = a - 0.5 and go to Step 8.

Step 5. Put a = a + 1 and calculate LDLT factorization of Ap - aBp.

Step 6. If all elements of the diagonal matrix D are greater than zero, then repeat Step 5 i.

Step 7. Put a = a - 1.5.

Step 8. End.

After using the above algorithm one should find the lower bound for the lowest eigenvalue, and always ei(z) - a < 1.5.

The initial approximation a(z) = a0 in the algorithm (1—8) is given a priori, for example, a0 = 0, or with known asymptotic estimates of a0 < i( ).

In carrying out this algorithm (1—8), we obtain a lower estimation of the smallest eigenvalue of (25), and always ei(z) - a* < 1.5, i.e. the resulting matrix Ap - a*Bp is positive defined. As a result, (25) reduces to

(Ap - a*Bp) ^h = eh Bp^h, eh = eh + a*, (40)

which has positive eigenvalues 0 < eh( z) < ... < eN(z).

Thus, the algorithm (1—8) allows one to automatically determine the shift a* of the spectrum of problem (25) and save the number of iterations in the calculation of the solution set by the iteration in the subspace using the routine SSPACE [33] for any value of the parameter z of a given interval Qz. To solve the sequence of BVPs

xIn this step misprint of paper [26,28] is corrected

(11)-(13) on the grid Qhi [Zmin, ¿max] = {Z! = Zmin, Zi+1 = Zi + hi, Zimax = ¿max} with a sufficiently small step variable hi, the values of a(zi+1) can be estimated by the formula a(zi+1) = e'h(zi) + (delh(zi)/dz — \delh(zi)/dz| - 2f5i)hi/2, pi > 1, calculated in the previous step zi, not using an algorithm (1-8). In this case, the algorithm (1-8) is used, for example, only when the value of z = zmin or z = zmax.

The above algorithms have been implemented in Fortran 77 as a program complex POTHEA [31].

4. Benchmark Calculation of the Eigenenergies of a Helium

Atom

Time-independent Schrodinger equation for a helium atom with zero total angular momentum in the body-fixed hyperspherical coordinates z = R £ [0, y = a £

[0, n], x = •& £ [0, n] can be written as a BVP for the following 3D-elliptic equation in atomic units [10]:

(

1 d . d 4 \

-R dRR5 dR + R (Ha] R) + ya;R)) + 2E\ *(R, a, #) = 0, (41)

H(â, a; R) =--2— TT" sin a——+ ~

sin2 a

1 ( d . 2 d 1 d .

[dasin ada + sînMntsin(â)w),

V(â,a;R) = R (----+ 1 1

V 7 2 \ sin(a/2) cos(a/2) - sin( a) cosâ)

The total wave function R,a,â) satisfies the following boundary conditions:

(42)

lim R5 9^(B:,'a, ^ = 0, lim R5^(R, a, â) = 0, R^o dR R^C

. 2 d^(R,a,â) . fa.dV(R,a,â) n

lim sin a---= 0, lim sin(â)-—-= 0,

a^o,-K da dâ

and is normalized by the condition

CO K k

J d RR5 y d a sin2 a J dâ sin â^2 (R,a,â) = 1. (43)

4.1. Reduction of the 3D BVP to the 1D BVP: Kantorovich Expansion

Consider formal expansion of the solution of Eqs. (41)-(43) over a set of two-dimensional parametric basis functions {Bi(a,^;R)}^=1 (N ^ to):

*(R,a,#) = ^2Bj (^,a;R)X3 (R). (44)

3 = 1

T

In Eq. (44) the functions x(R) = (%1( R),x2(R),... ,Xn(R)) are unknown, and the

adiabatic functions B( a; R) = (B1(a, R), B2(§, a; R),..., Bn($, a; R)) form an orthonormal basis for each value of hyperradius R which is treated here as a slowly varying adiabatic parameter.

After minimizing the Rayleigh-Ritz variational functional, and using the expansion (44), 3D BVP (41)-(42) is reduced to a finite set of N ordinary second-order differential

equations for x(R)

(

- R1 dRR5 dR + U(R) + Q< + R- dRdRR) - 1) *(*) = ».

lim R5 = 0, lim R5x(R) = 0.

R^o dR R^^

(45)

Here I, U( R) and Q(R) are the matrices of dimension N x N :

lij = Sij. Uij(R) = Uji(R) = 2 '( )R+2 j( )àij + Hi:j(R)

dBi(à,a;R) dBj(à,a;R)

Hij (R) = H^ (R) = da sin2 a d'à sin à

dR dR (46)

oo

K K

Qij(R) = ~Q,i (R) = -I da sin2 a I dà sin àBi(à, a; R)dBj ^R^ •

This problem can be solved by FEM at values R belonging to the Gauss nodes of a finite element grid Q R with help of KANTBP programs [29-31]. In the KM [24,25] the parametric basis functions Bi(°&, a; R) are determined as solutions of the following parametric BVP:

d d 1 d d \ -sin2 + ¡1^MsinV) +va; R) -

Bi(ê,a; R) = 0, (47)

. 2 dBi(à,a;R) n . dBi(à, a; R) lim sin a---= 0, lim sin à-—-= 0,

a^-o,' da dà

' ' (48)

/ da sin2 a I dà sin àBi (à, a; R)Bj (à,a;R) = %

oo

1

4.2. Reduction of the Parametric 2D BVP (47)-(48) to the Parametric

1D BVP

Consider the following expansion of adiabatic surface function in the C form: Bi(a, à; R):

J max

Bi(à,a;R) = (à) ef(a;R). (49)

3=1

Here ipj (à) = Pj-1(cosà) is the normalized Legendre polynomial:

sin à^^ = A^-(à), A, = j(j - 1), sin à d à d à

lim sin( à)d^Î(à) =0, / dà sin àipi(à)ipj (à) = tf^-o,' dà J

o

(50)

After minimization of the variational functional we get that the eigenfunctions

T

)(a;R) = R), ^(a; R),..., ^j^^(a;R)^ and the eigenvalues £i(R) satisfy

the following eigenvalue problem for a set of jmax ordinary differential equations:

/18 8 \

--2 Î78 sin2 a— + W(a; R) — e^R) Î £(i)(a; R) — 0,

V sin2 a 8a 8a )

.. . 2 8^(a;R) lim s.n a---

a^0,TT 8a

0.

Here I, W(a; R) are the symmetric matrices of dimension jmax x jmax

K

Sij = j da sin2 a (^(i) (a; R)) (a; R)

hj —

Wij(a; R)

Xi + Xj . R

+ ^ -

2

2

_ 2 sin2 a 2 V sin(a/2) cos(a/2)

R

% + R2^i/p(a),

W/?ep(a) — d' sin'

. Pi-1 (cos')P0-i(cos') _ [A Pi-i(n)Pj-i(ri)

V1 — sin a cos'

d

1

1/1 — (sin a)r]

(51)

(52)

(53)

Because of the symmetry of matrix elements, Wij(a; R) — Wij(i — a;R), with respect to a — i/2, problem (51)-(53) will be considered for a G [0,i/2] with the following boundary conditions for the ground and first exited states:

, 2 8tj(i)(a;R) 0 lim sin a---— 0.

8 a

The 1D weak singular integral (53) was conventionally calculated analytically using the Clebsch-Gordan coefficients [8,9]. But this approach gives big numerical errors at large numbers and because of calculations of the factorial of large numbers (required factorials of numbers up to 4jmax — 3). After the change of the variable in (53)

V ^ V(a, 0 — — C2) + C, CG [—1,1], a G [0,i/2],

we obtain the nonsingular integral, WirJep(a) — WirJep(i — a)

1

1

W£p(a) = cos-1 j dCPi-1 (v(a, 0) P-i (v(a, Q).

J cos( a/2) J

-1

The last ID integral calculated by means of the 96-order Gauss-Legendre quadrature and this approach give results with accuracy < 10-14 at i, j < 50, i.e. with the double precision accuracy.

The potential curves 4 R-2(ei(R) + 1), radial diagonal and nondiagonal matrix elements Hij (R), and radial matrix elements Qij (R) as functions of hyperradius R are displayed in Figure 1. As can be seen from Figure 1, our algorithm for continuity conditions from Sec. 2.2 for the eigenfunction Bi(a,rd;R works well.

The numerical experiments show a strict correspondence with the theoretical estimations for the eigenvalues, eigenfunctions of (29) and their derivatives with respect to the parameter (38). In particular, we calculated the values of the Runge coefficients

-H,:

--- H1: - — H„

-— H2:

........... H2.

------H

-0,3

6 8 R

10

12

14

1,0 0,8 0,6 0,4 0,2 0,0

0,4

0,20,0 -0,2-0,4-0,6-0,8-1,0

-H„

i j ■ - - - H22

! i — - H33

i ; > /i ; > / .....H44 -- H55

.'•' 1. / / 'À ;: i : \

68 R

10

12

14

68 R

Figure 1. Potential curves 4R-2(£i(R) + 1) (top-left), radial diagonal (top-right) and nondiagonal (bottom-left) matrix elements Hij(R), and radial matrix elements Qij (R) (bottom-right) plotted vs hyperradius R

2

4

0

2

4

Pi = l°g2

ai h/2 ~al

af'2 h/4

1=1 + 6.

(54)

with absolute errors on four twice condensed grids for their eigenvalues, their derivatives, and matrix elements, respectively:

h \h/8( n h\ h

= 1 £j (P) - £jh =

ah = \\C>;R) - tj(a-,R)\\o, ah =

de11/8 (R) de'

dR

„h/8

_¿_

dR

dC'8 (a;R) dÇh(a;R)

d R

d R

(55)

aj = lffjJ/8(R) - Hj(R)|, aj = |Qhf(R) - Qj(R)|,

From (54) we obtained numerical estimations of the convergence order of proposed numerical schemes, i.e. theoretical estimations equal to Pi = p+1 for 1 = 3,4 and Pi = 2p otherwise. For the chosen approximation order p = 4 for their eigenvalues, their derivatives, and matrix elements we obtain numerical estimations of Runge coefficients within 7.5^7.8, and for their eigenfunctions and their derivatives in the range 4.6+4.8, which corresponds to the theoretical error estimates at fixed number jmax of equations (51). The calculations (55) were performed with a specified accuracy of ~ 10-12 by means of POTHEA program at relative error tolerance e\ = 4 ■ 10-16 of calculated eigenvalues (30) on computer 2 x Xeon 3.2 GHz, 4 GB RAM, using the Intel Fortran

0

77 and the data type of real*8, which provides 16 significant digits. The running time for this example is 2 seconds for jmax = 12, N = 6 and 1000 seconds for jmax = 50, N = 50.

A convergence study of several matrix elements with respect to number max = 12, 28,40, 50, 60, 70,80,100,120 of the Legendre polynomials and the number of finite elements, Net = 6,12,18,24, 30, 36, of the grid Qa = {0(Ne^/2} and their order p = 4 are presented in Tables 1-4. One can see that the potential curves 2 R-2(sj(R) + 1) and the matrix elements Hij (R) converge monotonically from above, with the increasing of the numbers Nei and jmax. The absolute values of the matrix elements Qij (R) converge monotonically from above with increasing jmax and from below with increasing Nei.

Table 1

Convergence of potential curves 2R-2(£3(R) + 1), 2R-2(£5(R) + 1) at R = 7.65 as a function of maximum number of Legendre polynomials and numbers of finite elements Nei of order p = 4

J max Nei 2 R-2(e 3 (R) + 1), 10-1 2 R-2 (e 5(R) + 1), 10-1

12 6 —6.179 614 071 100 9 —3.717 091 841 269 3

12 12 —6.179 614 958 562 6 —3.717 092 590 328 6

12 18 —6.179 614 963 052 6 —3.717 092 593 788 5

12 24 —6.179 614 963 229 4 —3.717 092 593 920 9

12 30 —6.179 614 963 246 2 —3.717 092 593 932 9

12 36 —6.179 614 963 248 3 —3.717 092 593 934 6

28 36 —6.179 640 623 131 8 —3.717 093 236 000 2

40 36 —6.179 641 912 441 6 —3.717 093 288 608 1

50 36 —6.179 642 233 040 9 —3.717 093 302 610 3

60 36 —6.179 642 373 030 0 —3.717 093 308 887 9

70 36 —6.179 642 443 638 6 —3.717 093 312 101 5

80 36 —6.179 642 483 034 8 —3.717 093 313 911 2

100 36 —6.179 642 521 803 3 —3.717 093 315 706 2

120 36 —6.179 642 538 813 0 —3.717 093 316 498 0

As it is shown in Tables 1-4, the convergence of eigenvalues and matrix elements vs the number of Legendre polynomials Pj-1 (rj), is proportional to their order ~ j-3. It follows from estimations of values of the matrix elements Wi^j (a) ~ 1/V7 (in particular, for integral (53) at = 1, we have

W1rJep(a) = 2exp(-(j - 1/2)arch(sin-1 a))/(^2j — 1Vsina),

see, for example, [37]) and AJ-1 = (j — 1)j ~ j2, which leads to estimations for the correction of eigenvalues 5e ~ j-3 in second-order perturbation theory. It means that accuracy of an order of ~ 10-12 of calculations will be achieved at least jmax ~ 1500.

In the benchmark calculations the grids in R and a have been chosen as follows: Qr = {0(200)10(200)30} and Qa = {0(150)^/2}. Enclosed in parentheses are the numbers of finite elements of the order = 4 in each interval. The set of matrix elements including the eigenfunction with number N = 50 were calculated with an accuracy of an order of 10-8, using the number of finite elements Nei = 150 at e2 = 10-12. A banded system of (150*4+1)*50=30050 linear algebraic equations (25) with the mean bandwidth (4+1)*50=250 has been stable solved with relative error tolerance e 2 = 10-12 at each value of hyperradius R belonging to the Gauss nodes of the grid Qr. A convergence study of the ground and first exited state energies of a Helium

Table 2

Convergence of matrix elements Q35(R), H35(R), H55(R) at R = 7.65, like Table 1

J max Nel Q35(R), 10-1 H35(R), 10-2 H55 (R), 10-2

12 6 1.345 981 051 4077 2.269 823 448 371 8.333 864 894 452

12 12 1.345 984 463 2035 2.269 832 028 222 8.333 904 465 350

12 18 1.345 984 480 2338 2.269 832 074 795 8.333 904 678 183

12 24 1.345 984 480 8970 2.269 832 076 647 8.333 904 686 599

12 30 1.345 984 480 9591 2.269 832 076 831 8.333 904 687 440

12 36 1.345 984 480 9704 2.269 832 076 853 8.333 904 687 547

28 36 1.345 970 507 6275 2.270 529 070 029 8.335 180 574 435

40 36 1.345 970 198 8424 2.270 563 372 554 8.335 245 600 693

50 36 1.345 970 135 2041 2.270 571 881 381 8.335 261 797 875

60 36 1.345 970 109 7007 2.270 575 593 492 8.335 268 874 497

70 36 1.345 970 097 4981 2.270 577 464 966 8.335 272 444 815

80 36 1.345 970 090 9213 2.270 578 508 889 8.335 274 437 152

100 36 1.345 970 084 6454 2.270 579 535 991 8.335 276 397 958

120 36 1.345 970 081 9686 2.270 579 986 560 8.335 277 258 259

Table 3

Convergence of potential curves 2R-2(£43 (Д) + 1), 2R-2(s4b(K) + 1), like Table 1

J max Net 2 R-2(e 4s(R) + 1) 2 R-2 (e 45 ( R) + 1)

28 6 4.439 005 647 8840 4.879 922 636 3814

28 12 4.438 991 442 6147 4.878 939 387 2213

28 18 4.438 991 281 3574 4.878 936 678 0110

28 24 4.438 991 270 2793 4.878 936 575 2142

28 30 4.438 991 268 8569 4.878 936 565 5674

28 36 4.438 991 268 5908 4.878 936 564 0653

40 36 4.438 814 541 3791 4.878 929 789 5129

50 36 4.438 775 794 7833 4.878 928 117 4560

60 36 4.438 759 636 3135 4.878 927 388 1689

70 36 4.438 751 689 2542 4.878 927 020 3953

80 36 4.438 747 322 9410 4.878 926 815 1861

100 36 4.438 743 080 9327 4.878 926 613 2176

120 36 4.438 741 240 5373 4.878 926 524 5981

atom with the number N of radial equations (45) and the number jmax of the Legendre polynomials are presented in Tables 5 and 6.

One can see that the energy eigenvalues converge monotonically from above, with the N = 45, jmax = 50 - channel value being E1 = -2.903 72415 a.u. and E2 = -2.145 973 22 a.u. Tables show that the obtained results agree with an accuracy of an order of 10-6 at jmax ~ N with variational estimations [15,16] and have a more high accuracy than the previous ones [9,10]. So a similar accuracy can be achieved also in calculations of high exited states of a He atom, for which variational calculations were not applied, taking into account appropriate asymptotic behaviors of the matrix elements and solutions [9].

Table 4

Convergence of matrix elements Q4345(R), H4345(R), H4545(R) at R = 7.65, like

Table 1

J max Nel Q4345(R), 10-3 #4345 (R), 10-4 #4545 (R), 10 3

28 6 7.163 551 693 508 1.313 245 172 874 1.034 074 714 010

28 12 7.192 416 552 701 1.313 393 326 976 1.037 535 372 894

28 18 7.192 470 461 759 1.313 394 061 373 1.037 544 063 945

28 24 7.192 471 802 131 1.313 393 807 683 1.037 544 380 393

28 30 7.192 471 876 618 1.313 393 761 626 1.037 544 409 101

28 36 7.192 471 882 307 1.313 393 751 746 1.037 544 413 457

40 36 7.164 925 249 674 1.304 767 510 954 1.036 946 196 503

50 36 7.158 600 336 853 1.302 825 852 351 1.036 806 107 806

60 36 7.155 920 393 086 1.302 012 009 221 1.036 746 299 243

70 36 7.154 591 453 293 1.301 611 523 137 1.036 716 518 763

80 36 7.153 857 856 101 1.301 391 714 358 1.036 700 039 417

100 36 7.153 142 578 776 1.301 178 645 274 1.036 683 940 956

120 36 7.152 831 407 895 1.301 086 493 836 1.036 676 926 529

Table 5

Convergence of the ground state energy (in a.u.) for a Helium atom versus number N of basis functions and number jmax of the Legendre polynomials

N jmax = 12 [10] Jmax 12 Jmax 21 jmax — 28

1 -2.887 911 68 -2.895 539 19 -2.895 551 19 -2.895 552 76

2 -2.891 379 91 -2.898 631 57 -2.898 643 21 -2.898 644 74

6 -2.903 004 48 -2.903 644 06 -2.903 655 96 -2.903 657 52

10 -2.903 636 13 -2.903 702 86 -2.903 714 79 -2.903 716 36

15 -2.903 705 49 -2.903 708 67 -2.903 720 60 -2.903 722 16

21 -2.903 722 64 -2.903 722 99

28 -2.903 722 66

N jmax — 35 Jmax 40 Jmax 45 Jmax = 50

1 -2.895 553 32 -2.895 553 52 -2.895 553 63 -2.895 553 71

2 -2.898 645 28 -2.898 645 47 -2.898 645 58 -2.898 645 66

6 -2.903 658 08 -2.903 658 27 -2.903 658 39 -2.903 658 46

10 -2.903 716 91 -2.903 717 10 -2.903 717 22 -2.903 717 30

15 -2.903 722 72 -2.903 722 91 -2.903 723 03 -2.903 723 10

21 -2.903 723 54 -2.903 723 74 -2.903 723 85 -2.903 723 93

28 -2.903 723 55 -2.903 723 74 -2.903 723 85 -2.903 723 93

35 -2.903 723 91 -2.903 724 03 -2.903 724 10

40 -2.903 724 03 -2.903 724 10

45 -2.903 724 15

[9] -2.903 722 99

[15] -2.903 724 37

5. Conclusion

The main results presented in the paper were summarized in the abstract. The numerical analysis shows that the finite element discretization of the problem and developed numerical schemes and algorithms implemented in Fortran 77 as the program complex POTHEA provide stable calculations with a specified accuracy of ~ 10-12 of

Table 6

Convergence of the first exited state energy (in a.u.) like Table 5

N Jmax 21 jmax — 28 jmax — 35

1 -2.139 935 59 -2.139 935 68 -2.139 935 71

2 -2.141 664 27 -2.141 664 32 -2.141 664 34

6 -2.145 700 08 -2.145 700 17 -2.145 700 20

10 -2.145 914 95 -2.145 915 04 -2.145 915 07

15 -2.145 957 21 -2.145 957 30 -2.145 957 34

21 -2.145 968 71 -2.145 968 74

28 -2.145 970 24

N Jmax — 40 Jmax — 45 jmax — 50

1 -2.139 935 72 -2.139 935 72 -2.139 935 73

2 -2.141 664 35 -2.141 664 35 -2.141 664 36

6 -2.145 700 21 -2.145 700 21 -2.145 700 22

10 -2.145 915 09 -2.145 915 09 -2.145 915 10

15 -2.145 957 35 -2.145 957 36 -2.145 957 36

21 -2.145 968 76 -2.145 968 76 -2.145 968 77

28 -2.145 970 26 -2.145 970 26 -2.145 970 27

35 -2.145 972 10 -2.145 972 10 -2.145 972 11

40 -2.145 972 62 -2.145 972 63

45 -2.145 973 22

[9] -2.145 956 97

[16] -2.145 974 04

parametric eigenfunctions and their derivatives with respect to the parameter within an accuracy of the same order 0(hp+1), as well as its parameter eigenvalues, their derivatives and the matrix elements within an accuracy of the same order 0(h2p) in step h of the finite element grid, in accordance with theoretical estimations at fixed number jmax of equations (51). However, slow convergence of eigenvalues and matrix elements vs the number jmax of Legendre polynomials of expansion (5) which is proportional to the inverse cube of their order was shown. This fact gives restriction on the required accuracy higher than ~ 10-8 of solving the parametric 2D BVPs and calculations of the matrix elements by using the conventional expansion (5) at reasonable computer resources. The benchmark calculations of the matrix elements and eigenenergies of a Helium atom confirm applicability of the elaborated algorithms and program complexes taking into account the above restrictions. Solutions of this problem can be got by using the multistep Kantorovich expansion [12], for example, Eq. (6) in hyperspheroidal coordinates [13,14], which will be considered in the following papers.

Thus, the calculated parametric eigenvalues, eigenfunctions (parametric basis functions) and the matrix elements can be used for the numerical solution with the required accuracy of bound states and scattering problem for the three dimensional equation of the Schrodinger type, including long-range potentials of the Coulomb type or for various three-dimensional elliptic equations in partial derivatives, with the help of the program complexes POTHEA and KANTBP [28-31]. The generalization of the algorithm for solving a system of the parametric coupled 2D BVPs in the framework of the projection method and FEM, which can be applied for solving multidimensional boundary value problems for equations of Schrodinger type, will be given in further papers.

References

1. Macek J. Properties of Autoionizing States of He // J. Phys. B. — 1968. — Vol. 1. — Pp. 831-843.

2. Hornos J. E, MacDowell S. W, Caldwell C. D. Two-Electron Wave Functions in Hyperspherical Coordinates // Phys. Rev. A. — 1986. — Vol. 33. — Pp. 22122224.

3. Two-Dimensional Basis Functions for the Three-Body Problem in Hyperspherical Coordinates / A. G. Abrashkevich, S. I. Vinitskii, M. S. Kaschiev, I. V. Puzynin // Sov. Nucl. Phys. — 1988. — Vol. 48. — Pp. 602-608.

4. Convergence of the Hyperspherical Adiabatic Expansion for Helium-Like Systems / A. G. Abrashkevich, D. G. Abrashkevich, M. S. Kaschiev et al. // J. Phys. B: At. Mol. Opt. Phys. — 1989. — Vol. 22. — Pp. 3957-3963.

5. Multichannel Calculation of the Electric-Dipole Oscillator Strengths for the Discrete 1S'e - 1P° Transitions in Helium within the Hyperspherical Adiabatic Approach / A. G. Abrashkevich, D. G. Abrashkevich, M. I. Gaysak et al. // Physics Letters A. — 1991. — Vol. 152. — Pp. 467-471.

6. Adiabatic Hyperspherical Representation in Barycentric Coordinates for HeliumLike Systems / A. G. Abrashkevich, D. G. Abrashkevich, I. V. Puzynin, S. I. Vinit-sky // J. Phys. B. — 1991. — Vol. 24. — Pp. 1615-1638.

7. Doubly Excited States of H- and He in the Coupled-Channel Hyperspherical Adiabatic Approach / A. G. Abrashkevich, D. G. Abrashkevich, M. S. Kaschiev et al. // Phys. Rev. A. — 199. — Vol. 45. — Pp. 5274-5277.

8. Abrashkevich A. G., Abrashkevich D. G., Shapiro M. HSTERM - A Program to Calculate Potential Curves and Radial Matrix Elements for Two-Electron Systems within the Hyperspherical Adiabatic Approach // Comput.Phys. Comm. — 1995. — Vol. 90. — Pp. 311-339.

9. De Groote J. J., Masiliy M., Hornos J. E. Highly Excited States for the Helium Atom in the Hyperspherical Adiabatic Approach // J. Phys. B. — 1998. — Vol. 31. — Pp. 4755-4764.

10. Abrashkevich A. G., Kaschiev M. S., Vinitsky S. I. A New Method for Solving an Eigenvalue Problem for a System of Three Coulomb Particles within the Hy-perspherical Adiabatic Representation // J. Comp. Phys. — 2000. — Vol. 163. — Pp. 328-348.

11. Abrashkevich A. G., Puzynin I. V., Vinitsky S. I. ASYMPT: a Program for Calculating Asymptotics of Hyperspherical Potential Curves and Adiabatic Potentials // Comput. Phys. Commun. — 2000. — Vol. 125. — Pp. 259-281.

12. Algorithm for Reduction of Boundary-Value Problems in Multistep Adiabatic Approximation / A. A. Gusev, O. Chuluunbaatar, V. P. Gerdt et al. // arXiv:1005.2089. — 2010.

13. Kaschiev M., Vinitsky S. I. Schrodinger Equation for a Three Particle System in Spheroidal Coordinates // Sov. J. Nucl. Phys. — 1986. — Vol. 44. — Pp. 246-260.

14. Abramov D. I. Hyperspherical Coulomb Spheroidal Basis in the Coulomb Three-Body Problem // Phys. Atom. Nucl. — 2013. — Vol. 76. — Pp. 196-207.

15. Chuluunbaatar O, Puzynin I. V., Vinitsky S. I. Uncoupled Correlated Calculations of Helium Isoelectronic Bound States // J. Phys. B. — 2001. — Vol. 34. — Pp. L425-L432.

16. Drake G. W. F., Van Z.-C. Variational eigenvalues for the S states of helium // Chem. Phys. Lett. — 1994. — Vol. 229. — Pp. 486-490.

17. Calculation of a Hydrogen Atom Photoionization in a Strong Magnetic Field by Using the Angular Oblate Spheroidal Functions / O. Chuluunbaatar, A. A. Gusev, V. L. Derbov et al. // J. Phys. A. — 2007. — Vol. 40. — Pp. 11485-11524.

18. POTHMF: A Program for Computing Potential Curves and Matrix Elements of the Coupled Adiabatic Radial Equations for a Hydrogen-Like Atom in a Homogeneous Magnetic Field / O. Chuluunbaatar, A. A. Gusev, V. P. Gerdt et al. // Comput. Phys. Commun. — 2008. — Vol. 178. — Pp. 301-330.

19. The Application of Adiabatic Method for the Description of Impurity States in

Quantum Nanostructures / A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky et al. // J. Phys. Conf. Ser. — 2010. — Vol. 248. — Pp. 012047-1-8.

20. Adiabatic Description Of Nonspherical Quantum Dot Models / A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky et al. // Phys. Atom. Nucl. — 2012. — Vol. 75. — Pp. 1281-1297.

21. Channeling Problem for Charged Particles Produced by Confining Environment / O. Chuluunbaatar, A. A. Gusev, V. L. Derbov et al. // Phys. Atom. Nucl. — 2009. — Vol. 72. — Pp. 811-821.

22. Symbolic-Numerical Algorithms to Solve the Quantum Tunneling Problem for a Coupled Pair of Ions / A. A. Gusev, S. I. Vinitsky, O. Chuluunbaatar et al. // Lecture Notes in Computer Science. — 2011. — Vol. 6885. — Pp. 175-191.

23. Symbolic-Numerical Algorithm for Generating Cluster Eigenfunctions: Tunneling of Clusters Through Repulsive Barriers / A. Gusev, S. Vinitsky, O. Chuluunbaatar et al. // Lecture Notes in Computer Science. — 2013. — Vol. 8136. — Pp. 427-440.

24. Виницкий С. И., Гусев А. А., Чулуунбаатар О. Решение краевых задач шрё-дингеровского типа методом Канторовича // Вестник СПбГУ: Серия 4 «Физика. Химия». — 2010. — Т. 3. — С. 111-115. [Vinitsky S.I., Gusev A.A., Chuluunbaatar O. Solution of the Schrödinger-Type Boundary-Value Problems by Kantorovich Method // Vestnik SPbGU: Series 4. Physics, Chemistry. — 2010. — Vol. 3. — P. 111-115. ]

25. Kantorovich L. V., Krylov V. I. Approximate Methods of Higher Analysis. — New York: Wiley, 1964.

26. Чулуунбаатар О. Алгоритм численного решения параметрической задачи Штурма-Лиувилля и вычисления производных от решения по параметру методом конечных элементов // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2009. — Т. 2. — С. 54-65. [Chuluunbaatar O. Algorithm of Numerical Solving the Parametric Sturm-Liouville Problem and Calculation of Solution Derivatives with Respect to the Parameter via the Finite-Element Method // Bulletin of Peoples' Friendship University of Russia. Series "Mathematics. Information Sciences. Physics". — 2009. — No 2. — P. 54-65. ]

27. A Symbolic-Numerical Algorithm for the Computation of Matrix Elements in the Parametric Eigenvalue Problem / S. I. Vinitsky, V. P. Gerdt, A. A. Gusev et al. // Programming and Computer Software. — 2007. — Vol. 33. — Pp. 105-116.

28. 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 // Comput. Phys. Commun. — 2009. — Vol. 180. — Pp. 13581375.

29. KANTBP: A Program for Computing Energy Levels, Reaction Matrix and Radial Wave Functions in the Coupled-Channel Hyperspherical Adiabatic Approach / O. Chuluunbaatar, A. A. Gusev, A. G. Abrashkevich et al. // Comput. Phys. Commun. — 2007. — Vol. 177. — Pp. 649-675.

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

31. Chuluunbaatar O, Gusev A. A., Vinitsky S. I., Abrashkevich A. G. A Program Package for Solution of Two-Dimensional Discrete and Continuum Spectra Boundary-Value Problems in Kantorovich (Adiabatic) Approach. — 2013. — http://wwwinfo.jinr.ru/programs/jinrlib/kantbp/indexe.html.

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

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

34. Chuluunbaatar O. The Scientific Doctoral Thesis. — JINR, 2010.

35. Gantmacher F. R. The Theory of Matrices. — USA: AMS, Providence, 2000.

36. Golub G. H., Van Loan C. F. Matrix Computations. — Johns Hopkins Univ. Press, 1996.

37. Whitakker E. T., Watson G. N. A Course of Modern Analysis. — Cambridge: Univ. Press, 1927.

УДК 517.958:530.145.6

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

элементов

Л. Л. Гусев

Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, г. Дубна, Московской обл., 141980, Россия

Представлены эффективные и стабильные алгоритмы численного решения с заданной точностью параметрической двумерной краевой задачи на собственные значения (КЗСЗ). КЗСЗ формулируется для самосопряженного эллиптического дифференциального уравнения в частных производных с краевыми условиями Неймана и/или Дирихле в конечной двумерной области. Исходная задача редуцируется к параметрической однородной одномерной КЗСЗ для системы обыкновенных дифференциальных уравнений второго порядка (ОДУ). Редукция производится разложением искомого решения по подходящему набору ортогональных собственных функций вспомогательной задачи Штурма—Лиувилля по одной из переменных. Производные по параметру от собственных значений и соответствующих собственных вектор-функций редуцированной задачи определяются как решения параметрической неоднородной одномерной КЗСЗ, полученной дифференцированием по параметру редуцированной задачи. Полученные КЗСЗ решаются методом конечных элементов с автоматическим выбором сдвига спектра. Алгоритм, реализованный на Фортране 77 в виде программы РОТНЕА, вычисляет с заданной точностью набор ~ 50 собственных значений (потенциальных термов), собственных функций и их первых производных по параметру, а также матричных элементов — интегралов от произведения собственных функций и/или первых производных собственных функций по параметру. Вычисленные потенциальные термы и матричные элементы можно использовать для формирования матрицы переменных коэффициентов системы ОДУ, которая возникает при редукции трёхмерной КЗСЗ в рамках многоканального адиабатического подхода или метода Канторовича. Эффективность и стабильность алгоритма продемонстрирована численным анализом собственных решений параметрической двумерной КЗСЗ и вычисленных матричных элементов которые применяются при решении с помощью программы КАМТБР трёхмерной КЗСЗ для уравнения Шрё-дингера для атома гелия с нулевым полным угловым моментом в гиперсферических координатах.

Ключевые слова: параметрическая двумерная задача на собственные значения, эллиптическое уравнение второго порядка, метод конечных элементов, метод Канторовича, гиперсферические координаты, атом гелия.