Algorithms and programs for solving boundary-value problems for systems of second-order ODEs with piecewise constant potentials: multichannel scattering and eigenvalue problems Текст научной статьи по специальности «Математика»

UDC 519.632.4

Algorithms and Programs for Solving Boundary-Value Problems for Systems of Second-Order ODEs with Piecewise Constant Potentials: Multichannel Scattering and Eigenvalue Problems

A. A. Gusev*, O. Chuluunbaatar*, S. I. Vinitsky*t, L. L. Hai**,

V. L. DerbovS, A. Gozdz32

* Joint Institute for Nuclear Research, Dubna, Russia ^ RUDN University, Moscow, Russia * Belgorod State University, Belgorod, Russia S Saratov State University, Saratov, Russia 32 Institute of Physics, University of M. Curie-Sklodowska, Lublin, Poland

The new algorithms and programs, implemented in Maple, for solving waveguide-type multichannel scattering and eigenvalue problems for systems of the second-order ODEs with N x N matrix piecewise constant coefficients on the axis are proposed. New algorithm and program for solving the boundary-value problems by method of matching the fundamental solutions (MMFS) of the system of ODEs at the points of discontinuity of potentials are elaborated. In each of the subintervals of an axis the general solution of the system of ODEs are sought in the form of linear combination of 2N fundamental solutions with unknown coefficients. Each fundamental solution explicitly dependent on spectral parameter and eigenvalues and eigenvectors of algebraic eigenvalue problems with N x N matrix of constant potentials. From the condition of continuity for the solutions and their derivatives at the discontinuity points of the potentials, the system of algebraic equations is followed. In the case of bound or metastable state problem the obtained system of algebraic equations contains nonlinear dependence of unknown spectral parameter. For solving such nonlinear problem symbolic-numerical algorithm is formulated. The benchmark calculations of bound, metastable and scattering states of BVPs for systems of the second-order ODEs obtained using program of the MMFS are compared with those obtained using program of the finite element method.

Key words and phrases: multichannel scattering problem, eigenvalue problem, system of second order ordinary differential equations, method of matching the fundamental solutions

1. Introduction

The boundary-value problems for systems of N second-order ordinary differential equations (ODEs) of the waveguide type with the matrix of piecewise constant potentials arise in mathematical modelling of quantum-dimensional nanostructures and optical multilayer systems [1]. For example, application Kantorovich method to the solution of Maxwell's equations in an integrated optical waveguide with an irregular change of parameters along two horizontal directions is considered in Ref. [2]. In this paper the BVP for a system of differential equations for the coefficient Kantorovich functions using computer algebra system Maple was derived. For solving such type of the BVPs the method of matching the fundamental solutions (MMFS) at each boundary between the adjacent axis subintervals was applied for single second-order ODE, while application to system of ODEs is required [3].

The generalization of MMFS and creation of effective algorithms and programs for solving the BVPs for N second order ODEs with piecewise constant potentials in axis is the aim of the present paper. The efficiency of the algorithm and program is demonstrated by the example of calculating the resonance and metastable states of the multichannel scattering problem in the axis for the set of N second order ODEs with piecewise constant

Received 12th July, 2016.

The work was supported by the Russian Foundation for Basic Research, grant No. 14-01-00420 and Bogoliubov-Infeld JINR program. The reported study was partially funded within the Agreement No. 02.a03.21.0008 dated 24.11.2016 between the Ministry of Education and Science of the Russian Federation and RUDN University.

potentials. The results of solving the boundary-value problem by program of the MMFS in axis are compared with those obtained using program KANTBP 4M [4] of the finite element method (FEM) in finite interval with appropriate boundary conditions.

The structure of the paper is following. In Section 2 we present the algorithm of MMFS for solving the BVPs for systems of N second-order ordinary differential equations with piecewise constant potentials and calculation of bound, metastable and scattering states. In Section 3 we compare the benchmark calculations MMFS of bound, metastable and scattering states with results obtained with the program KANTBP 4M that realizes the Finite Element Method (FEM).

2. Algorithm of matching the fundamental solutions

Let us consider the boundary-value problems for systems of N second-order ordinary differential equations (ODEs) of the waveguide type

(-I ^ + V(z) - E l)$(z) = 0, *(z) = ($1(z),..., $n (z))T, (1)

with the matrix of piecewise constant potentials V(z) = {Vji(z)}, i,j = 1, ...,N:

VJi(z) = {Vji;1,Z < Zl,Vji;2,Z < Z2,...,Vji;k-1,Z < Zk-1,Vji-,k,z>zk-i}. (2)

Here we suppose that the matrix of potentials is symmetric Vji(z) = Vij(z) and real-valued Im Vji(z) = 0, the case of complex valued potentials is considered in [5].

Step 1. In each of the subintervals z G (zm-1,zm), m = 1,...,k, the system of N ODEs (1) is a system of ODEs with constant coefficients

(-I^ + Vm - E i) ^m(z) = 0,

that has the general solution, explicitly depending upon the spectral parameter E

N

*m(z)= £ (c{2m-m+tFt](E, -z)^^m)+C{2m-1)N+tF(m)(E,z)^^m)) . (3)

i=1

Here Fj;m)(E,z) = exp ^X^-E^j, C2mN+i are unknown coefficients, X^ and = {^l^,..., ^Nl)}T are the solutions of the algebraic eigenvalue problem

Vm^[m) = x(T)^<T), (4)

calculated numerically. For the symmetric real matrix Vm the eigenvectors obey

the condition of orthogonality and normalisation (^¿m))T= ôij.

Step 2.1. For bound states ,E<^max= min^l^,..., A^, Alk),..., A^) the asymptotic conditions at z ^ describe the exponentially decreasing solutions &m(z) of

(3) N

(z^-œ)= £ Cn+ expf/\(1)-W) *(1), (5)

i=i ^ '

$fc(z^+œ)= £ C(2k-2)N + exp f-^A(k)-W) *(k). (6)

¿=1 ^ '

The wave functions of bound states satisfy to the normalization conditions

+œ

J dz&(z)$(z) = 1, (7)

—œ

where t means a Hermitian conjugation.

Step 2.2. For metastable states EM =Re EM +i Im EM, Im EM<0 the asymptotic conditions at z^-œ describe outgoing waves in the open channels and exponentially decreasing solutions in the closed ones

n ( Cn + expf A(1)-£MÀ Re EM < X^,

*i(z ^-œ) = £*(1) - { (1) (8)

¿=i I C,expi -%\JEM-X^z), ReEM > A,(1),

N f C(2k—2)N+iexJ-fifT^z), ReEM < A(k),

$k(z ^ +œ) = V > ,-{ (9)

-1 I C(2k—1)N+iexJ i \JEM-\(k)zj, ReEM > X^.

The wave functions of metastable states satisfy to the normalization conditions

+ œ

J dz®T(z)&(z) = 1 (10)

—œ

where T means a transposition.

Step 2.3. For multichannel scattering problem with fixed energy E = Re E the desired matrix of solutions is calculated with the asymptotic form "incident wave + outgoing waves" (see Fig. 1)

, , / xfcax (Z)t^+xmL(z)T% *> o,

^ (z ^ ±œ) = { , , < ^ (11)

i xmnw+xan)(z)^+X^I^R:, * < o, ^ , f xmx(z)+xmtx(^)r^+xmu)Rt, *> o,

(z ^ ±œ) = { ^ ) ( : (12)

i xmn(,)T :+xmi^T t, o.

Here R: and R: are desired square matrices of reflection amplitudes with the dimension Nq x Nq and N^ x N^ and T: and T : are the rectangular matrices of transmission amplitudes with the dimension N^ x N% and xN^, where and N^ are the numbers of open channels in the asymptotic regions. Here the desired scattering matrix S with the dimension N0 x N0 is calculated

S=(R: R:), sts=sst=i, (13)

which is unitary and symmetric for real effective potentials.

Figure 1. Left and right asymptotic solutions "incident wave + outgoing waves",

respectively

Let us write the expressions (11)-(12) in the explicit form

( (

exp l

N

(1)

exp i » JEM-X^z

hem _x(1)

-Su L +

=1

+ Rii L -

ex^-% JE^z)

hem-xp

N

R^l exp(\f\(1)-EMzj,

E € X

(i)

(k)

Tii t/EM-xW

E>\

( k)

=1

T'l. expi-^A<k)-EMA E € A((k), )

N

(1)

T-

ex^ -% JEM-Xli1)z

^(z ^ +<») = X) *

N

( k)

T^ exp(yj\^-EMz^, E €

' expf— ^EM-X[k)z"]

-Su R +

)

=1

4 EM-x\k) +Ri iR

ex^ » VEM-X\k)z

RriR exp ( -A/\^-EM z),

4/em-x[k)

, E>^i k), E € A(k).

(14)

(15)

(16)

(17)

Step 3. From the condition of continuity for the solutions and their derivatives lim $m-1(z) - $m(z) = 0, lim - m = 2,...,k, (18)

dz

dz

at the discontinuity points z = zm-\, m = 2,...,k of the potentials (2), the system of algebraic equations is followed.

m

m

Step 3.1 In the case of a multichannel scattering problem for the wave incident from the left (see Fig. 1), for each value i% = 1,...,N% we have a system of 2N(k-1) inhomogeneous linear equations with 2 N( k-1) unknowns, C2N+i;s,...,C2n(k-1);s, and RiiL, i = 1..Ng, R^, i = 1 + N^,..., N, TUL, i = 1N, T^, i=\ + Ng ,..., N:' from condition of continuity for the solutions:

RilLF^)(E, Z1 ), E ^X^ f) *

^((x^^,* z1) +|

" )=o-

-C2N^2)(E, -Z1)^2)-C3N^2)(E, Z1)^2) 1=0, (19)

N

£(C(2m—4)N +Ft—1)(E, -Zm—1)^m—1)+C{2m—3)N+iF^m—1)(E, Zm—1)^m—1)

=1

-C{2m—2)N +iFjT)(E, -Zm—1)^m)-C(2m—1)N+iF$m) (E, Zm—1)*[=0, (20)

N

£ 2k —4)N+*FÎk — 1)(E, -Zk—1)^(ik—1)+C{2k —3)N^F}k—1) (E, Zk—1)*(k—1)-i=1\

Tu èx(k) (E,i Zk—1), E>X(k T^F(k)(E, - zk—1), E ^X(k)

from condition of continuity for the derivatives of solutions:

*( k)) =0;

(21)

Nff ( Ru ^(E, -Z1), E>x(k 1 \_(1)_

^^(E, ,1), E <A(k) N 4

(E,z Z1) Su l +|

- C2n+iG(^)(E, -Z1 )^(t)-C3N+iG^i)(E, ¿n)*?^ =0,

(22)

N

C(2m—4)N +^^(E, -Zm—1)*(r—1)+C(2m—3)N+iG^1)(E, Zm—1)*^^

=1

G(2m—2)N+lG<T)(E, - Zm—1)*(T)-C(2m—1)N +iGf) (E, =0, (23)

N (

M

=1

C(2k—4)N+iG(k—1 (E, -Zk—1)*(k—1)+C(2k—3)N+iG(k—1 (E, Zk—1)*(k—1)-

f TuLXC(ik)(E,zZk—1), E>X(k)

T^iLG(k)(E, - zk—1), E ^^ '=0, (24)

}*( =0,

where the following notations are used:

G(m)(E, ±zm) = \(T) - E F(m)(E, ±zm), m = 3,...,k - 1, (25)

x(m)(E, ±izm) = F(m)(E, ±izm)/4E - \t\ m = 1,k, Wm 4-,Zm) = F(m)(E, ±izm)t/n-

X(m)(E, ±izm) = F(m)(E, ±izm) 4 E - \\ ', m = 1,k.

In the case of multichannel scattering problem for the wave incident from the right (see Fig. 1) for each value = 1,...,N^ we have the system of 2N(k-1) inhomogeneous linear equations with 2N(k-1) unknowns, C2n+1;s,..., G2n(fc-i);s, and R^ r , i — 1..N^, RIk, i = 1 + N...,N, TUR, i = 1..N^, T<rtR, i = 1 + N..., N: °

from condition of continuity for the solutions:

(i) (i)

^ (J R (E, ~*Z1), E>\(1' 1(1) M+l Fz(1)(E,Z1), E € \(1) J * -

-C2N+f(2) (E, -Z1)*(2)-C3N+F(2) (E, Z1)*i2^ =0,

(26)

N

E (c(2—4)N + F(m~ 1)(E, -zm- 1)*(m" 1)+C(2m_ 3)N + F(m~1) (E,Zm. 1)*f"1)

=1

-C(2m-2)N+F(m)(E, -zm. 1)*(m)-C(2m_ 1)n+F(m)(E, Zm. 1)*^) =0, (27)

N (

M

=1

C(2k-4)N+F(k-1)(E, -Zk-1)*t1]+C(2k-3)N+F(k-1) (E,Zk-1)*(k-1)-Rii R x( k)(E,iZk-1), E>\( k)

f Rur^(E^z^) E>\\ ' 1 (k) _ (k) \

from condition of continuity for the derivatives of solutions:

(1) (1)

V(J R^ (^ E>X^) 1*(1)_

kv\ TURG(1)(E,Z1), E € A,(1) / * -

R 0-

-C2N+Gf\E, -Z1 )*t(2)-C3N+Gf)(E,Z1)*f) 1=0, (29)

N

E (c(2—4)N + Gim-1)(E, -zm-1)*m-1)+C{2TO_3)N+G\m-1)(E,zm-1)*\m-1)

=1

-C(2m-2)N+GiT)(E, -zm-1)*\m)-c{2m_1)N+G^ (E, zm-1)*\m)) =0, (30)

N

£(C{2k-4)N+iG(k-1)(E, - zk-l]^f-l)+C(2k-3)N +iG(k-1)(E, zk-l)^f-l)-

i=1

(

I R^gXf(E,lZk-t), E>A</k) \ ,k) (k),„ ,s ^ „

where the notations from (25) are used. Step 3.2 or Step 3.3 for the problem of bound or metastable states the nonlinear system of 2N( k-1) + 1 equations with 2N( k-1) + 1 unknowns follows

N

£(C(2m-4)N+Ft-1)(E, -zm-l)^t-l)+C{2m-,)N+Ft-1)(E, zm-l)^t-l)-

i=1

-C(2m-2)N+lF^m)(E, -zm-i)^ln)-C(2m-1)N^F^(E, zm-i)^<fn)) =0, (32)

N

£(C(2m-4)N+G(rl\E, -zm-1)^(m-1)+C(2m-3)N^(^(E,zm-1)*f-1)-

¿=i

-C(2m-2)N+iG(m)(E, -Zm-i)*(m)-C(2m-1)N(E, =0, (33)

where G((m)(E, ±zm) = ±\J\(-m) -EF(m)(E, ±zm), m = 2,..., k. The equation following from the normalisation condition (7) or (10) will be taken into account in final steps.

Step 3.2 From the condition of exponential decay of the bound-state eigenfunctions (5), for z ^ ±ro the limitation for the eigenvalues arises

E<Emax= min(A11) ,...,A(N), A{k) ,...,A(N)) (34)

and the system of equations (32) is completed with the additional conditions

C1 =0,...,Cn = 0, C( 2k-1)N+1 =0,...,C2kN = 0. (35)

The limitation of the parameter E from below

E>Emin= min(A12),..., A{N\..., A1k-1),..., A{k-1)) (36)

follows from the condition of boundedness of the self-adjoint operator, corresponding to the problem (1)—(2) in the inner subintervals z G (zm-1, zm), where m takes the values m = 2, ... , k - 1.

Step 3.3 For metastable states the system of equations (32) instead of (35) is completed with the additional conditions

(C\ = 0, Re EM <A(i1) 1 iCn(2k-1)* = 0, Re EM < ^k 1 , = 12 ^ 1 Cn+г = 0, Re EM > A(1) J , \ Cn {2k-2)i = 0, Re EM > A(k)\ ' ' ,

that follow from the asymptotic conditions (8).

Step 4. Writing the system of equations (32) and (33) in the matrix form

M( E) = 0

and calculating the determinant detM(^) of the obtained matrix, we arrive at the secular equation for the spectral parameter E,

detM(£) = 0. (37)

Remark 1. The determinant of the matrix detM(E^) in the analytical form as a function of the unknown spectral parameter E can be calculated in Maple during the reasonable time only for small number of subintervals k and small number of equations N, kN < 19. For example, for N = 3, k = 4 one expansion of the determinant of the matrix with the dimension 24 x 24 takes more than 3600 seconds, while the construction of the matrix M(E^) and the calculation of its determinant for one given value of E takes about 0.6 second. Therefore, we calculate the determinant detM(^) numerically using the appropriate grid E: Emin,..., Emax from the interval Emin < E < Emax, where the boundary values Emin and Emax are given by Eqs. (34) and (36). One of the possible ways of constructing the suitable grid is implemented using the following algorithm. In the exponential functions are replaced with their truncated expansions in the vicinity of the chosen value E0, i.e., with the polynomials of the order Q with numerical coefficients

, - . , - . exp( ±\J\(m)-Eoz)

exp WA|m)-Ez) = exp iv/A^-^ ^-V , J AE + ..., (38)

V ^ V J 2^^-Eo

where AE = E-E0. Substituting the obtained expressions into the determinant detM(^), we get a polynomial of the order 2N(k - 1)Q of the variable AE, which is again truncated to the order Q, i.e.,

detM(£) = detM(^o) + ai (AE) + ... + aQ(AE)Q. (39)

This leads to the limitation of the grid step A£"max < |edetM(^0)/aQ|1/Q , where £ > 0 is the preassigned number, the accuracy of the determinant detM(^) series expansion. Similar expansion is used to correct the roots of Eq. (37), calculated using a certain method which will be published in a more detail elsewhere.

Remark 2. Near the singular points E = A(m), i = 1,...,^, m = 2,...,k - 1, for which the subradical expressions in the exponents in Eq. (38) turn into zero, the domain of convergence of the expansions (38) and (39) is not large. In this case instead of Eq. (38) the following expansion is used

exp(± sj \(T) -Ez} = 1 ± * ^ X(T)-E + 2(X(T)-E). (40)

After calculating the set of solutions E1,..., Et of the equation (37), i.e., all values of the spectral parameter E for which detM(E^) turns into zero, for each Es, s = 1,...,i we get the degenerate algebraic system of equations (32), linear and homogeneous with respect to the unknown coefficients C1;s, ...,C2Nk-s. To calculate the desired eigenvectors Cs = (C1;s, ...,C2Nk;s)T, we add an additional condition, e.g., C2N+1;s+-. +^2N(k-2);s =1. As a result, the inhomogeneous system of algebraic equations obtained for each Es, s = 1,...,i has the unique solution Cs = (C1;s, ...,C2Nk;s)T. Substituting it into the left-hand side of the normalisation condition (7) or (10), we get the expressions for the

normalising coefficients BS2, s = 1, ...,t:

k N

EE Î dz]\<Ô{2m-2)N+l-Am)(E1 -z) + C{2m-l)N+i E^ (E, z)

^-1 _1 J

ra=1 ¿=1

= Bl

Zm — 1

for bound states and

k N

ÊË / dz (A 2™-2)N+i;s^M(£, 2m-1)N+F^m)(E, z)) ' = B2g

m=1 i=1 J

— 1

for metastable states.

Table 1

M

The eigenvalues of bound E{ and metastable Ef = Re Ef + i Im Ef , Im < 0 states of the BVP (1) with effective potentials (43) at N = 4: (*) calculations by algorithm of MMFS. The same results are given using FEM with 3rd-type BC in intervals z e (-2, 2) and z e (-8, 8). (I) and (II): the estimations SEi = SfEM - Ei of eigenvalues Ei calculations using FEM in interval z e (-8, 8) with Dirichlet (I) and

Neumann (II) boundary conditions

m

2

m

(*) E1 E2 -0.3260959460 0.6358528584 ßM ßM 2.8027364541 3.7392597199 -0.0000002130« -0.0263330505«

(I) (*)+1.69e-7 (*)+2.09e-4

(II) (*)-1.69e-7 (*)-2.10e-4

As a result, we calculate the desired coefficients Cs : Ci;s = Ci;s/Bs, i = 1,..., 2Nk, s = 1,..., t, providing the fulfilment of the normalisation condition (7) or (10) for the desired orthogonal set of eigenfunctions, = 1, ... , :

$s( z) = {$S1)( z), Z < Z1, < ^2,..., z < Zk-1, &Sk),z> Zk-l}. (41)

3. Benchmark Calculations of Bound, Metastable and

Scattering States

For example let us consider the problem of bound, metastable and scattering states for the Schrodinger equation in the 2D domain Qyz = {y £ (°,n),z £ (-to, +to)},

( d2 d2 \

(-^2 - M (y,z) -E) *(y, z) = °, *) = *(*, *) = °, (42)

with piecewise constant potential

V (y, z) = {°,z < -2; - y, Izl < 2;2,z > 2} ,

presented in Fig. 2. We seek the solution of Eq. (42) in the form of expansion y, z)= Bi(y)$i(z) over the set of basis functions Bi(y)=sin(«y), which leads

to the system of ODEs (1) with the matrix of effective potentials (2) given in the form

Vij = {i2Sij,z « -2; i2Sij - t y^L sin(^)^L sin(jy)dy, z « 2; (i2 + 2)z > 2}. (43)

J v^ v^

Figure 2. 2D potential and some diagonal (solid curves) and nondiagonal (dashed

curves) effective potentials (44)

This system has two bound and a set of metastable states. The eigenvalues were calculated at N = 4, with matrix of potentials

Vij(z)--

1 0 0 0

0 4 0 0

0 0 9 0

\0 0 0 16/

/1--1 2

16 90 32

1- 0 9 -

4--4 2

48

25-

48 259—-9 2

96

32 2250

96 49-

\ 32 n 96 i g _-

\ 225- 0 49- 16 2/

3 0 0 0

0 6 0 0

0 0 11 0

0 0 0 18

< -2,

kl « 2,

> 2,

(44)

using the algorithm implementing the MMFS in the Maple system, and presented in Table 1.

So, the parameters A(m) and matrices of general solutions (3) calculated by

Step 1 by solving algebraic eigenvalue problem (4) take the form

A(1) = diag ( 1 4 9 16), ^(1) =

( 1 0 0 0 \ 0100 0010 0001

(45)

A(2) = diag ( -0.675679 2.459337 14.484720 7.448436 ) / -0.983202 0.182314 -0.003169 -0.008059 \

*(2) =

0.181948 -0.013991 0.003524

Ap = diag ( 3 6 11 18), ^(3) =

0.975788 -0.004632 -0.121288 -0.120697 -0.088442 -0.988642 0.005598 -0.996065 -0.088372

/ 1 0 0 0 \ 0100

0010 0001

(46)

(47)

(48)

From here we have desirable restrictions of Step 3.2 for the bound-state eigenvalues

-0.675679 = min( A

(2) 1,

..,Af})<E< min(A1

(1)

A1)

Af,

) = 1. (49)

These values were used to test the algorithm and programm KANTBP 4M, implementing the FEM in the Maple system [4]. The FEM calculations were performed with the Hermite elements of the seventh order (j/=7, p=3, nmax=2) on the finite element grid fiz = {-8(4) — 2(4)2(4)8}, where in the parentheses the number of elements is indicated. The calculations of the bound states with the Neumann and Dirichlet conditions yield

FEM —7

the upper and lower estimates of 5Ei = Ef — Ei with the accuracy 5E1 & 1.7 ■ 10 7

and 5E2 & 2.1 ■ 10-4. Note, the results, obtained by MMFS and FEM with 3rd-type BC

(III) on both intervals z G (—2, 2) and z G (—8,8) coincide with an accuracy of the order

10-10. In this table we show also upper and lower estimations of eigenvalues of bound

states obtained by FEM on interval z G (—8, 8) with Dirichlet (I) and Neumann (II) BC,

respectively.

The results of calculations of eigenvalues of bound and metastable states and the matrix elements of reflection and transmission amplitudes for the multichannel scattering problem, obtained by matching the fundamental solutions and using FEM on the grid fiZ = {—2(4)2} with the boundary conditions of the third kind and the asymptotic solutions of Step 2, coincide with the accuracy to the order of 10-10. Eigenfunctions y, z) of bound and metastable states of the 2D boundary-value problem and the components of the eigenfunctions ( z) (the imaginary parts shown by dotted lines) of the corresponding system of ODEs (1) are shown in Fig. 3.

Figure 3. Eigenfunctions ^(y, z) of bound and metastable states of the 2D boundary-value problem and the components of the eigenfunctions (z) (the imaginary parts shown by dotted lines) of the corresponding system of ODEs (1)

The considered system has the set of threshold energies E: E^ =1, E[n =3, Epn=4, Epn=6, ... that are different for the left- and right-hand asymptotic regions of axis. At given E £ (1, 3] we have only one open channel N% = 1, N^ = 0 and we have only incident and reflected waves in left-hand side. The S matrix (13) contains

only one element S = (En). At given E £ (3,4] and E £ (4, 6] for the wave incident from the right there is one open channel N^ = 1 and for wave incident from the left the there are one and two open channels N% = 1 and N% = 2, respectively. The dependence of the elements of scattering S matrix (13) is shown in Fig. 4.

E E E

Figure 4. The dependence of real and image parts, and the square of absolute values of elements of S—matrix (13) versus the scattering energy E. The threshold

• 771th; R o j TTfth; L *

energies ^ =3 and b2 =4

For example at E = 2.005, E = 3.305 and E = 4.805 the S-matrix takes form

S( E = 2.005) = (0.21+0.98« ), S(E = 3.305) =

0.53 0.32

—0.57—0.53«

—0.57—0

0.36+0

.53 .51

S( E = 4.805)

-0.14-0.17« 0.10-0.28« 0.10-0.28 -0.72+0.52

-0.88+0.29« 0.07+0.34«

-0.88+0.29 0.07+0.34

0.06 0.13

where the R^, R^, T^ and T^ submatrices are separated by lines. Corresponding typical wave functions of the multichannel scattering problems are presented in Fig. 5.

LR(1): E=2.005

LA

Op

ai - e -V A. 2 yf 6 8

-p.5 I / z

J-1

VA 1.5

RL(1):E=3.305

'A :

-8/-& V.1 /-2-.W 4 6; k

J ' XJ z ; 1

-1

-1.5

Figure 5. The real (solid lines) and the imaginary (dotted lines) components of the solution of the scattering problem for the wave incident from the right, RL(1), and the waves, incident from the left from the first, LR(1), and the second LR(2)

open channels

One can see from Fig. 4 we have the resonance behavior of transmission and reflection amplitudes near E « 3.7 that corresponds to the resonance state at E = 3.748 presented in Fig. 6 while the other resonance state with the energy in vicinity of E = 2.803 is not resolved in this scale. The scattering energy E = 2.803 and E = 3.748 corresponds to real part of energy Re Ef = 2.803 and Re E2f = 3.739 of metastable states (see Table 1), while half-width of peaks corresponds to imaginary parts of the energy Im Ef1 = -2.1 ■ 10-7 and Im Erf = -0.026 of these metastable states. So, the wave functions of metastable states in Fig. 3 have similar behavior with wave functions of resonance scattering states in Fig. 6.

LR(1): E=2.8027364541 LR(1): E=3.748 RL(1):E=3.748

Figure 6. The the same as Fig. 5 but for resonance solutions of multichannel scattering problem that will be compared with the eigenfunctions of the metastable states in Fig. 3 and Table 1

So, the elements of S-matrix at different values of energy E, obtained by MMFS and FEM with 3rd-type BC (III) on both intervals z e (-2,2) and z e (-8, 8) coincide with an accuracy of the order 10-10.

4. Conclusion

Effective Symbolic-Numerical Algorithms of Matching the Fundamental Solutions is elaborated and implemented in Maple. At large number k of subintervals or large number N of equations we have a cumbersome secular equation for bound and metastable states solving by author's algorithm, which will be published in a more detail elsewhere. It was shown that the results of benchmark calculations obtained by algorithms of Matching the Fundamental Solutions at each boundary between the adjacent axis subintervals and Finite Element Method in finite interval with third type boundary conditions coincided with an accuracy of the order of 10-10. The proposed algorithms and programs can be applied for solving boundary-value problems for systems of N second-order ordinary differential equations (ODEs) of the waveguide type with the matrix of piecewise constant potentials arisen in mathematical modelling of optical multilayer systems and quantum-dimensional nanostructures.

References

1. C.-C. Huang, C.-C. Huang, J.-Y. Yang, An Efficient Method for Computing Optical Waveguides With Discontinuous Refractive Index Profiles Using Spectral Collocation Method With Domain Decomposition, Journal of Lightwave Technology 21 (10) (2003) 2284.

2. A. L. Sevastyanov, L. A. Sevastianov, A. A. Tyutyunnik, Analytical Calculations of Derivation Partial Differential Equations for Coefficient Kantorovich Functions, Matem. Mod. 27 (2) (2015) 103.

3. M. N. Gevorkyan, D. S. Kulyabov, K. P. Lovetskiy, A. L. Sevastyanov, L. A. Sev-astyanov, Waveguide Modes of a Planar Optical Waveguide, Mathematical Modelling and Geometry 3 (1) (2015) 43.

4. A. A. Gusev, L. L. Hai, O. Chuluunbaatar, S. I. Vinitsky, Program KANTBP 4M for Solving Boundary-Value Problems for Systems of Ordinary Differential Equations of the Second Order.

URL http://wwwinfo.jinr.ru/programs/jinrlib/kantbp4m

5. A. A. Gusev, L. L. Hai, O. Chuluunbaatar, V. Ulziibayar, S. I. Vinitsky, V. L. Derbov, A. GoZdZ, V. A. Rostovtsev, Symbolic-Numeric Solution of Boundary-Value Problems for the Schrödinger Equation using the Finite Element Method: Scattering Problem and Resonance States, Lecture Notes in Computer Science 9301 (2015) 182.

УДК 519.632.4

Алгоритмы и программы решения краевых задач для систем ОДУ второго порядка с кусочно-постоянными потенциалами: многоканальная задача рассеяния и задача на собственные

значения

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

* Объединённый институт ядерных исследований, г. Дубна

^ Российский университет дружбы народов, г. Москва

* Белгородский государственный университет, г. Белгород S Саратовский государственный университет, г. Саратов

32 Институт физики, университет им. М. Кюри-Склодовска, г. Люблин, Польша

Предложены новые алгоритмы и программы, реализованные в системе Maple для решения многоканальной задачи рассеяния и задачи на собственные значения волноводно-го типа для систем ОДУ второго порядка с матрицей кусочно-постоянных коэффициентов размерностью N х N на оси. Разработаны новые алгоритм и программа для решения краевой задачи методом сшивки фундаментальных решений (МСФР) системы ОДУ в точках разрыва потенциалов. На каждом из подынтервалов оси общее решение системы ОДУ ищется в виде линейной комбинации 2N фундаментальные решений с неизвестными коэффициентами. Каждое фундаментальное решение явно зависит от спектрального параметра и собственных значений и собственных векторов алгебраических задач на собственные значения с матрицей постоянных потенциалов размерностью N х N. Из условия непрерывности решений и их производных в точках разрывов потенциалов следует система алгебраических уравнений. В случае задачи на связанные или метастабильные состояния полученная система алгебраических уравнений содержит нелинейную зависимость от неизвестного спектрального параметра. Для решения такой нелинейной задачи сформулирован символьно-численный алгоритм. Дано сравнение эталонных расчётов связанных, метастабильных состояний и состояний рассеяния краевых задач для систем ОДУ второго порядка, выполненных с помощью программ, реализующих алгоритмы МСФР и метода конечных элементов.

Ключевые слова: многоканальная задача рассеяния, задача на собственные значения, система ОДУ второго порядка, методом сшивки фундаментальных решений

Литература

1. Huang C.-C., Huang C.-C., Yang J.-Y. An Efficient Method for Computing Optical Waveguides With Discontinuous Refractive Index Profiles Using Spectral Collocation Method With Domain Decomposition // Journal of Lightwave Technology. — 2003. — Vol. 21, No 10. — P. 2284.

2. Sevastyanov A. L., Sevastianov L. A., Tyutyunnik A. A. Analytical Calculations of Derivation Partial Differential Equations for Coefficient Kantorovich Functions // Matem. Mod. — 2015. — Vol. 27, No 2. — P. 103.

3. Waveguide Modes of a Planar Optical Waveguide / M. N. Gevorkyan, D. S. Kulyabov, K. P. Lovetskiy, A. L. Sevastyanov, L. A. Sevastyanov // Mathematical Modelling and Geometry. — 2015. — Vol. 3, No 1. — P. 43.

4. Gusev A. A., Hai L. L., Chuluunbaatar O., Vinitsky S. I. Program KANTBP 4M for Solving Boundary-Value Problems for Systems of Ordinary Differential Equations of the Second Order. — http://wwwinfo.jinr.ru/programs/jinrlib/kantbp4m.

5. Symbolic-Numeric Solution of Boundary-Value Problems for the Schrödinger Equation using the Finite Element Method: Scattering Problem and Resonance States / A. A. Gusev, L. L. Hai, O. Chuluunbaatar, V. Ulziibayar, S. I. Vinitsky, V. L. Der-bov, A. GoZdZ, V. A. Rostovtsev // Lecture Notes in Computer Science. — 2015. — Vol. 9301. — P. 182.

© Gusev A.A., Chuluunbaatar O., Vinitsky S.I., Hai L. L., Derbov V.L., GoZdZ A.,

2016