9 9 9 1 X-ray Diffraction by Time-Dependent Deformed Crystals: n a Theoretical Model and Numerical Analysis J 0 2 Svetlana Sytova ] h Abstract p - The objective of this article is to study the behavior of electromagnetic field under p X-ray diffraction by time-dependent deformed crystals. Derived system of differential m equations lookslike the Takagiequationsin the caseof non-stationarycrystals. This is a o systemofmultidimensionalfirst-orderhyperbolicequationswithcomplextime-dependent c coefficients. Efficient difference schemes based on the multicomponent modification of . s the alternating direction method are proposed. The stability and convergence of devised c i schemesareproved. Numericalresultsareshownforthecaseofanidealcrystal,acrystal s heated uniformly according to a linear law and a time-varying bent crystal. Detailed y h numerical studies indicate the importance of consideration even small crystal changes. p [ MS Classification: 65M06 (Primary), 78A45 (Secondary) 1 Keywords: X-ray Diffraction, Takagi Equations, PDEs, Hyperbolic Systems, Finite v Differences, Scientific Computing 3 3 0 Contents 1 0 9 9 1. Introduction 2 / s c 2. Physical and mathematical model of X-ray dynamical diffraction by time- i s dependent deformed crystals 2 y h 3. Numerical analysis 6 p 3.1. Difference schemes for solving hyperbolic system in two space dimensions 6 : v i 3.2. Stability and convergence of difference schemes 7 X 4. Results of numerical experiments 10 r a 4.1. Diffraction by ideal crystal 10 4.2. Diffraction by time-dependent heated crystal 13 4.3. Diffraction by time-dependent bent crystal 17 5. Summary 18 6. Acknowledgements 19 7. References 19 1 1. Introduction Mathematical modeling of X-ray diffraction by time-dependent deformed crystals refers to physical problems of intensive beams passing through crystals. So, a relativistic electron beam passes through the crystal target and leads to its heating and deformation. Thesystem of differential equations describing X-ray dynamical diffraction by non-stationary deformed crystals was obtained in [1]. This system looks like the Takagi equations [2]– [3] in the case of non-stationary crystals. Up to now in many ref. (e.g. [4]– [7]) the theory of X-ray dynamical diffraction by stationary crystals for differentdeformations was developed. Proper systemsofdifferentialequationsarestationaryhyperbolicsystemsfortwoindependentspatial variables. In [4]– [7] the solutions of these systems were obtained analytically for some cases ofdeformations. [8]iddevotedtonumericalcalculationofpropagationofX-raysinstationary perfect crystals and in a crystal submitted to a thermal gradient. In [9] the theory of time- dependent X-ray diffraction by ideal crystals was developed on the basis of Green-function formalism for some suppositions. The exact analytical solution of the system being studied in this work is difficult if not impossible to obtain. That is why we propose difference schemes for numerical solution. To solve multidimensional hyperbolic systems it is conventional to use different componentwise splitting methods, locally one-dimensional method, alternating direction method and others. They all have one common advantage, since they allow to reduce the solving of complicated problemtosolving ofasystem of simplerones. Butsometimes they donotgive sufficientpre- cision ofapproximation solutionunderratherwidegridspacings andlow solution smoothness because the disbalancement of discrete nature causes the violation of discrete analogues of conservationlaws. Thealternatingdirectionmethodisefficientwhensolvingtwo-dimensional parabolic equations. We use the multicomponent modification of the alternating direction method [10]whichis devoid of suchimperfections. Thismethodprovidesacomplete approx- imation. Itcanbeappliedformulticomponentdecomposition, doesnotrequiretheoperator’s commutability. It can be used in solving both stationary and non-stationary problems. Difference schemes presented allow the peculiarities of initial system solution behavior. The problem of stability and convergence of proposed difference schemes are considered. We present results of numerical experiments carried out. We compare efficiency of suggested schemes in the case of diffraction by ideal stationary crystal. Tests and results of numerical experiments are demonstrated in the case of heated crystal. In our experiments it is assumed that thecrystal was heated uniformlyaccording toa linear law. Thesource of crystal heating was not specified. It may be the electron beam passing through the crystal. In [1] we have given the formulae which allow to determine the crystal temperature under electron beam heating. We present also results of numerical modeling of X-ray diffraction by a time-varying bent crystal. The source of crystal bending is not discussed too. 2. Theoretical Model of X-ray Dynamical Diffraction by Time-Dependent Deformed Crystals We will use the physical notation [11]. Let a monocrystal plane be affected by some time- varying field of forces, which cause the crystal to be deformed. At the same time let a plane electromagnetic wave with frequency ω and wave vector k be incident on this monocrystal plane. We consider two different diffraction geometry which are depicted in Figure.1. In the caseof Bragg geometry thediffracted wave leaves thecrystalthroughthesameplanethatthe direct wave comes in. In Laue case the diffracted wave leaves the crystal through the back plane of the crystal. The electromagnetic field inside the crystal in two-wave approximation 2 Figure 1: Diffraction geometry: a) Bragg case, b) Laue case. is written in the form: D(r,t) = D(r,t)exp(i(kr−ωt))+D (r,t)exp(i((k+τ)r−ωt)), τ where D and D are the amplitudes of electromagnetic induction of direct and diffracted τ waves, respectively, and τ is the reciprocal lattice vector. Let us examine a weakly distorted region in the crystal, where for the deformation vector u(r,t) the following inequalities are correct: ∂u 1∂u ≪ 1, ≪ 1, (cid:12)∂r(cid:12) (cid:12)c ∂t(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) where c is the velocity of light. (cid:12) (cid:12) (cid:12) (cid:12) We can write :τ = τ(1−u). Here τ (r,t) is the reciprocal lattice vector in deformed d d crystal. u(r,t) is the crystal deformation tensor. u = 1/2(∂u /∂x +∂u /∂x ). Let us call ij i j j i b considered system of coordinates S. b Toobtain anexpansioninseriesofthereciprocallattice vector letuspasstoanewsystem of coordinates S′: r = r−u(r,t). Here in each fixed instant of time the Bravais lattice of d deformed crystal is coincident with one of undistorted crystal in the system S. So, in the ′ system S the crystal structure is periodic. And it is disturbed in the system S. ′ Now in S for electric susceptibility ǫ we can write: ǫ(r ;ω) = ǫ(τ ;ω)exp(iτ r ), d d d d Xτd ǫ(r−u;ω) = ǫ(τ ;ω)exp(iτ (r−u)). d d Xτd Or, finally restoring in the system S, we obtain: ǫ(r;ω) = ǫ((1−u)τ;ω)exp(iτ r), d Xτd b 3 where ǫ(0;ω) = 1+g0, ǫ((1−u)τ;ω) = 1+gτ(r,t), ǫ(−(1−u)τ;ω) = 1+g−τ(r,t). LetusassumethattheamplitudesD andD arechanging sufficiently slowly inthespace τ b b and time: 1∂D 1∂D ≪ |D|, τ ≪ |D |, i= 1, 2, 3. τ (cid:12)k ∂x (cid:12) (cid:12)k ∂x (cid:12) (cid:12) i(cid:12) (cid:12) i (cid:12) (cid:12)(cid:12) (cid:12)(cid:12)1 ∂D (cid:12)(cid:12) (cid:12)(cid:12)1 ∂D ≪ |D|, τ ≪ |D |. τ (cid:12)ω ∂t (cid:12) (cid:12)ω ∂t (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) ThenfromMaxwell’s equatio(cid:12)ns the(cid:12)following sy(cid:12)stem of(cid:12)differential equations was derived [1]: 2i∂D 2i ω ∂t + k2kgradD+χ0D+χτDτ = 0, 2i∂D 2i ω ∂tτ + k2kτgradDτ +χ−τD+(χ0−α(r,t)−s(r,t))Dτ = 0, (2.1) where 2k grad(τu) (τ2+2kτ) 2 ∂u ω α(r,t)= α0− τ k2 , α0 = k2 ; s(r,t)= ω(cid:16)τ ∂t (cid:17); k = c. χ0, χ±τ are the zero and ±τ Fourier components of the crystal electric susceptibility. The difference between our system (2.1) and the Takagi equations [3] is in the term α, which depends on time now, and in the appearance of the term s. Let us rewrite the system (2.1) in the generalized form having picked out vectors of σ- polarization from amplitudes of electromagnetic induction D and D and having specified τ three independent variables t, z, x. The spatial variable y is a parameter. One can write a full three-dimensional system. ∂D ∂D ∂D +A11 +A12 +Q11D+Q12Dτ = 0, ∂t ∂z ∂x ∂D ∂D ∂D τ τ τ +A21 +A22 +Q21D+Q22Dτ = 0, (2.2) ∂t ∂z ∂x where ck ck ck ck z x τz τx A11 = , A12 = , A21 = , A22 = ; (2.3) k k k k Q11 = −0.5iωχ0, Q12 = −0.5iωχτ, Q21 = −0.5iωχ−τ, Q22 = −0.5iω(χ0 −α(z,x,t)−s(z,x,t)). (2.4) Initial and boundary conditions are written in the domain G ={(z,x,t),0 ≤ z ≤ L ,0 ≤ z x≤ L ,0 ≤ t ≤ T}. In the Bragg case the boundary conditions are written as follows: x D(0,x,t) = D0, D (L ,x,t) = 0, 0 ≤ x≤ L , t > 0. (2.5) τ z x In Laue geometry, where the diffracted wave leaves the crystal through the crystal back plane, the boundary condition for amplitude D should be written at z = 0. τ 4 As is known [11], the exact solution of the stationary X-ray diffraction problem has the following form: D = c1exp(ikδ1z)+c2exp(ikδ2z), Dτ = c1s1exp(ikδ1z)+c2s2exp(ikδ2z), (2.6) where δ1,δ2 are the solutions of the dispersion equation: (2δγ0 −χ0)(2δγ1 −α0−χ0)−χτχ−τ = 0; 2δiγ0−χ0 s = , i= 1,2; i χ τ γ0 and γ1 are the cosines of the angles between k and kτ, respectively, and the z axis; −D0s2e2(Lz) D0s1e1(Lz) c1 = , c2 = . s1e1(Lz)−s2e2(Lz) s1e1(Lz)−s2e2(Lz) Here and below the following designations are used: e1(z) = exp(ikδ1z), e2(z) = exp(ikδ2z). Let us impose initial conditions corresponding to the exact solution of the stationary X-ray diffraction problem in an ideal crystal: D(z,x,0) = c1e1(z)+c2e2(z), Dτ(z,x,0) = c1s1e1(z)+c2s2e2(z), 0 ≤z ≤ Lz, 0 ≤ x≤ Lx. In the X-ray range the amplitudes (2.6) oscillate with sufficiently high frequency. For large thickness of crystal it is complicated to obtain good numerical solutions of system (2.2) with coefficients (2.4). So, let us find solution of (2.2) for functions D(z,x,t) and D (z,x,t) τ which vary more slowly than e1(z) or e2(z): D(z,x,t) = D(z,x,t)(c1e1(z)+c2e2(z)), Dτ(z,x,t) = Dτ(z,x,t)(c1s1e1(z)+c2s2e2(z)). (2.7) Then the coefficients (2.4) have to be presented by the formulae: χ0 −2kz/k(c1δ1e1(z)+c2δ2e2(z)) Q11 = −0.5iω , c1e1(z)+c2e2(z) c1s1e1(z)+c2s2e2(z) Q12 = −0.5iωχτ , c1e1(z)+c2e2(z) c1e1(z)+c2e2(z) Q21 = −0.5iωχ−τ , c1s1e1(z)+c2s2e2(z) χ0 −α(z,x,t)−s(z,x,t)−2kτz/k(c1s1δ1e1(z)+c2s2δ2e2(z)) Q22 = −0.5iω . (2.8) c1s1e1(z)+c2s2e2(z) The boundary conditions (2.5) take the form: D(0,x,t) = 1, D (L ,x,t) = 1, 0 ≤ x≤ L , t > 0. (2.9) τ z x For this case the initial conditions have to be equal to 1 too. 5 3. Numerical Analysis The original differential problem is a system of multidimensional first-order differential equa- tions of hyperbolic type with complex-valued time-dependent coefficients. The numerical schemes employed in this work are based on the multicomponent modification of the alter- nating direction method. This method was originally developed in [10]. It turned out to be an effective way for efficient implementations of difference schemes. This method is econom- ical and unconditionally stable without stabilizing corrections for any dimension problems of mathematical physics. It does not require spatial operator’s commutability for the validity of the stability conditions. This method is efficient for operation with complex arithmetic. Its main idea is in the reduction of the initial problem to consecutive or parallel solution of weakly held subproblems in subregions with simpler structure. That is why it allows us to perform computations on parallel computers. The main feature of this method implies that the grids for different directions can be chosen independently and for different components of approximate solution one can use proper methods. Let introduce the Hilbert space of complex vector functions: H = L2(G). In this space the inner product and the norm are applied in the usual way: (u,v) = u(x)v(x)dx, kuk = (u,u)1/2. Z G In H our system (2.2) is hyperbolic. 3.1. Efficient Schemes for Solving Hyperbolic System in Two Space Dimensions We use the following notation [12]: yx = (yi+1−yi)/hx — right difference derivative, yx = (yi−yi−1)/hx — left one, yi = y(xi); yt = (y−y)/ht, y = y(tk+1), y = y(tk). Let us replace the domaibn G of the contbinuous change of variables by the grid domain Gzxt = {(zi,xj,tk);zi = ihz,i = 0,1,...,N1,N1 = [Lz/hz],xj = jhx,j = 0,1,...,N2, N2 = [Lx/hx],tk = kht,k = 0,1,...,N3,N3 = [T/ht]}. The following system of difference equations approximates on G the system (2.2) with zxt coefficients (2.3)–(2.4): ∗ ∗ 1 1 2 1 1 Dt +A11Dz +A12Dx+Q11D +Q12Dτ = 0, ∗ ∗ Dτ1t +A21Dτ1zb+A22Dτ2x +Q21Db1 +(Q22Dbτ1) = 0, (3.10) b b d ∗ ∗ 2 1 2 1 1 Dt +A11Dz +A12Dx+Q11D +Q12Dτ = 0, ∗ ∗ Dτ2t +A21Dτ1zb+A22Dτ2xb+Q21Db1 +(Q22Dbτ1) = 0, (3.11) b ∗ b b d ∗ where D = 0.5(Di +Di−1), Q22 = 0.5(Q22(zi−1,xj,tk+1)+Q22(zi,xj,tk+1)) for the first ∗ ∗ equations of (3.10) and (3.11) dand D = 0.5(Di + Di+1), Q22 = 0.5(Q22(zi,xj,tk+1) + Q22(zi+1,xj,tk+1)) for the last ones. d 6 For coefficients (2.3), (2.8) the system (2.2) can be approximated by the system of differ- ence equations of the following form: 1 1 2 1 1 Dt +A11Dz +A12Dx +Q11D +Q12Dτ = 0, Dτ1t+A21Dτ1bz +A22Dτ2x+Q21Db1+Q22Db1τ = 0; b b d 2 1 2 1 1 Dt +A11Dz +A12Dx +Q11D +Q12Dτ = 0, Dτ2t+A21Dτ1bz +A22Dbτ2x+Q21Db1+Q22Db1τ = 0. (3.12) b b b d In cited schemes the directions of difference derivatives with respect to x (left or right) are selected in dependence of waves directions. D1, D2, D1 and D2 are two components τ τ of approximate solutions for D and D , respectively. One can choose any of these two τ components or its half-sum as a solution of (2.2). The boundary and initial conditions are approximated in the accurate form. In Laue case where the diffracted wave moves on a positive direction of the z axis, for D τ we should write left difference derivatives with respect to z. The schemes (3.10)–(3.11) and (3.12) are completely consistent. The consistency clearly follows from the manner in which these schemes were constructed. On sufficiently smooth solutions they are of the first order approximation with respect to time and space. We can give the difference scheme of the second order approximation with respect to z. In this case it should be rewritten (3.10): 1∗ 1 2 1∗ 1∗ Dt +A11Dz +A12Dx+Q11D +Q12Dτ = 0, Dτ1t∗+A21Dτ1zb+A22Dτ2x+Q21Db1∗+(Q22Dbτ1)∗ = 0. (3.13) b b d The scheme for the second component (3.11) is not changed. One can write a scheme of the second order approximation with respect to time. But as has been shown in numerical experiments it does not lead to sensible changes in solution pattern. For the difference schemes presented, the stability relative to initial data and also the convergence of the difference problem solution to the solution of differential problem (2.2) can be proved. This follows from the properties of the multicomponent modification of the alternating direction method [10]. Let us prove the corresponding Theorems. 3.2. Stability and Convergence of Difference Schemes We use the energy inequalities method [12]. Let rewrite the system (3.10)–(3.11) in the form: D1t +Λ1(D1)+Λ2(D2) =0, (3.14) D2t +bΛ1(D1)+Λ2(D2) = 0, (3.15) b b where ∗ ∗ D = D , Λ1(D)= A11Dz +Q11D∗ +Q12Dτ ∗ , Λ2(D) = A12Dx . (cid:18)Dτ (cid:19) (cid:18)A21Dτz +Q21D +(Q22Dτ) (cid:19) (cid:18)A22Dτx(cid:19) Let us introduce the following notation: ′ ′′ ′ ′′ y =Re(y), y = Im(y), y = y −iy . 7 We use the inner products: N−1 (y,v) = hy v , ω i i Xi=1 where ω = {x = ih,i = 0,1,...,N,Nh =L} is a one-dimensional grid; i N1−1N2−1 (Y,V) = (Y,V) = h h Y V . Gzxt z x ij ij Xi=1 jX=1 In addition to this let us introduce the norm: 2 ′ 2 ′′ 2 kYk = (Y,Y), |Y| = kY k +kY k . q Lemma. If y(x0) = 0 then (y,yx)ω ≥ 0, if y(xN)= 0 then (y,yx)ω ≤ 0. Proof. Let us write the following transformation chain: N−1 N−1 2 2 (y,yx)ω = yi(yi−yi−1)= (yi −yi−1−yi−1(yi−yi−1)) = Xi=1 Xi=1 N−1 N−2 2 2 2 2 yN−1−y0 − yi−1(yi−yi−1) = yN−1−y0 − yi(yi+1−yi)= Xi=1 Xi=0 N−2 N−2 N−2 N−2 2 2 2 2 2 2 yN−1−y0 + yi −0.5 (yi+1+yi) +0.5 yi+1+0.5 yi = Xi=0 Xi=0 Xi=0 Xi=0 N−2 N−2 N−2 2 2 2 2 2 0.5yN−1 −0.5y0 + yi + yi+1−0.5 (yi+1+yi) . Xi=0 Xi=0 Xi=0 The first term in the last expression is greater than or equal to 0, the second is equal to 0. Taking into consideration the inequality: 2 2 2 (yi+1+yi) ≤ 2(yi+1+yi), we obtain: (y,y ) ≥ 0. Second Lemma’s inequality is proved similarly. 2 x ω Teorem 3.1. The difference scheme (3.10)–(3.11) is unconditionally stable relative to the initial data. For its solution the following estimates hold: |Di|2 ≤ M |Di(t0)|2+|Λ1(D1(t0))+Λ2(D2(t0))|2 , (3.16) (cid:16) (cid:17) b where M is a bounded positive constant independent of grid spacings, i= 1,2. Proof. Multiply (3.14) by (Λ1(D1))t, (3.15) by (Λ2(D2))t and sum : (D1t)′,(Λ1(D1))′t + (D1t)′′,(Λ1(D1))′t′ + (cid:16) (cid:17) (cid:16) (cid:17) (D2t)′,(Λ2(D2))′t + (D2t)′′,(Λ2(D2))′t′ + (cid:16) (cid:17) (cid:16) (cid:17) (Λ1(D1))′,(Λ1(D1))′t + (Λ1(D1))′′,(Λ1(D1))′t′ + (cid:16) (cid:17) (cid:16) (cid:17) (Λb2(D2))′,(Λ1(D1))′t + (Λb2(D2))′′,(Λ1(D1))′t′ + (cid:16) (cid:17) (cid:16) (cid:17) (Λ1(D1))′,(Λ2(D2))′t + (Λ1(D1))′′,(Λ2(D2))′t′ + (cid:16) (cid:17) (cid:16) (cid:17) (Λ2b(D2))′,(Λ2(D2))′t + (Λ2b(D2))′′,(Λ2(D2))′t′ = 0. (3.17) (cid:16) (cid:17) (cid:16) (cid:17) b b 8 Let us multiply (3.17) by h and take into account the form of time derivatives. Then we t obtain: Φ(D)+0.5 |Λ1(D1)+Λ2(D2)|2+|Λ1(D1)+Λ2(D2)|2 + (cid:16) (cid:17) 0.5 |Λ1(Db1)|2+|Λ2b(D2)|2+|Λ1(D1)|2 +|Λ2(D2)|2 − (cid:16) (cid:17) (bΛ1(D1))′,(Λb1(D1))′ − (Λ2(D2))′,(Λ2(D2))′ − (cid:16) (cid:17) (cid:16) (cid:17) (Λ1(Db1))′′,(Λ1(D1))′′ − (Λ2(bD2))′′,(Λ2(D2))′′ = 0, (3.18) (cid:16) (cid:17) (cid:16) (cid:17) b b where Φ(D) = ht (Dt1)′,(Λ1(D1))′t +ht (Dt1)′′,(Λ1(D1))′t′ + (cid:16) (cid:17) (cid:16) (cid:17) ht (Dt2)′,(Λ2(D2))′t +ht (Dt2)′′,(Λ2(D2))′t′ . (3.19) (cid:16) (cid:17) (cid:16) (cid:17) Let us re-arrange (3.18): Φ(D)+0.5|Λ1(D1)+Λ2(D2)|2−0.5|Λ1(D1)+Λ2(D2)|2 + b 0.5h2t|b(Λ1(D1))t|2+0.5h2t|(Λ2(D2))t|2 = 0. Or, more properly, Φ(D)+0.5|Λ1(D1)+Λ2(D2)|2 ≤ 0.5|Λ1(D1)+Λ2(D2)|2. (3.20) IfweprovethatΦ(D) ≥b0thenthebestimate (3.16) willbeobtained fori= 2. We consider the first term in (3.19) (Dt1)′,(Λ1(D1)′t + (Dt1)′′,(Λ1(D1)′t′ = (cid:16) (cid:17) (cid:16) (cid:17) 1 ′ 1 ′ 1 ′ 1∗ ′′ 1 ′ 1∗ ′′ (Dt),(A11Dz)t − (Dt),Q11(D )t − (Dt),Q12(Dτ )t + (cid:16) (cid:17) (cid:16) (cid:17) (cid:16) (cid:17) 1 ′ 1 ′ 1 ′ 1∗ ′′ 1 ′ 1∗ ′′ (Dτt),(A21Dτz)t − (Dτt),Q21(D )t − (Dτt),(Q22Dτ )t + (cid:16) (cid:17) (cid:16) (cid:17) (cid:16) (cid:17) 1 ′′ 1 ′′ 1 ′′ 1∗ ′ 1 ′′ 1∗ ′ (Dt) ,(A11Dz)t + (Dt) ,Q11(D )t + (Dt) ,Q12(Dτ )t + (cid:16) (cid:17) (cid:16) (cid:17) (cid:16) (cid:17) 1 ′′ 1 ′′ 1 ′′ 1∗ ′ 1 ′′ 1∗ ′ (Dτt) ,(A21Dτz)t + (Dτt) ,Q21(D )t + (Dτt) ,(Q22Dτ )t ≥ 0, (3.21) (cid:16) (cid:17) (cid:16) (cid:17) (cid:16) (cid:17) as since from Lemma we have: 1 ′ 1 ′ 1 ′′ 1 ′′ ′ (Dt),(A11Dz)t ≥ 0, (Dt) ,(A11Dz)t ≥ 0, A11 = (A11) > 0; (cid:16) (cid:17) (cid:16) (cid:17) 1 ′ 1 ′ 1 ′′ 1 ′′ ′ (Dτt),(A21Dτz)t ≥ 0, (Dτt) ,(A21Dτz)t ≥ 0, A21 = (A21) < 0. (cid:16) (cid:17) (cid:16) (cid:17) In (3.21) we took into account that the coefficients Q from (2.4) were pure imaginary, Q12 = Q21. Also we have written out corresponding inner products. In the same fashion as Lemma we obtain for the second term in (3.19): (Dt2)′,(Λ2(D2)′t + (Dt2)′′,(Λ2(D2)′t′ ≥ 0. (3.22) (cid:16) (cid:17) (cid:16) (cid:17) b b So, expressions (3.21) and (3.22) mean that operators Λ1 and Λ2 are positive definite. If we repeatedly apply the expression (3.20) to the right and if we take into account that D2t = −Λ1(D1)−Λ2(D2), (3.23) b b 9 we have the estimate (3.16) for i = 2. Now let us obtain the estimate (3.16) for i= 1. Multiply D1 by D1 and take into account the expression which follows from our designations: b b D1 = D1−ht Λ1(D1)+Λ2(D2) . (cid:16) (cid:17) b b Using ǫ-inequality [12] |(u,v)| ≤ ǫkuk2 +1/(4ǫ)kvk2, (ǫ > 0), we have : |D1|2 = D1−ht Λ1(D1)+Λ2(D2) ,D1−ht Λ1(D1)+Λ2(D2) = (cid:18) (cid:16) (cid:17) (cid:16) (cid:17)(cid:19) b |D1|2+h2t|Λ1b(D1)+Λ2(D2)|2−2ht (D1)b′,(Λ1(D1)+Λ2(D2))′ (cid:16) (cid:17) −2ht (D1)′′,(Λ1(Db1)+Λ2(D2))′′ ≤ M1|D1|2+M2b|Λ1(D1)+Λ2(D2)|2. (cid:16) (cid:17) b b Let us consider the second term in previous inequality: |Λ1(D1)+Λ2(D2)|2 = |Λ1(D1)+Λ2(D2)−ht(Λ2(D2))t|2 = |Λ1(D1)+Λ2(D2)|2 +b h2t|(Λ2(D2))t|2−2htb (Λ1(D1)b+Λ2(D2))′,(Λ2(D2))′t − (cid:16) (cid:17) b b 2ht (Λ1(D1)+Λ2(D2))′′,(Λb2(D2))′t′ b≤ |Λ1(D1)+Λ2(D2)|2. (cid:16) (cid:17) b b b b This was obtained by taking into account (3.23) and (3.22). Finally we have: |D1|2 ≤ M3 |D1|2+|Λ1(D1)+Λ2(D2)|2 . (cid:16) (cid:17) b Repeated application of this inequality to the right yields the estimate (3.16) for i= 1. 2 We denote the discretization error as Zi = Di−D, i= 1,2, where D is the exact solution of initial differential problem. Teorem 3.2. Let the differential problem (2.2)–(2.5) have a unique solution. Then the solution of the difference problem (3.10)–(3.11) converges to the solution of the initial differential problem as h ,h ,h → 0. The discretization error may be written as t z x |Zi| ≤ O(h +h +h ). t z x Proof follows immediately from consistency of the scheme (3.10)–(3.11), Theorem 3.1 and Lax’s Equivalence Theorem [13]. 2 The stability and convergence of schemes (3.12), (3.13)–(3.11) can be proved in an anal- ogous way. 4. Results of Numerical Experiments 4.1. Diffraction by Ideal Crystal TheproblemofstudyingofelectromagneticfieldsunderX-raydiffractioninsidethecrystal target is constituent of the problem of modeling intensive beams passing through crystals, X- rayfreeelectronlaserandothers. Thereforeletanalyzetheoperationofschemes(3.10)–(3.11) and (3.12) in the case of ideal absorbing crystals. Figure 2–5 display results of numerical experiments in the crystal of LiH. This crystal was chosen because of small absorption. The design parameters were the following. The 10