CC BY
30
UDC 517.958:530.145.6
Description of a Program for Computing Eigenvalues and Eigenfunctions and Their First Derivatives with Respect to the Parameter of the Coupled Parametric Self-Adjoined Elliptic Differential Equations
A. A. Gusev*, O. Chuluunbaatar^, S. I. Vinitsky*, A. G. Abrashkevich*
* Joint Institute for Nuclear Research 6, Joliot-Curie, Dubna, Moscow region, Russia, 141980 ^ School of Mathematics and Computer Science National University of Mongolia, Mongolia 1 IBM Toronto Lab, 8200 Warden Avenue, Markham, ON L6G 1C7, Canada
Brief description of a FORTRAN 77 program is presented for calculating with the given accuracy eigenvalues, eigenfunctions and their first derivatives with respect to the parameter of the coupled parametric self-adjoined elliptic differential equations with the Dirichlet and/or Neumann type boundary conditions on the finite interval. The original problem is projected to the parametric homogeneous and nonhomogeneous 1D boundary-value problems for a set of ordinary second order differential equations which is solved by the finite element method. The program calculates also potential matrix elements - integrals of the eigenfunctions multiplied by their first derivatives with respect to the parameter. Parametric eigenvalues (so-called potential curves) and matrix elements computed by the POTHEA program can be used for solving the bound state and multi-channel scattering problems for a system of the coupled second-order ordinary differential equations with the help of the KANTBP programs. As a test desk, the program is applied to the calculation of the potential curves and matrix elements of Schrodinger equation for a system of three charged particles with zero total angular momentum.
Key words and phrases: boundary value problem, finite element method, Kantorovich method.
1. Introduction
In this work we present a brief description of a POTHEA program for calculating with a given accuracy eigenvalues, eigenfunctions and their first derivatives with respect to the parameter of the coupled parametric self-adjoined elliptic differential equations with the Dirichlet and/or Neumann type boundary conditions on the finite interval [1]. The original problem is projected to the parametric homogeneous and nonhomogeneous 1D BVPs for a set of ordinary second order differential equations which is solved by the finite element method [2]. The program calculates also potential matrix elements - integrals of the eigenfunctions multiplied by their derivatives with respect to the parameter.
Potential curves and matrix elements computed by the POTHEA program can be used for solving the bound state and multi-channel scattering problems for a system of the coupled second-order ordinary differential equations with the help of the KANTBP programs [1,3].
As a benchmark, we present calculation with a given accuracy of potential curves and matrix elements which is applied for calculation of ground state energy and first exited state energy of an helium atom in the framework of the Kantorovich method implemented like the close-coupled hyperspherical adiabatic approach [4]. The numeric results show that the program developed is very efficient and allows to obtain numerical solutions of the above problems with the required accuracy using very little computational resources.
Received 27th September, 2013. The work was supported partially by grants 13-602-02 JINR, 11-01-00523 and 13-01-00668 RFBR..
2. Statement of the Problem
Let us consider a boundary problem for a parametric two dimensional self-adjoined second order ordinary differential equation on the region Q = [xmin, xmax] x [ymin, ymax]
1 d d 1 1 d d \
h(y)^d h(x) — + U (x, y;z) —e(z))B (x, y;z) = 0. (1)
/i( y)dy dy f?,(y) f4(x)dx dx
with the Dirichlet and/or Neumann type boundary conditions (t = min, max)
lim f2(y)dyB(x, y;z) = 0 or B(x, yt;z) = 0, x G (xm\n,xmax), (2)
y^yt
lim f5(x)dxB(x, y;z) = 0 or B(xt, y;z) = 0, y G (ymin, ymax).
Here z G [zmin, Zmax] is a parameter, functions fi(y) > 0, f2(y) > 0, f3(y) > 0, f 4(x) > 0, f 5(x) > 0, and dy f2(y), dxfs(x), U(x, y;z), dzU(x, y;z) are continuous on the (x, y) G Q/dQ. Also assume that the parametric boundary value problem (BVP) (1), (2) has only discrete spectrum.
The program executes the following steps.
In Step 1 program calculates a set of jmax smallest eigenvalues s1 ( z) < e2(z) < ... < en( z), and £1(z) > a(z), and the corresponding eigenfunctions {Bj(x, y\z)}^^=1 G Fz ~ ~L2(Qx>y), satisfying the orthogonality and normalization conditions
^max
J dy fi(y) j dx f4(x)Bi(x, y-,z)Bj(yar^ y;z) = Sij, (3)
^min ^min
where Sij is the Kronecker symbol, and a(z) > —to is the lower bound of the smallest eigenvalue of e1(z).
In Step 2 program 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 program computes matrix elements defined by the integrals
ymax Xmax
Hij (z) = Hji(z) = J dy ¡1(y) J dx f4(x)дzBi(x, y;z)dzBj(x, y■,z), (4)
ymin ^min
ymax ^max
Qij (z) = —Qji(z) = — J dy h (y) J dx f4(x)Bi(x, ¥■ z)dzBj(x, y;z).
ymin ^min
with an accuracy of the same order achieved for the corresponding eigenvalues of the BVP (1)-(3).
2.1. Reduction of the 2D BVP to the 1D BVP
Step 1.1. The partial wave function Bi(x, y; z) is expanded over the orthonormal basis functions {^j (x)}£=T:
Jmax
Bi(x, y;z) = ^, (x) $\y;z). (5)
3 = 1
In Eq. (5), the vector-functions £(i}(y; z) = (^ ( y; z),..., ( y; Z))T are unknown. The functions (x) are determined as solutions of the following eigenvalue problem (t = min, max):
(-¿)dx /s(x)dx + Uo(x)) (x) = ^(x), (6)
/5(x)ddxxl = 0 or 1>j(xt) = 0,
where U0(x) is a known function and
dx /4 (x)%pi(x)%pj (x) = 5ij. (7)
Note, this problem can be numerically solved with a given accuracy by mens of the ODPEVP program [23].
Step 1.2. After minimizing the Rayleigh-Ritz variational functional, and using the expansion (5), Eq. (1) is reduced to a finite set of jmax ordinary second-order differential equations (t = min, max)
1 B B
(D(y; z) - Ei(z) I) £(i}(y; z) = 0, D( y; z) = - j-^I-f2(y)— + W(y; z), (8)
lim h(y)ByH(i}(y;^) = 0 or ; z) = 0. (9)
y^yt
Here I, W(y; z) are symmetric matrices of dimension jmax x jmax
ymax
hj = ôij = i dy h (y) {n(i)(y;z))T Ç(j}(y;z), (10)
7ij(y;z) = yz+^f5*j + j dxÎ4(x)iKx)[UU(x,y;z)- (
y
Step 2. Taking a derivative of the boundary problem (8)-(9) with respect to parameter z, we get that Bz£(t\y,z) can be obtained as a solution of the following boundary problem (t = min, max)
(D(y;z) - £l ( z) I)
BÇ(Z} (y;z) B
— (W(y;z) - £i(z) I)
^ }(y;z),
lim f2(y)By(BzH(%}(y;z)) = 0 or Bz^(€}(yt; z) = 0.
y^yt
The BVP (11)-(12) has a unique solution, if and only if:
(11) (12)
Be j(z) B
f dy h {y) ^}(y;z))T B-zWf^^}(y;z), (13)
y
j dvMv) (i%;r))T d^dp)=0. (14)
^min
Step 3. The required matrix elements (4) are represented by
Umax / / .x \ T (
di( %\y;,z)\ d^(y;;z)
¡/max / ( \ .
Hl3 (z) = Hji(z)= J dy h(y)l 9^(f^J
dz
(15)
Td Î (j)(y;z)
Qij (z) = -Qn(z) = -J dy h(y) (t(i)(y;z)) dz
3. Test Desk
For a Helium atom with zero angular momentum in hyperspherical coordinates x = a, y = z = R, by using weight functions f1(x) = f2(x) = f3(x) = sin2 a, fi(y) = h( ?/) = sintf and potential function
U(a,0;R) = R - . --+ 1 ) ,
2 y sin(a/2) cos(a/2) - sin( a) costf J
we reduce boundary value problem (1)-(3) using expansion (5) with analytical solution of problem (6)-(7) in the form of Legendre polynomials Pj (cos d) to boundary problems (8)-(14). The later have been solved by the POTHEA program on finite element grids = {0.(150)^/2} with fourth-order Lagrange elements with accuracy eps = 10-12 and the corresponding matrix elements H( z), Q(-z) have been calculated with accuracy eps = 10-6, at run time is 4 seconds.
The following values of numerical parameters and characters described in [1] have been used in the test run via the supplied input file POTHEA.INP:
&PARAMS TITLE=' PARAMETRIC 2D DIFFERENTIAL EQUATION
ICOUN=0,PARAM=10D0,NROOT=6,MDIM=12,NPOL=4,RTOL=1.D-12, NITEM=2000,SHIFT=-1.1D0,ICHK=1,IPRINT=0,IPRSTP=15, NMESH=3,RMESH=0.0D0,150.D0,1.570796326794896D0, NDIR=1, NDIL=12, NMDIL=0,IBOUND=4,
FNOUT='3DNGSS.LPR,,IOUT=7,POTEN=,3DNGSS.PTN,,IOUP=10, FMATR='3DNGSS.MAT',IOUM=11,EVWFN=,3DNGSS.WFN,,IOUF=1
&END
All calculation details of this problem were written into file POTHEA.LPR.
Test Run Output
PROBLEM: PARAMETRIC 2D DIFFERENTIAL EQUATION
CONTROL INFORMATION
y
NUMBER OF DIFFERENTIAL EQUATIONS NUMBER OF ENERGY LEVELS REQUIRED NUMBER OF FINITE ELEMENTS . . .
(MDIM ) = (NROOT ) = (NELEM ) =
12 6
150
NUMBER OF GRID POINTS..........(NGRID ) = 601
ORDER OF SHAPE FUNCTIONS........(NPOL ) = 4
ORDER OF GAUSS-LEGENDRE QUADRATURE . . . (NGQ ) = 5
NUMBER OF SUBSPACE ITERATION VECTORS. . . (NC ) = 12
BOUNDARY CONDITION CODE ......... (IBOUND) = 4
SHIFT OF EIGENVALUE...........(SHIFT ) = -1.10000
CONVERGENCE TOLERANCE .......... (RTOL ) = 0.100000E-11
VALUE OF PARAMETER............(PARAM ) = 10.0000
SUBDIVISION OF RHO-REGION ON THE FINITE-ELEMENT GROUPS:
NO OF NUMBER OF BEGIN OF LENGTH OF GRID END OF GROUP ELEMENTS INTERVAL ELEMENT STEP INTERVAL
1 150 0.000 0.01047 0.00262 1.571
T O T A L S Y S T E M D A T A
TOTAL NUMBER OF ALGEBRAIC EQUATIONS. . . . (NN ) = 7212
TOTAL NUMBER OF MATRIX ELEMENTS......(NWK) = 262878
MAXIMUM HALF BANDWIDTH..........(MK ) = 60
MEAN HALF BANDWIDTH..........(MMK) = 36
NDIM, MDIM= 12 12
THERE ARE 0 ROOTS LOWER THEN SHIFT CONVERGENCE REACHED FOR RTOL 0.1000E-11 ITERATION NUMBER 68 RELATIVE TOLERANCE REACHED ON EIGENVALUES 0.0000E+00 0.4385E-15 0.5146E-14 0.1811E-13 0.7383E-15 0.1124E-11
R O O T N U M B E R E I G E N V A L U E D E R I V A T I V E
1 -106.1449119429294 -20.49786509377458
2 -32.40954538140649 -5.355570130015161
3 -30.37792481031590 -5.456241767021648
4 -21.97501254692271 -3.414051588534342
5 -19.24789115244545 -3.025292223400815
6 -15.24989121658901 -2.766774467691610
POTENTIAL MATRICES H(I,J) AND Q(I,J):
H-MATRIX 0.7530D-02 0.7277D-02 -.4429D-02 -.2551D-02 0.1761D-02 -.6160D-03
AT THE PARAMETER = 0.7277D-02 -.4429D 0.2831D-01 -.4826D-02 0.8100D-02 -.7355D-03 0.9388D-03
02
-.4826D-02 0.1766D-01 -.1222D-02 0.5248D-02 -.2412D-02
10.00000 -.2551D-02 0.8100D-02 -.1222D-02 0.2799D-01 -.4505D-02 0.1145D-02
0.1761D-02 -.7355D-03 0.5248D-02 -.4505D-02 0.1618D-01 -.7024D-02
-.6160D-03 0.9388D-03 -.2412D-02 0.1145D-02 -.7024D-02 0.1048D-01
Q-MATRIX AT THE PARAMETER
0.1003D-14 -.4752D-01 0.2558D-01 -.2470D-01 0.1301D-01 -.6895D-02
0.4752D-01 0.2947D-15 -.1519D-01 -.1397D+00 0.2136D-01 -.5194D-02
-.2558D-01 0.1519D-01 -.2203D-15 0.1770D-02 -.9586D-01 0.5406D-01
10.00000 0.2470D-01 0.1397D+00 -.1770D-02 0.5867D-16 0.2047D-01 -.3502D-02
-.1301D-01 -.2136D-01 0.9586D-01 -.2047D-01 -.9482D-16 0.8956D-02
0.6895D-02 0.5194D-02 -.5406D-01 0.3502D-02 -.8956D-02 0.3119D-16
References
1. A Program Package for Solution of Two-Dimensional Discrete and Continuum Spectra Boundary-Value Problems in Kantorovich (Adiabatic) Approach / O. Chuluun-baatar, A. A. Gusev, S. I. Vinitsky, A. G. Abrashkevich // JINR Lib. — 2013. — http://wwwinfo.jinr.ru/programs/jinrlib/kantbp/indexe.html.
2. Gusev A. A. 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 // Bulletin of PFUR. "Series Mathematics. Information Sciences. Physics". — 2013. — No 1. — Pp. 10-31.
3. 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. Abrashke-vich // Comput. Phys. Commun. — 2008. — Vol. 179. — Pp. 685-693.
4. 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 Hyperspherical Adiabatic Representation // J. Comp. Physics. — 2000. — Vol. 163, No 2. — Pp. 328-348.
УДК 517.958:530.145.6
Описание программы вычисления собственных значений и
собственных функций и их первых производных по параметру для параметрической самосопряжённой системы эллиптических дифференциальных уравнений
А. А. Гусев*, О. Чулуунбаатар+, С. И. Виницкий*, А. Г. Абрашкевич*
* Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, г.Дубна, Московская область, Россия, 141980 ^ Факультет математики и компьютерных наук Монгольский государственный университет, Монголия * Лаборатории IBM в Торонто 8200 Ворден авеню, Маркхэм, ON L6G 1C7, Канада
Представлено краткое описание программ на языке Фортран 77 для расчёта с заданной точностью собственных значений, собственных функций и их первых производных по параметру для параметрической самосопряжённой системы эллиптических дифференциальных уравнений на конечном интервале с граничными условиями Дирихле и/или Неймана. Исходная задача проецируется на параметрические однородные и неоднородные одномерные краевые задачи для системы обыкновенных дифференциальных уравнений второго порядка, решаемые методом конечных элементов. Программа рассчитывает также потенциальные матричные элементы — интегралы от собственных функций, умноженные на их первые производные по параметру. Собственные значения, зависящие от параметра (так называемые потенциальные кривые) и матричных элементов, рассчитываемые программой POTHEA, могут быть использованы для решения с помощью программы KANTBP задач на связанные состояния и многоканальные задачи рассеяния для системы второго порядка обыкновенных дифференциальных уравнений. В качестве теста программа использована для расчёта потенциальных кривых и матричных элементов уравнения Шрёдингера для системы трёх заряженных частиц с нулевым полным угловым импульсом.
Ключевые слова: краевая задача, метод конечных элементов, метод Канторовича.
