CC BY
27
UDC 530.145
WKB-BASED SCHEMES FOR TWO-BAND SCHRODINGER EQUATIONS IN THE HIGHLY OSCILLATORY REGIME
J. Geier1, A. Arnold1
1 Institute for Analysis and Scientific Computing, Vienna University of Technology, Wien, Austria
[email protected], [email protected] PACS 02.60.Lj, 02.60.Lh, 85.30.De
An efficient and accurate numerical method is presented for the solution of highly oscillatory differential equations in one spatial dimension. While standard methods would require a very fine grid to resolve the oscillations, the presented approach uses first an analytic WKB-type transformation, which filters out the dominant oscillations. The resulting ODE-system is much smoother and can hence be discretized on a much coarser grid, with significantly reduced numerical costs.
Here we are concerned with stationary two-band Schrodinger equations employed in quantum transport applications. We focus on the Kane-model and the two band k ■ p-model. The accuracy of the presented method is illustrated on a numerical example.
Keywords: Schrodinger equation, Kane-model, two-band k ■ p-model, highly oscillating wave functions, higher order WKB-approximation, asymptotically correct finite difference scheme.
1. Introduction
This paper deals with an asymptotically correct scheme for the numerical solution of stationary, vector valued Schrodinger equations in the highly oscillatory regime. To fix the ideas we first recall this procedure from [2], where the scalar Schrodinger equation in 1D is dealt with. We consider highly oscillating differential equations of the type
e2^'(x) + a(x)^(x) = 0 , (1)
where 0 < e ^ 1 is a very small parameter and a(x) > a0 > 0 a sufficiently smooth function. For very small e > 0, the wave length A = ~¥= is very small, such that the solution tp becomes
V a-(z)
highly oscillating. In a classical ODE-scheme such a situation requires a very fine mesh in order to accurately resolve the oscillations. Hence, standard numerical methods would be very costly and inefficient here. The goal of the numerical WKB-scheme is to provide a method that uses a coarse spatial grid with step size h > A.
The strategy in [2] is to first analytically rewrite (1) as a "smoother" (i.e. less oscillatory) problem that is then easier to handle numerically. In this analytic preprocessing, the dominant oscillations are eliminated. This transformation is closely related to using a second order WKB-approximation for tp:
V{x) « C exP (^)) , fa) := f (y^) _ dT , (2)
V a(x) J o v 7
with /3 :=-^(a-V*)" .
Standard numerical methods for (1) (like in [8]) used to require a step size h = O(e) for accurate resolution. The WKB-based schemes from [2] or [18] reduce that limitation to h = 0(y/e) (when using a forth order quadrature rule for the integration of the phase 4> in (2)). If this phase is explicitly integrable (for piecewise linear a(x), e.g.), the mesh size h can even be chosen independently of e in the scheme from [2]. This scheme is then asymptotically correct w.r.t. £.
The goal of this paper (and of [5]) is to generalize this procedure to vector valued Schrodinger equations in one spatial dimension. Here we shall consider two-band Schrodinger systems that are employed in quantum transport models, e.g. for simulating resonant interband tunneling diodes (see [3,20]). In particular we shall focus on the Kane-model [13,14] and the k ■ p-model [19].
We remark that numerical integrators for oscillatory ODE-systems have been addressed by several authors in the last decade (cf. Sec. 5 of [9], [11,18], Sec. XIV of [6]). Those methods are closely related to using the zeroth order WKB-approximation, tp{x) & C exp J(j' y/adr^j to eliminate the dominant oscillations in (1). But the higher order transformations employed below provide a refined asymptotic, with error estimates of higher order in
This paper is organized as follows. In Sec. 2 we introduce the Kane-model and k ■ p-model, along with their open boundary conditions that are needed in quantum transport applications. These boundary value problems (BVPs) are then transformed into equivalent initial value problems (IVPs). In Sec. 3 we make an analytic preprocessing of these oscillatory IVPs. This transforms them into "smoother" ODEs. In Sec. 4 we develop for these ODEs a second order, asymptotically correct numerical scheme. The error of this method is of the order 0(mm(eh2,£3)). To achieve this, particular care has to be taken for the discretization of the involved oscillatory integrals. In Sec. 5 we illustrate the effectiveness of our method in a numerical example for the Kane-model and on an example from the literature [18].
2. Two-band Schrodinger models
For several novel semiconductor devices (like resonant interband tunneling diodes, see [3,20]) single-band effective mass models become insufficient to simulate the quantum transport through such a device. Hence, it is getting ever more important to include realistic band structures in quantum transport models. In this section we shall discuss two independent, stationary two-band Schrodinger models (Kane-model and two-band k ■ p-model) that are used for simulations of one dimensional semiconductor devices like a resonant tunnel diode. We assume that the considered semiconductor structure occupies the interval [a, b] and is connected to reservoirs at x = a and x = b. We also assume that both reservoirs are homogeneous and extend to x = So all material (and energy) parameters are constant in each reservoir, which is hence populated only by traveling plane waves.
In analogy to (1) we shall discuss in this paper only the numerically challenging oscillatory regime of traveling waves. The evanescent regime with tunneling is of course equally important for quantum transport, but the numerical treatment of those smooth wave functions is simple, and will not be discussed here. However, WKB-based discretizations of the coupled oscillatory-evanescent situation are currently under investigation, and will be the topic of an upcoming work. In [17], transparent boundary conditions (TBCs) for the Kane-model and k ■ p-model were derived, as well as discrete TBCs for finite difference schemes. However, such classical schemes are numerically expensive in the oscillatory case. So it is the goal of this paper to develop a more efficient alternative.
2.1. The two-band Kane-model
The simplest multi-band approach is the two-band Kane-model (cf. [15]). This is an inter-band model, describing the coupled electron transport in the conduction and the valence bands. Here the "wave function" ^(x) e C2 is a solution of the ODE
H^ = E^, (3)
with
H
/ V{x) -iep{x)j-x \ \-isp(x)j-x V(x)-EB(x)J
We denote by E > 0 the prescribed (constant-in-x) injection energy of the electrons, and e > 0 is a small constant. Here, the potential V is the band edge of the conduction band, and Ea > 0 is the band gap between the conduction and the valence band. The function p > 0 is related to the Kane-parameter. The dispersion relation of this Kane model is discussed in detail in Sec. 3.1 of [17].
In order to have unique solutions we assume:
Assumption 1. The functions V,Eg ,p are piecewise1 Lipschitz continuous on the non-empty bounded interval [a, b] and constant on the exterior domains (—m, a] and [b, <x>). Further, there exist open neighborhoods of a,b where they are continuous.
We rewrite the ODE (3) in the more convenient form
,,, , if 0 a(x) \
^ = e[m 0 )
with
E-V(x)+Eg(x) E - V(x)
a(x) :=-^——r-and pi:x) :=-. , .
p(x) p(x)
We shall consider here a scattering model, subject to a given, incoming plane wave. Hence, we shall need transparent boundary conditions (TBCs) at both (artificial) boundary points a, b (as derived in Sec. 3.2 of [17]). We denote the system matrix of (4) by A(x), with the eigenvalues
A(x) = ±-y/a(x)/3(x)
= (E ~ V(x) + Eg{x)){E - V{x)) =: ±ik(x). (5)
If the given injection energy E is larger than the conduction band edge, i. e. if E — V(x) > 0 holds on the whole interval [a, b] (and thus on the whole real line), the eigenvalues A = ±ik are purely imaginary and hence we only have traveling waves2. Let v±(x) be (real valued) eigenvectors of A(x) corresponding to the eigenvalues A(x) = ±ik(x). From the right exterior
1Here, piecewise Lipschitz continuous means that there are finitely many (non trivial) pairwise disjoint intervals 11 ,...,In such that [a, b] c |J"=1 Ij and the functions are globally Lipschitz continuous on each single interval Ii,..., In.
2One also gets purely imaginary eigenvalues if the energy is smaller than the valence band edge, i. e. E — V(x) + Eg (x) < 0. This energy corresponds to holes in the valence band, and the situation would be analogous to the case discussed here.
domain [b, m) we now inject a left traveling electron wave with prescribed amplitude db = 0. Then the resulting scattering state has the following form in the exterior domains:
(x) = dae_lk(a)(x_a)v-(a) for x < a ,
M%) = cbelk(b)(x-b)v+(b) + dbe_lk(b)(x_b)v-(b) for x > b, (6)
with constants cb,da,db E C. We denote the solution in the interior domain [a,b] by ty. Due to Assumption 1 the solution on the whole real line is continuously differentiable in certain open neighborhoods of the boundary points a, b. Hence we get the (homogeneous) left TBC
= (a) E span[i—(a)] ^ {A(a) + ik(a) Id)= 0
^ ty'(a) + ik(a)ty(a)=0 .
At the right boundary we combine the first derivative ty' with ty to eliminate the reflection constant cb. From (6) we obtain the (inhomogeneous) right TBC
ty'(b) - ik(b)ty(b) = -2ik(b) db v-(b). The resulting (inhomogeneous) BVP then reads
ty'(x) - A(x)ty(x) = 0, x E [a,b] , (7)
ty' (a) + ik(a)Idty(a) = 0 (8)
ty'(b) - ik(b)Idty(b) = -2ik(b)dbv-(b). (9)
Its unique solvability was established in Sec. 3.3 of [17]. Thus we state without proof:
Lemma 2.1. The BVP (7)-(9) has a unique solution in (W2'^(a, b))2.
Remark 2.2. Since ty satisfies the ODE (7) even on the boundary, the TBCs can be reformulated as
(A(a) + ik(a)Id)= 0 ,
(A(b) - ik(b)Id)ty(b) = -2ik(b)dbv-(b). (10)
We shall now reformulate the BVP (7)-(9) as an IVP, which is easier to solve numerically (particularly for our highly oscillatory regime). To this end we first solve (using the left boundary condition) the IVP
4>'_(x) = A(x)ty_(x) , x E [a,b]; ty_(a) = v-(a) E C2 . (11)
Since ty has to fulfill (8) (which is equivalent to E span[i—(a)] ), there exists a constant c_ E R such that = c_v_(a). Hence we get ty = . Using the remaining boundary condition (10) we get
c_ (.A(b) - ik(b) Id)ty_(b) = -2ik(b)dbv-(b). The inner product of this equation with v_ (b) yields
-2%k{b)db\\v.{b)\\2 _ -2ik(b)db\\v-(b)\\2 C~ ~ v_(b)T(A(b) - ik(b) Id)V-(&) ~ v_(b)T(4>'_(b) - ik(b)ip_(b)) '
Analogously, we get for a right traveling plane wave in the left exterior domain (-ro,a] with prescribed amplitude ca = 0:
ty'(x) - A(x)ty(x) = 0 ty' (a) + ik(a)Idty(a) = 2ik(a)cav+ (a) ty'(b) - ik(b)Idty(b) = 0 .
It holds ^ = c+V>+, where ^+ solves
(x) = A(x)^+(x), Mb)= v+(b) e C2
and
c+
2ik(a)ca ||i>+(a)|
2ik(a)ca ||i>+(a)|
v+ (a)T(A(a) + ik(a) Id)^+(a) v+(a)T(-0+ (a) + ik(a)vp+ (a))
Recall from (4) that the system matrix A(x) is proportional to 1/e and off-diagonal. We now aim to transform out the dominant oscillations in the IVP (11). To this end we first diagonalize A as A(x) = 1Q(x)~1L(x)Q(x) with
Q(x)-1 :=
(
_i _ ( xAM ~Va(x)
vW) vW)
)
and L(x)
(
ek(x) 0 0 -ek(x)
Note that the matrix valued functions Q and L are real valued and actually independent of e. The ansatz u := and Qv- = (0, 1)T yields for x E [a, b]:
u' = -Lu + Bu , u{b) =
e
(12)
with B := Q'Q-X. The same transformation also works for the other case of a right traveling plane wave with prescribed amplitude in the exterior domain (-m, a]. We only have to replace the initial condition by u(a) = (1,0)T.
2.2. The two-band k ■ p-model
In this section we discuss a slightly more involved inter-band model for the coupled quantum transport in the conduction and the valence bands (cf. [3,12]). A different inter-band k ■ p-model is analyzed in Sec. 4 of [17]. And for extended multi-band k ■ p-models (including the intra-band coupling of heavy and light holes within the valence band) we refer to [13], Sec. 6 of [3], and Sec. 5 of [17]. In our two-band model, the "wave function" ^ = tyiE C2 solves a 2 x 2 Schrodinger boundary value problem
H(x)^(x) $ (a) - Ka(E )^(a) ^(b) - Kb(E)^(b)
E^(x) 0
re C2 ,
x e (a, b)
(13)
(14)
(15)
with the Hamiltonian
H
and
— £
dx2
0
V (X) - Eg (X)
P (x)
( 0 ip(x) \ y ip(x) 0 J
The real parameter e > 0 is a small constant, which is often (depending on the model) proportional to the reduced Planck constant h. By E we denote the given injection energy of the particles. The potential V(x) is the band edge of the conduction band, and Eg (x) > 0 is the energy gap between the conduction and valence bands. Further, p(x) > 0 is the coupling coefficient (related to the Kane-parameter) between the two bands. As in [3] we assume p(x)2 > Eg(x) Vx E [a,b], which implies that the valence band is concave close to k = 0. We remark that (13)-(15) is a scattering problem with given E, and not an eigenvalue problem.
The matrices Ka, Kb E C2x2 and the vector r in (14), (15) constitute the TBCs for the k ■ p-model. (15) models the injection of plane waves at x = b. Their derivation follows the same strategy as for the Kane model in Sec. 2.1. But for the more involved details we refer to [3]. The self-consistent potential V is the solution of the following Poisson problem:
V" (x) = n(x) , x E (a,b), V(a) = Vi > 0 , V(b) = 0 . The charge density n is defined by
/><x
n(x) = f (E) \ip(x,E)|2dE,
J 0
for a prescribed function f that models the semiconductor statistics. Well-posedness of this nonlinear BVP is established in Th. 2.2 of [3].
If one is interested in a numerical approximation of n(x) or of the current density
/><x
J(x) = f (E)( - Im^',^) + 2p Re(^2))(x,E) dE , 0
one has to use an iterative scheme, like the Gummel method, to solve the nonlinear problem. In each step one has to calculate a suitable approximation for the charge density n, and hence one has to solve (13) for a large number of (discrete) energies. Since 0 < e ^ 1 is a small constant the wave function ^ is highly oscillatory for E > \\V^. In order to speed up the calculations it is very useful to have a numerical scheme that produces an accurate approximation for without having to resolve all oscillations of the wave function.
It is often more convenient to solve, instead of a BVP, an equivalent initial value problem. As we will see in a moment, it is possible to determine the solution of (13) from just four (vector valued) IVP-solutions. Let the functions p, V, and Eg be Lipschitz continuous on [a, b]. Hence the IVP for a matrix-valued wave function $(x) e C2x2,
H(x)$(x) = E$(x), x E (a, b), (16)
&(a) = Id, (17)
$'(a) = Ka, (18)
has a unique solution. Further, it holds for every vector valued solution 0 of
H(x)0(x) = E0(x) , x E (a, b), (19)
0(a) - Ka0(a) = 0 ,
that 0(x) = &(x)4>(a). Since the solution ^ of the BVP (13)-(15) solves (19), we can write ^(x) = &(x)^(a). Hence we get from the remaining right boundary condition (15)
(&(b) -Kb$(b))^(a) = r.
Since the BVP (13)-(15) is well-posed (cf. [3], Prop. 2.1), the above equation has a unique solution which yields ^(a) and consequently
As we have seen, the solution ^ of the BVP (13)-(15) is (uniquely) determined by the solution $ of the IVP (16)-(18). Thus, in the sequel we shall derive and discuss an efficient numerical method to solve the IVP for the (vector valued) equation (19). Equation (19) for 0 E C2 takes the form
e V -eP(x)$ + A(x)0 = 0 , (20)
with A(x) := diag(E — V(x),E — V(x) + Eg(x)). For E > \\Vthe matrix A(x) is positive definite for all x E [a, b]. Now we want to rewrite (20) as a first order IVP, with the same form as (12). To this end we set
which yields for v(x) = (vi(x), v2(x))T E C4:
£ V P J V 0 0 /
= №) •
£Ka
Let us denote by L the first matrix of (21). Since P(x) is skew-symmetric for all x E [a, b], the same holds for L. Hence there exists a matrix function Q: [a, b] ^ c4x4, such that for all x E [a, b] it holds (cf. the Appendix Sec. 6)
L(x) = iQ*(x)L(x)Q(x),
with L(x) real and diagonal. Finally we set
u(x) := Q(x)v(x) E C4 ,
which yields
u' = -Lu + Bu , u{a) = Uq , (22)
with
B(x) = Q'(x)Q*(x) + Q(x)^Ah(-x^~2W °)Q*0r).
Of course, the above transformation procedure is not limited to the special case (13). One can apply it to any ODE of type (20) with P(x) skew-symmetric and A(x) positive definite. Hence we shall continue our discussion for general equations of the form (22).
Since L is diagonal and real valued, the solution u is highly oscillatory. But \\u\\ is bounded by a constant independent of . In order to prove this, we introduce a smoother, "adiabatic" variable r] (which coincides with the "r/" from [11,18]). This can be interpreted as the lowest order WKB-type transformation (see Sec. 3 for details).
In the sequel we denote by \\-\\ the Euclidean norm on Cra and also its subordinated matrix norm. Further we define for continuous vector valued functions f: [a, b] ^ Cd the w-norm by
\\f\L = sup \\f(x)\\ ,
x£[a,b]
and analogously for matrix valued functions. With these definitions we can establish Lemma 2.3. Let u be the unique solution of the IVP (22), and let
E£(x) := exp ^ jT L(s) ds^j E C4 .
Then the new quantity rj(x) := E*(x) u(x) E C4 solves the IVP
7]' = (E*BE£) 'q , 'q(xo) = ^ := uo . There exists a constant c > 0 independent of e such that it holds for all x E I:
\\u(x)\\ = \\v(x)\<c.
Proof. Since L(x) is diagonal and real valued, E£(x) is unitary. Differentiation of the ansatz ^ = E*u yields the IVP. A standard Gronwall argument yields
\\V(x)\\ < e\x-a^BE^ H^ll .
Since E£ is unitary, it follows
\\«(®)H = \\v(x)\\< elx-ami- M .
Since B does not depend on e, the desired constant is c = e1'3-'1^^^ . □
3. Analytic preprocessing
The imaginary part of the system matrix of the IVPs (12) and (22) (i.e. ^L) causes oscillations of the solution u, whose wave length is proportional to e. Thus, if e is very small, the function u is highly oscillatory. Hence, directly solving these IVPs with standard integrators is computationally very costly.
In order to devise an efficient numerical method, we first apply an analytic preprocessing to the IVPs. While the system matrix in (22) has the leading order 0(1/e), the transformed problem will have the order O(e). This does not mean that the oscillations disappear. In fact, one can observe that the frequencies are even (approximately) doubled. But the amplitude of the oscillations in the new variable y will only be O(e), down from the 0(1) oscillations in the original variable u (cf. Lemma 2.3). Clearly, we can expect that the new "smoother" IVP for y is much easier to solve numerically than that for u, particularly in the regime of small e.
The following discussion focusses on the IVPs (12) and (22). But it can be extended to a larger class of problems and refined to higher e-order asymptotics (see Corollary 3.3, and for more details [5]). We formulate the essential requirements on the functions diag(Al; ...,\d) = L from both IVPs in
Assumption 2. There exists a 5 > 0 such that for all x E [a,b] and i = j it holds
|A,(x) - A,(x)| > 8. I. e. L(x) has only simple eigenvalues, which stay separated over the whole interval [a, b]. For the transformation we make the following ansatz
y(x) := E(x)T-l(x) u(x) , (23)
where we set
Ee(x) := exp (j-J^ L(s) ds^J , (24)
T£(x) := To(x)+ eTi(x). (25)
Since L(x) is a real valued diagonal matrix for all x E [a,b], the matrix E£(x) is unitary and diagonal, i. e. it holds
E-\x) = E*(x) = exp ( - - I L(s) ds
Before we can write down the matrix Tl we need some notation. Let A, B E Cdxd. The entry-wise product (Hadamardproduct) of A and B is defined by
A © B := (aijbij)i<i,<d.
Moreover, for a diagonal matrix L = diag(Ai, . ..,Ad) with pairwise distinct simple eigenvalues Ai,..., A<i, we define the following two matrices:
Dl := (Xi - Xj)i<itj<d; D~l := (y^A , (Dl) := 0. (26)
We remark that it holds V A E Cdxd:
LA — AL = Dl © A and DL © D- © A = A — diag A. (27)
With this notation we can establish
Lemma 3.1. We make Assumption 2. Let the matrix valued functions L,B E Cr ([a, b], Cdxd) for some r E N. Further, let T0 be the unique solution of the (matrix valued) IVP
TO(x) = diag(B(x))To(x) , To(a) = Id E Cdxd, (28)
and let
Ti(x) := iD~l{x) © (B(x)To(x)).
Then there exists a constant £0 > 0 such that: For all 0 < e < £0, the matrix T£(x) is regular for all x E [a, b], and the variable y defined in (23) solves the IVP
y' = £E*sSiE£y, y(a) = yo. (29)
The matrix valued function S1 is given by
Si := T-1(BTi — Ti) ,
and hence is of class Cr-i([a, b], Cdxd).
Proof Clearly, the IVP (28) has a unique solution To E Cr+i([a, b], Cdxd). By Assumption 2,
D-{x) and Ti(x) are well defined for all x E [a, b], with D-{x), Ti(x) E Cr([a, b], Cdxd).
**dxd\
LU DL(X), Ti(x) E C Ua' u]> C
Since To(x) is regular for all x, we can write
T£(x) = To(x)( Id+£To(x)-iTi(x)).
We set
1
£o := mm
xe[aM \\To(x)LiTi(x)\\ •
With the von Neumann series we immediately get that T£(x) is regular for all 0 < e < eo, and
T£ E Cr([a, b], Cdxd). Differentiating (23) yields
y' = E:[l(T-1L-LT-1)+T-1' + T-1B](E:T-1)-1y.
Using T-iL — LT-i = T£l1[L, T£]T£l1 and 0 = T£l1'T£ + T-iT"s we get
y' = E*eT~l T£] - T'£ + BTe)Eey .
With T£ = To + eTi we deduce
y' = E*eT~l(^[L, To] + i[L, Ti] + BT0 — Tq + eBTx - £T[)Eey . (30)
Since To is diagonal the first commutator in (30) is zero. Using (27) we find
t[L,Ti] = IDl © Ti
= —Dl ©D- © (BTo) = —(BTo — diag(BTo)) ,
which yields i[L, Ti] + BTQ - TO = 0. Hence (30) reduces to
y' = sE*£T~l(BTi - Ti)E£V.
B,Ti E Cr([a,b], Cdxd) implies Si = T~l(BTl - Ti) E Cr-i([a,b], €dxd) which completes the proof. □
Corollary 3.2. Let the assumptions of Lemma 3.1 hold. Then there exists a constant c > 0, independent of e, such that
(x) - yo\\ < c£ and \\y'(x)\\ < c£ hold for all x G [a, b], 0 < e < y. Proof. Integration of (29) yields
r-x r-x
\\y(x) - yo\\ < e W(E:SiE£)(s)W\\y(s) - yo\\ ds + e \\(E*£\ \ ds .
J a J a
Since E£(s) is unitary for all s E [a,b] we obtain
r-x
(x) - yo\\ < e (snUy^) - V0\\ ds + elb - «US\L \M -
A Gronwall argument implies
(x) - yo\\ < elb - alUS\\00\\yo\\ e£ №
which yields the first estimate, since T£ 1 (and hence also Si) is uniformly bounded for 0 < e < e0/2 and x E [a, b]. The second estimate then follows from the ODE (29). □
For the numerical scheme it turns out that one can achieve better results (i.e. higher order e-asymptotics of the error) if the diagonal of the system matrix in (29) would be zero.
Corollary 3.3. Let the assumptions of Lemma 3.1 hold and let (the diagonal matrix function) R: [a, b] ^ <Cdxd be the unique solution of the IVP
R' = e diag(Si) R, R(a) = Id . (31)
If y is the solution of (29), then z := R-1 y solves
z' = £ (E*e SE£) z , z(a) = z0 (32)
with S := R-1 (Si - diag(Si))R. Furthermore, z admits the improved estimates
\\z(x) - Z0\\< c £2 , \\z' (X)H<C£,
with a constant c > 0 independent of e.
Proof. See [5]. □
In typical quantum transport applications, the coefficient functions often have jump discontinuities. This can be dealt with by restarting a new IVP at the x-point of discontinuity.
4. Numerical scheme
In this section we derive an explicit one step scheme with convergence order two for the IVP (32). To be more precise we shall derive local error estimates of order3 0(e 1hl), where hn is the local step size of the spatial grid. Let a = xQ < x1 < • •• < xN = b be a grid and let hn := xn+1 — xn be the local step size. Further we define the maximum step size by h := max0<n<N-1(xn+1 — xn). Integration of (32) over the subinterval [xn,xn+1] yields iteratively
fXn+1 fXn+1 i'S
z(xn+i) = z(xn) + £ / M(s) ds z(xn)+ e2 M(s) M(r)drds z(xn)
fXn + 1
+£ 3 / M (s) M (r) M (t)z(t)dtdrds,
(33)
where we use the abbreviation M := E*SEe e Cdxd. The last term is of order 0(£3hn) and, hence, can be neglected for the purpose of our error estimate. Let us define
I\x) :
M( )
E*(t)S (t)Ee(t)dt.
with x e [xn,xn+\] and I2(xn+\) :=
fXn + 1
M(x) M(t) dt dx
r-Xn+1
M(x)Il(x) dx.
(34)
(35)
So the first and second integral of (33) are I1 (xn+\) and I2(xn+\), respectively. Since they both have highly oscillatory integrands, we need specially designed quadratures in order to obtain sufficiently good approximations (for a review of highly oscillatory quadrature we refer to [10]). Here we shall use an ansatz, which can be found in [21] and also in [5].
Recall from Corollary 3.3 that the diagonal components of S are zero. So, we now consider the off-diagonal elements of I1. Its ¿jth-component reads
with
Neglecting the indices,
f j(x)
I(x) : =
Sijity-i*"® dt
[Ai(s) — \j (s)] ds
dt
(36)
(37)
is a prototype oscillatory integral for the off-diagonal elements of I1. Due to Assumption 2 we have \'Pij\> S > 0 for all x e [a, b]. Hence, we also assume for the phase in (37) |ip'\> 5 > 0 for all x e [a, b]. To derive an accurate quadrature for e small we shall leave the oscillatory factor unchanged and approximate only the e-independent part s of the integrand. We shall choose an approximation, such that the resulting integral is exactly computable. This is the case when replace s in (37) e.g. by p'pk, with some k e N. Hence, the space spanned by these generalized monomials is well suited to find an appropriate approximation of s. Now let
p(x) := p' (x)( Cq + cip(x)) (38)
3We say f: [0, h0] x (0, e0) ^ Cd is of order 0(enhm), iff there exists a constant c > 0 independent of e and h, such that \\f (h,e)\\ < cenhm for all (h,e) e [0,ho] x (0,eo).
n
u
u
S
n
n
n
X
X
u
n
X
u
u
u
X
n
X
n
be our ansatz function, with some c0,c1 E C.
In order to determine c0, cl we integrate by parts in (37):
I(x) = zs^-e-i^)
-%£ \ —7-z I e £
t=xn Jxn\ <f
(39)
Repeated integration by parts then yields the asymptotic expansion of I as discussed in [10]. For a moment let us replace s by s — p, which shall be the pointwise approximation error. Thus, if we require (s — p)(xn) = (s — p)(xn+\) = 0, the integral I(xn+\) is of order 0(ehn), provided
^tis well defined and bounded. And if (--/) is even continuously differentiate, a further integration by parts yields I(x) = 0(e2). Since this asymptotic behavior is desirable, we now choose p as the solution of the generalized interpolation problem
p(xn) = s(xn), p(xn+i) = s(xn+i).
This implies
Sn + 1 s.
n
-
fn + l fn 1 °n iACW
c 1 = —1- and Co = — - CiWn , (40)
<Pn+i - Vn fin
where we use the abbreviation sn = s(xn), and analogously for Lp,^. Replacing s by p in I yields the quadrature
r-x
I(x) « q(x) := / p(t) dt = i£e~^{t) (c0 + dip(t) - t£Cl) . (41)
JXn n
Lemma 4.1. Let E C2([a,b]) and let | fi| > 8 > 0 for all x E [a,b]. Then there exists a constant c > 0 independent of e and n, such that for all x E [xn, xn+\] C [a, b] it holds
ll(x) - q(x)l < ch2n min(e,hn) , (42)
ll(xn+i) - q(xn+i)l < chn min(hn-k£k) . (43)
k=0,2
Proof. Let us start with the integral
r-x
3{x) := I(x)-q{x) = / (s - p){t)e~^{t) dt. (44)
J Xn
By assumption, l^'(x)l> 8 > 0 for all x E [xn, xn+i\. Hence, is invertible (without restriction of generality let be increasing) and we substitute £ = tp(t) in (44), yielding
J{x) = / = / (g - d£. (45)
Jf(Xn) V V / Jf(Xn)
Here we set g(£) := (</?_1(£)) and 7r(£) := c0 + ci£. By construction, ir interpolates g at the points tp(xn), ip(xn+i). Since sE C2([a,b]), so is g. The theory of polynomial interpolation [16,22] thus yields the following error representation for £ e [p(xn), 'p(xn+i)]:
9(0-AO = ^(<(0)«-^))«-^)),
for some ((£) e [f(xn), '•p(xn+i)]. This yields the uniform estimate
|»(0 - ir(0 < ll^ll J^-f^f < cfcJ, ,46)
x
x
and c is independent of n. From [16] we get an estimate for the first derivative of the interpolation error. There exists a ^ e [p(xn), p(xn+1)], such that it holds for all £ e [p(xn), p(xn+1)\.
\g'(0 — n'(£)| < ||g(2)\\x — < ch.
Hence it holds for all x e [xn, xn+1] (using (45) and (46))
(47)
\J(x)\ < \p(xn+i) - p(xn)\ch2n < ch.
bn
Moreover, integration by parts in (45) gives
J(x) = i£(g -n)(£)e
"H
p(x)
rp(x)
- IE I (g-n)'(£)e
Thus, the estimate (42) follows from (47). Since n interpolates g at the boundary points p(xn), p(xn+1), we find (with another integration by parts)
J(xn+1)
f<P(xn +1)
-ie /
Jp(Xn)
v(xn+l)
-(Ief(g -*)'(0e
-I?
r-p(x„ + l)
£ = p(xn) Jp(xn)
This yields
\J(xn+i)\ < cmm(e0hn, £1h2n, £2hn)
which is equivalent to (43).
□
To get a quadrature for 11 from (34) we apply the scalar quadrature given by (40), (41) to each (off-diagonal) component of I1. Let us denote the matrix valued phase function of E£ by
$ = diag(01,.. .,0d), i.e.
X
Ee{x) = e^x) with $(x) = / L(t)dt.
Hence it holds (see (36))
pij(x) = fa(x) — fa(x) = {D^(x))ij and p'ij(x) = K(x) — Xj(x) = (Dl{x))ij.
We denote with P = P(x) e Cdxd the matrix valued analogue of the interpolation polynomial p ^om (38). With the notation of Sec. 3 it reads:
P (x) = Dl(x) q(C1 + C1 ®Dhx)) .
From (40) we obtain
Cl = (Sn+1 &D-n+i -Sn &D-J © ((D^+1 -D*n +Id)®- Id)
Cl
Sn ©DT - Cl ©D* .
(48)
(49)
(50)
Here denotes the matrix with componentwise reciprocal elements. Since (48)-(50) is just a compact form to write the (componentwise) interpolation function i we get the quadrature
t—xn
(51)
Componentwise application of Lemma 4.1 yields
x
Corollary 4.2. Let L,S E C2([a,b], Cdxd) and let Assumption 2 hold. Then there exists a constant c > 0 independent of e and n, such that for all x E [xn, xn+\] C [a, b] it holds
||/i(x) - Qi(x)|| < ch2n min(e,hn) , (52)
||/\xn+i) - QWi)|| < chn min(h2n-kek). (53)
k=0,2
It remains to find a discretization for the second oscillatory integral 12(xn+i). As a first step we replace Ii(x) in (35) by Qi (x). This yields an error of at most 0(eh2n) in the integrand. Hence, there exists a constant c > 0 independent of e (and n), such that
PXn + l
12(xn+i) - M(x)Qi(x) dx
< C£ h
n
In order to make the following formula more readable we define the matrix function
Pl(x) := Ci + Ci © DHx) - teCi . (54)
Hence it holds
X
t=Xn
which yields
rXn + l i'Xn + 1
M{x)Q\x)dx = ze M{x)e--^l-x]P\x)e^l-x] dx
Cn J %n
fXn + 1
-te / M(x)e-i^x")P\xn)ei^x") dx
f'Xn + 1
£
ie / e-i^x\S(x)Pl(x)ei^x) dx
JXn
For Il(xn+i) we already computed an approximation. Since the remaining integral is of the same type as Ii (replace S by SPi), we can use the same quadrature for the off diagonal elements and shall get the same error estimates (with a different constant c, however). Thus we set
S2 := SPi , (55)
C20
© D-n+1 - Sn © D-J © ((D^n+1 - D$n + Id)®- Id) , (56)
S2n © D-n - Ci © D$n (57)
and the quadrature for the off diagonal elements reads
Qla(xn+i) := + (58)
But, in general, the diagonal of S2 is not zero. Hence we also have to find an appropriate quadrature for
rXn + 1 . . i-Xn + 1
/ diag (e-WS^eW) (ir = / diag (S2(x))dx.
JXn JXn
Since this integral is independent of e we can use a standard (polynomial based) quadrature. Due to the desired order with respect to hn we choose the trapezoid rule, which yields
fXn + 1
' diag (S2(x)) dx - Q%ag(xn+i) < h ,
3
n
n
< cehn. (61)
with a constant c > 0 independent of e and n and
hn
Qdiag(^n+i) := y(diag(S,2+1) + diag(S'2)) . (59)
Finally we define
Q2(Xn+i) := Qlag (Xn+i)+ Qlff (Xn+i) (60)
and hence get the quadrature error estimate
I2{xn+l) - ie[Q2(xn+1) - Q\xn+l)e-^^P\xn)e^] Thus we have proven
Proposition 4.3 (Local error). Let S,L E C2([a,b], Cdxd) and let Assumption 2 hold. Further let 0 < n < N - 1 and let Pl(xn), Ql(xn+\), Q2(xn+i) (and the related quantities) be given by (49)-(60). We define
An := Q1(Xn+i) J
Bn := ie[Q2{xn+l)-Q\xn+l)e~^P\xn)e^^].
Than there exists a constant c > 0 independent of e and n, such that
11z(xn+i) - (id+eAn + e2Bn)z(xn)^ < cehn min h2~kek .
k=0,2 n
Remark 4.4. Let the assumptions of Proposition 4.3 hold. Then it follows from (53) in Corollary 4.2:
\\^J = HQ^Xn+i)! < p^Xn+i)! + chn min(h2n-kek).
k=0,2
Due to (39) every component of Il(xn+i) is at most of order O(e). Further \\Ii(^n+i)\\ is of order 0(hn) which yields
WAnW < cmin(hn, s).
It is easy to see, that ci, c0 from (40) are bounded as hn ^ 0. Hence the same holds for Cj, C^ and consequently for P1, and finally for S2 = SP1. Thus we can use the same argument as before and deduce from (61)
HBnH < cmin(hn, e).
This yields the estimate
||M+e^„ + e2Bn|| < 1 + cehn ,
which guarantees stability and hence convergence of the one step scheme (OSS) from Proposition 4.3.
4.1. Pseudo-code of the scheme
Here we shall write the numerical scheme from Sec. 4 in a pseudo-code. Let a = x0 < xi < ••• < xN = b be our grid and let zn denote the approximation of z at the grid point xn. We start with z° = z0. Assume we have already computed the quantities4 Sn,Ln, D-n, D<s>n, E£,n, and zn.
(1) compute Sn+i, Ln+i, §n+i, and E£,n+i
(2) compute D~-n+x and D^n+1 by (26)
4Here the lower index n denotes the exact quantity evaluated at the grid point xn.
(3) compute C1 and C1 by (50) and (49), i. e.
C1 = (Sn+1 QD-n+i — Sn ©D—) © ((D„n+i — D*n + Id)® — Id) , Cq1 = Sn ©D- — C1 © D$n .
(4) compute P;1 and P1+1 by (54):
P1 = CQ + C\ ©D^n — ieC\ , P1+1 = CQ + C\ ©D„n+i — zeC1.
(5) set Sn = SnP,n, S2+1 = Sn+1Pln+l and compute (see (56), (57))
C2 = (S2n+1 ©D-+ — Sn © D—) © ((D^+i — D$n +Id)® — Id) , CQ2 = S2n © D— — C2 ©D*n
(6) compute (see (51))
Q1 (xn+1) := \_E*£,n+1Pn+1E£,n+1 — E*£,nPnE£,A ■
(7) compute (see (58))
Qln(xn+l) = te e"**' (C02 + C2 © D*, - .
(8) compute (see (59))
hn
Qdiagfcn+i) = yidiagi^n+i) +diag(S'2)) •
(9) set Q (xn+1) = Qlmg (xn+1) + QOff (xn+1) (see (60))
(10) compute An, Bn as defined in Proposition 4.3:
An = Q1 (xn+1),
Bn = ie[Q2{xn+l)-Q\xn+l)e-^Plne^].
(11) compute approximation of z at xn+1:
zn+1 = (Id+eAn + e2Bn) zn
(12) update quantities for the next interval, i. e.
Sn = Sn+1 , Ln = Ln+1 , D-n = D-n+1 , $n = $n+1 , D$n = D$n+1 , E£,n = E£,n+1 .
5. Numerical results
In this section we illustrate the convergence behavior (stated in Proposition 4.3) of our numerical approximation to the solution z of the IVP (32). The results are derived with the scheme from Sec. 4. The procedure how to approximate T£,R£,S,$ is not discussed in this work. For a detailed discussion of the algorithm to compute these variables we refer to [5]. Anyhow, the numerical integration of the phase $ usually incures an additional error for the original, oscillatory function u from (12) or (22). This situation is the same also for scalar ODEs (cf. Th. 3.1 in [2]).
We shall compare our one step scheme (OSS) to the Adiabatic Midpoint Rule (AMPR) from [18]. That integrator is a space-symmetric two-step scheme, which yields a convergence error of order O(e0h2) for the function r] defined in Lemma 2.3. If we would want to have the same error behavior for the original function u, we also have to impose the step size restriction h < \fe (when using the Simpson rule to approximate the matrix valued phase $). Using a higher order quadrature rule for $, would weaken this restriction on h.
Let us choose a family of equidistant grids. Let g E N, and define for n = 0,...,2a =: Na the grid points
x9n := a + nha with ha
b — a
For integers gi < g2, the grid corresponding to gi is a (coarser) subgrid of the g2-grid. Hence no interpolation is needed when comparing solutions on two grids. To generate error plots we fix a finite number of indices, e. g. g = 2,..., 14, and use the numerical solution on the finest grid as reference solution. To illustrate the convergence behavior of the OSS (w.r.t. the step size h, and in dependence of e) we shall give the relative Li-error.
Figures 1-3 show the relative Li-error of z for the example, which is already used in [18] to illustrate the performance of the AMPR (for the details of this example see Sec. 5.1). The graphs in Figure 4 are the relative Li-errors of the variable z, computed with the Kane model of Sec. 2.1. We used the following data: a = 0, b = p(x) = 1, E = 2, V(x) = 10x(|-x), Eg(x) = |sin2(27tx) +
In Figure 1 we plot the theoretical error prediction of Proposition 4.3, with a fitted leading constant. This behavior is reflected quite well in Figure 2, where we used almost exact values for the coefficients appearing in the IVP (32). I. e. we use the approximation of S, L, $ from the finest grid also for the coarser ones. For the simulation of Figure 3 the coefficients S, L, $ were approximated -as it will be done in practice- on the same grid that was used for the solution of the IVP. In the Figures 2, 3 we also observe the error threshold at about 10-i4, resulting from the Matlab computations in double precision.
The numerical experiments confirm the theoretical results. We observe the 0(e°h2) convergence behavior for the AMPR as discussed in [18]. So, the error of that scheme (for the variable ^ from Lemma 2.3) is uniform in e, but it does not decrease as e — 0. However, our OSS shows an even better error behavior than predicted in Proposition 4.3. While for large step sizes h the graphs of the z-error behave like 0(e3h°) (which coincides with the theoretical estimate), they seem to turn to an 0(e2h2) behavior, if h gets small enough (see Fig. 3, 4). This is a "better" convergence property than the predicted 0(elh2) behavior from Proposition 4.3. This behavior was also described in Sec. 3.3 of [2], and it is due to cancellation effects in successive integration steps. The Figures 2-4 also illustrate the asymptotic correctness of our OSS as e — 0, even for rather large values of h.
Both methods, the OSS and the AMPR are subject to the fact that (in general) the transformation back to the original variable u introduces an error of the order 0(e-i). This is due to the multiplication of z (and r]) with the highly oscillatory matrix Ee{x) = exp(^$(x)). Since $ is approximated with the Simpson rule (which yields an error of 0(h4) for $) we get an transformation error of 0(e-ih4). This explains the step size restriction mentioned in the beginning of this section. But if the matrix valued phase function $ is exactly known5, the error behavior of z, ^ carries over to u. In this situation our OSS yields much better results for u than the AMPR - with approximately the same numerical effort.
5.1. The example from [18]
The following example is not related to two-band Schrodinger models. But since it has been studied in the literature [18], it serves as a convenient test and comparison for our method. We
5E. g., piecewise linear functions V, Eg,p in the Kane model lead to an exactly integrable phase.
Fig. 1. Plot of the functions 20min(5e2h2, e3) (solid lines) and 8h2 (dashed lines) for different values of .
relative L*-error in n (AMPR) and z (OSS)
Fig. 2. Relative L1 -error of the OSS for z (solid lines) and the AMPR [18] for r] (dashed lines) for different values of e. "Exact" evaluation of S via interpolation is used.
Fig. 3. Relative Li-error of the OSS for z (solid lines) and the "adiabatic midpoint rule" from [18] for ^ (dashed lines) for different values of e. The function S is numerically approximated as discussed in [5].
Fig. 4. Relative Li-error of the OSS for z related to the Kane model of Sec. 2.1."Exact" evaluation of S via interpolation is used.
shall solve the IVP
(x) + A(x)fax) = 0 , x G I:=[-1,1] , ^(a) = ^o G C2 , fa (a) = fa G C2 ,
with
A(X := (X + 3 2x + ^ '
5 2x + 3
and 8 > 0 some fixed parameter. The diagonalization A = U*AU is given by
A(r) f ¡x + :i + IV^TW o_\2
w V 0 ¡x + 3-^Vx^TWJ '
U{x) = ( COSgXl with ^) = ? + -axctanf4V
v ; ^ — sm^(x) cos^(x) J ; 4 2 \28)
Since A is positive definite on [—1,1], we transform (as for the k ■ p-model in Sec. 2.2)
V\ := A^tfj, v2 := efa, which yields the first order IVP for v(x) = (v]_(x), v2(x))T E C4:
v = - ,i ) v, (62)
e V -A2 0 J V 0 0 J ' v 7
„<«) = (^).
The first matrix of (62) (which we denote by L) has the decomposition L = iQ*LQ with6 Q(x) = ^ ( j ] ) ® U(x), L(x) = ( J ) ® Afcr)* .
1 i J ^ 0 —1
Thus, the equivalent first order IVP for u = Q reads
i L 2
u'ix) = ( ^ ix) ) u(x) + B(x) u(x), (63)
e \ 0 — A 2 (x) J
«(.) (64)
with
B = (J + [^(UAl'A-kr).
In the numerical example we use 8 = 1. If 8 gets smaller, the eigenvalues of A(x = 0) approach each other. This situation is a so called avoided eigenvalue crossing and needs a special numerical treatment. Both schemes, OSS and AMPR, yield poor results for small 8.
6Here <g> denotes the Kronecker (or tensor) product for matrices as defined in [7].
6. Appendix: The diagonalization matrix Q for the two-band k ■ p-model
For the two-band k ■ p-model we can explicitly compute the transformation matrix Q and the eigenvalues of L. The matrix L is given by
/0 o Je[ o \
L
V
0
— \[E~\ 0
0 0
0 0
0
where we set E1 = E — V and E2 = E — V + Eg. The characteristic polynomial x reads
x(A) = A4 + (p2 + Ei + E2)A2 + E1E2 .
Hence, the four (imaginary) eigenvalues of L are
A = ±^r\jp2 + El + E2± v^2 + E1 + E2f - 4 EXE2 .
An eigenvector corresponding to the eigenvalue A is
i\[E~\pE2 \[E~2
v~x =
(
E2
A2 (A2 + p2 + E1) A A(A2 + p2 + E1)
7: ■')
T
Since A is a root of x we get
A2 (A2 + p2 + E1) = —E2(A2 + E1)
which yields
v~x =
(
ip\f~E\ \f~E2 ip\
A2 + Ei A A2 + E1
■1
T
Let iA1iA4 be the four pairwise different eigenvalues of L. Hence it is
L = iQ*LQ
with L = diag(A1,..., A4) and
Q
(
vi
)
with corresponding eigenvectors v1,...,v4.
Acknowledgment. The authors were supported by the FWF (Wissenschaftskolleg "Differentialgleichungen" and project I395-N16).
References
[1] Arnold A. Mathematical properties of quantum evolution equations // In: Quantum transport — modelling, analysis and asymptotics, Lecture Notes in Mathematics, 1946, Allaire G., Arnold A., Degond P., Hou Th.Y. — Berlin: Springer, 2008.
[2] Arnold A., Ben Abdallah N., Negulescu C. WKB-based schemes for the oscillatory 1D Schrodinger equation in the semi-classical limit // SIAM J. Numer. Anal. — 2011. — V. 49. — P. 1436-1460.
[3] Ben Abdallah N., Kefi-Ferhane J. Mathematical analysis of the two-band Schrodinger model // Math. Meth. Appl. Sci. — 2008. — V. 31. —P. 1131-1151.
[4] Ben Abdallah N., Pinaud O. Multiscale simulation of transport in an open quantum system: Resonances and WKB interpolation // J. Comput. Phys. — 2006. — V. 213. — P. 288-310.
[5] Geier J. Efficient integrators for linear highly oscillatory ODEs based on asymptotic expansions // PhD thesis, TU Wien, 2011.
[6] Hairer E., Lubich C., Wanner G. Geometric numerical integration: Structure-preserving algorithms for ordinary differential equations, 2nd Ed. — Berlin, Heidelberg: Springer-Verlag, 2006.
[7] Horn R.A., Johnson C.R. Topics in matrix analysis. — Cambridge: Cambridge Univ. Press, 1991.
[8] Ihlenburg F., Babuska I. Finite element solution of the Helmholtz equation with high wave number. II. The h-p version of the FEM // SIAM J. Numer. Anal. — 1997. — V. 34. — P. 315-358.
[9] Iserles A. On the global error of discretization methods for highly-oscillatory ordinary differential equations // BIT. — 2002. — V. 42. — P. 561-599.
[10] Iserles A., Noersett S.P., Olver S. Highly oscillatory quadrature: The story so far // In: Bermudez de Castro A. (ed.), Proceeding of ENuMath, Santiago de Compostella. — Springer Verlag, 2006. — P. 97-118.
[11] Jahnke T., Lubich C. Numerical integrators for quantum dynamics close to the adiabatic limit // Numer. Math. — 2003.—V. 94.— P. 289-314.
[12] Kefi J. Analyse mathematique et numerique de modeles quantiques pour les semiconducteurs // PhD thesis, Universite Toulouse III, 2003.
[13] Kane E.O. Energy band structure in p-type germanium and silicon // J. Phys. Chem. Solids. — 1956. — V. 1. — P. 82-99.
[14] Kane E.O. Band structure of indium antimonide // J. Phys. Chem. Solids. — 1957. — V. 1. — P. 249-261.
[15] Kane E.O. Energy Band Theory // In: Paul W. (ed.), Handbook on Semiconductors 1. — Amsterdam, New York, Oxford: North-Holland, 1982. — P. 193-217.
[16] Kansy K. Elementare Fehlerdarstellung fur Ableitungen bei der Hermite-Interpolation // Numer. Math. — 1973.—V. 21. —P. 350-354.
[17] Klindworth D. Discrete transparent boundary conditions for multiband effective mass approximations // Diploma thesis, TU Berlin, 2009.
http: / /www.math, tu-berlin. de/^ehrhardt/papers/diplorrLdklindworth.pdf
[18] Lorenz K., Jahnke T., Lubich C. Adiabatic integrators for highly oscillatory second-order linear differential equations with time-varying eigendecomposition // BIT. — 2005. — V. 45. — P. 91-115.
[19] Luttinger J.M., Kohn W. Motion of electrons and holes in perturbed periodic fields // Physical Review. — 1955.—V. 94.— P. 869-883.
[20] Magno R., Bracker A.S., Bennett B.R., Nosho B.Z., Whitman L.J. Barrier roughness effects in resonant interband tunnel diodes // J. Appl. Phys. — 2001. — V. 90. — P. 6177-6181.
[21] Olver S. Moment-free numerical integration of highly oscillatory functions // IMA J. Numer. Anal. — 2006. — V. 26. — P. 213-227.
[22] Stoer J., Bulirsch R. Introduction to numerical analysis. — New York: Springer-Verlag, 1993.
