ON A CLASS OF CONSERVATIVE, HIGHLY ACCURATE GALERKIN ...
ON A CLASS OF CONSERVATIVE, HIGHLY ACCURATE ¨ GALERKIN METHODS FOR THE SCHRODINGER EQUATION GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
Abstract. We construct and analyze efficient fully discrete Galerkin type methods that are of high order of accuracy and conservative in the L2 sense for approximating the solution of a form of the linear Schr¨odinger equation with a time-dependent coefficient, found e.g. in underwater acoustics. The time stepping procedures are based on the class of implicit Runge–Kutta methods known as the q−stage Gauss– Legendre schemes. L2 error estimates are proved that are of optimal order in space and of temporal order q + 2. An iterative procedure at each time step for the efficient implementation of the two-stage scheme is proposed and analyzed. ` ` . On construit et analyse des m´ RESUM E ethodes totalement discr`etes du type Galerkin, qui sont L2 −conservatives et d’order arbitraire, pour approcher la solution d’une forme de l’´equation lin´eaire de Schr¨odinger avec un coefficient qui d´epend du temps, trouv´ee par exemple dans l’acoustique sous-marine. La proc´edure de discr´etisation en temps est bas´ee sur la classe des m´ethodes implicites de Runge–Kutta connues comme les sch´emas de Gauss–Legendre `a q pas interm´ediaires. On obtient des estimations dans L2 pour les erreurs, qui sont d’order optimal en espace es d’order q + 2 en temps. On propose et analyse aussi une proc´educe it´erative pour r´esoudre les syst`emes lin´eaires ` a chaque pas de temps pour une application efficace de sch´ema Gauss– Legendre ` a deux pas interm´ediaires.
Dedicated to Professor Dr. G. H¨ ammerlin on the occasion of his 60th birthday, July 31, 1988.
1. Introduction In this paper we shall study conservative numerical methods of high order of accuracy for approximating the solution of the following initial- and boundary-value problem for a partial differential equation of the Schr¨odinger type. Let Ω be a bounded domain in RN with smooth boundary ∂Ω and let 0 < T < ∞ be given. We seek a complex-valued ¯ × [0, T ], satisfying: function u = u(x, t), (x, t) ∈ Ω ut = iL(t)u := i α∆u + β(x, t)u in Ω × [0, T ],
(1.1)
u=0
on ∂Ω × [0, T ], ¯ in Ω,
u(x, 0) = u0 (x)
where α is a given nonzero real number, β is a given smooth real-valued function on ¯ × [0, T ] and u0 is a given smooth complex-valued function on Ω. ¯ We shall assume Ω Work supported by the Institute of Applied and Computational Mathematics of the Research Center of Crete-FORTH. 1
2
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
that the data of (1.1) are smooth and compatible enough to ensure that the problem possesses a unique and smooth enough for our purposes solution; cf. e.g. [15, Chapter 5, Section 12] and [4] for relevant existence, uniqueness and regularity results. This form of the Schr¨odinger equation occurs, for example for N = 1, in underwater acoustics as “parabolic approximation” to the Helmholtz equation, cf. e.g. [20], posed with a variety of types of boundary conditions. For simplicity we consider here the case of homogeneous Dirichlet boundary conditions. In (1.1) the Laplacian ∆ could have as well been replaced by a second-order, symmetric, uniformly positive definite elliptic ¯ with space-dependent coefficients with no complications in the error operator on Ω estimates. We shall discretize (1.1) in space by a Galerkin-finite element type method and in time by a class of implicit Runge–Kutta schemes of arbitrary order, known as the Gauss–Legendre collocation type methods. We shall estimate the error of the fully discrete approximations in L2 and point out efficient ways for implementing the methods. Many finite difference and spectral schemes, usually of second-order temporal accuracy, have been proposed in the literature for problems such as (1.1); see e.g. the survey [13], the collections of papers in [19] and [14] and their references. Galerkin-finite element methods have also been considered. For the linear Schr¨odinger equation with time-independent coefficients early error estimates for semidiscrete approximations may be found in [21], [23]. Semidiscrete and fully discrete Galerkin approximations with Runge–Kutta time stepping have been analyzed in [4] in the general case where the Laplacian in the right-hand side of the p.d.e. is replaced by a second-order elliptic operator with space- and time-dependent coefficients. For alternative approaches based on separating real and imaginary parts, cf. [10], [16]; in this paper we shall discretize (1.1) directly using complex arithmetic. Among the growing literature on Galerkin type methods for nonlinear Schr¨odinger equations, some error estimates are shown in [17], where second-order time stepping procedures coupled with finite difference or Galerkin type space discretizations are analyzed; for computations with such schemes cf. e.g. [18]. We next introduce notation that will be used in the sequel. For integral s ≥ 0 let H s = H s (Ω) denote the usual, complex (Hilbert) Sobolev spaces with corresponding norm k · ks . For f, g ∈ L2 = H 0 , let (f, g) =
Z
f (x)g(x) dx Ω
be their L2 inner product; here the overbar denotes complex conjugation. Let k · k denote the associated L2 norm and | · |∞ be the norm of L∞ = L∞ (Ω). As usual, H01 will consist of the elements of H 1 that vanish on ∂Ω in the sence of trace. We shall discretize (1.1) in space by the standard Galerkin method, as follows. For 0 < h < 1, let Sh be a family of finite-dimensional subspaces of H01 in which approximations to the solution u(·, t) of (1.1) will be sought for given t ∈ [0, T ]. We assume that Sh satisfies the approximation property that there exists an integer r ≥ 2 and a constant c > 0
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
3
independent of h such that for v ∈ H s ∩ H01 (1.2)
inf kv − ϕk + hkv − ϕk1 ≤ chs kvks ,
ϕ∈Sh
1≤s≤r
and the inverse property that for some c > 0 independent of h, kϕk1 ≤ ch−1 kϕk,
(1.3)
∀ϕ ∈ Sh .
As in (1.2), (1.3), in the sequel the symbols c, C, ci etc. will denote generic constants independent of the discretization parameter h and the time step. Such constants may also depend on the solution and the data of (1.1). We define now the semidiscrete approximation of the solution u(t) = u(·, t) of (1.1) in Sh in the customary way as the map uh : [0, T ] → Sh satisfying (with β(t) = β(·, t)): (1.4)
(uht, ϕ) = −iαa(uh , ϕ) + i(β(t)uh , ϕ),
uh (0) = u0h ,
where, for ϕ, χ ∈ H 1 , a(ϕ, χ) =
N X
∀ϕ ∈ Sh , 0 ≤ t ≤ T,
(∂j ϕ, ∂j χ).
j=1
In addition, we shall henceforth assume that u0 ∈ H r ∩ H01 and that u0h is an element of Sh such that (1.5)
ku0 − u0h k ≤ chr ku0 kr ;
for example u0h = P u0, where P : L2 → Sh is the L2 −projection operator onto Sh . If we introduce the linear operators ∆h : Sh → Sh and Bh (t) : Sh → Sh , Lh (t) : Sh → Sh , 0 ≤ t ≤ T, by (1.6) (1.7) (1.8)
(∆h ϕ, χ) = −a(ϕ, χ), (Bh (t)ϕ, χ) = (β(t)ϕ, χ),
∀ϕ, χ ∈ Sh ,
∀ϕ, χ ∈ Sh ,
Lh (t) = α∆h + Bh (t),
i.e.,
Bh (t)ϕ = P [β(·, t)ϕ],
0 ≤ t ≤ T,
and note that, since α and β are real, ∆h , Bh (t) and Lh (t), 0 ≤ t ≤ T, are Hermitian operators on {Sh , (·, ·)}, we may write (1.4) as (1.9)
uht = iLh (t)uh ,
0 ≤ t ≤ T,
uh (0) = u0h .
The equations (1.4) or (1.9) represent an initial-value problem for a system of ordinary differential equations that obviously has a unique solution uh (t) for 0 ≤ t ≤ T ; they will only be used in the sequel in order to motivate the fully discrete approximations. Let us only remark here that it is not hard to show, by comparing e.g. uh to the elliptic projection of u in the standard way, cf. [22], [4] that if u0h satisfies (1.5), then for t ∈ [0, T ] there holds Z t n o 0 r k(u − uh )(t)k ≤ ch ku kr + ku(s)kr + kut (s)kr ds , 0
4
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
i.e. that uh (t) satisfies an optimal-order L2 error estimate if u is smooth enough. We shall discretize (1.9) in time using the well-known class of implicit Runge–Kutta procedures of collocation type known as the q−stage (q ≥ 1 integer) Gauss–Legendre methods, [5], [6], [9]. The methods are defined by constants A = (aij ) ∈ Rq×q , b = (b1 , . . . , bq )T ∈ Rq , τ = (τ1 , . . . , τq )T ∈ Rq that are constructed in the standard way. Specifically, the τi are the —distinct, in (0, 1)— zeros of the shifted Legendre polynomials (d/dx)q xq (1 − x)q , the weights bi are defined as the solution of the q × q linear system of equations represented by (3.1.1) below for 1 ≤ ℓ ≤ q, while the coefficients aij are defined for each i, 1 ≤ i ≤ q, as solutions of the q × q linear system of equations represented by (3.1.2). Let k > 0 be the (constant) time step, let tn = nk, n = 0, 1, . . . , M, where T = Mk, and tn,i = tn + τi k, 1 ≤ i ≤ q. Then, the q−stage Gauss–Legendre methods applied to the system of ordinary differential equations represented by (1.9) yield the following full discretization of (1.1): For 0 ≤ n ≤ M, we seek U n ∈ Sh , approximating u(tn ), and U n,m ∈ Sh , 1 ≤ m ≤ q, that satisfy (a) U 0 = u0h ,
for n = 0, 1, . . . , M − 1 : (1.10)
(b) U
n,m
n
= U + ik
q X
n,j amj Ln,j , h U
j=1
(c) U
n+1
n
= U + ik
q X
1 ≤ m ≤ q,
n,j bj Ln,j , h U
j=1
Lh (tn ), Ln,j h
where Lnh = = Lh (tn,j ). In Section 2 we shall show that, for each n, the system (1.10b) has a unique solution {U n,m }, 1 ≤ m ≤ q, in (Sh )q and therefore that U n+1 is defined uniquely in Sh for 0 ≤ n ≤ M − 1, by (1.10c). We shall also verify that (1.10) is unconditionally stable in L2 and, in particular, conservative, i.e. that it satisfies kU n k = kU 0 k, 0 ≤ n ≤ M; thus it mimicks the behavior of (1.1) and (1.9) for which ku(t)k = ku0k and kuh (t)k = ku0h k hold, respectively, for 0 ≤ t ≤ T. In Section 3 we shall estimate the error of the approximation U n in the L2 norm; specifically, we shall show in Theorem 3.1, which is the main result of the paper, that (1.11) max kU n − u(tn )k ≤ c k min(2q,q+2) + hr . 0≤n≤M
This error bound is of optimal order in space. As far as the temporal rate of convergence is concerned, it is well-known that the q−stage Gauss–Legendre methods have (classical) order 2q when applied to nonstiff ordinary differential equations. Therefore our proof certainly gives optimal-order temporal convergence, resp. 2, 4, for the one- and two-stage, resp., methods that seem to enjoy current practical importance. For q > 2 our result — O(k q+2) in time — shows the effect of “reduction of order due to stiffness”. This result is no worse than analogous estimates proven in the literature for Runge– Kutta full discretizations of initial- and boundary-value problems for p.d.e.’s with timedependent coefficients or nonlinear terms, posed with Dirichlet boundary conditions: In
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
5
his thesis Broc´ehn, [4], considers full discretizations of the Schr¨odinger equation with a general second-order elliptic operator with time-dependent coefficients using some semi-implicit Runge–Kutta methods. (The class of Gauss–Legendre schemes considered here is not semi-implicit with the exception of the one-stage scheme.) For his schemes, using different estimation techniques, he proves error bounds with temporal order of accuracy equal to min(p, q + 1), where, in his notation, q is the order of the quadrature rule associated with the intermediate stages of the Runge–Kutta method and plays an analogous role to the q used here, and p is the classical order (equal to 2q for the Gauss–Legendre methods). For the special elliptic operator of the right-hand side of (1.1) he obtains a temporal rate of q + 2, if p = q − 2, and under the mesh condition that kh−2 remain bounded as k, h → 0. In Theorem 2 of [8], the temporal discretization, by a class of Runge–Kutta schemes (disjoint from the Gauss–Legendre methods but suitable for parabolic problems) of an abstract semilinear parabolic equation is shown to have a rate of convergence exhibiting an analogous limitation to our q + 2 result. In [3] the Gauss–Legendre methods are applied to a nonlinear p.d.e., the generalized Korteweg–de Vries equation, posed in one space dimension with periodic boundary conditions and discretized in space with smooth periodic splines on a uniform mesh. The exact analog to (1.11) is proved then by a different technique from the one at hand, with the details of the space discretization and the periodicity of the exact and discrete solutions playing a crucial role. (After the completion of the original version of the paper at hand, we learnt that Karakashian and McKinney, [12], proved the optimal order 2q for the Korteweg–de Vries equation; their remarkable proof again relies heavily on the periodic boundary conditions.) Finally, in Section 4 we confine attention to the 2−stage Gauss–Legendre method and devise a scheme that avoids solving the 2 dim Sh × 2 dim Sh linear system represented by (1.10b) for q = 2. A suitable decoupling strategy and an iteration scheme enables us to produce stable and optimal-order accurate approximations to U n by solving a number of linear systems of size dim Sh × dim Sh at each time step; these systems will have sparse matrices if Sh is furnished with a finite element basis with elements of small support. Acknowledgement. The authors record their thanks to a referee for pointing out reference [4] to them.
2. Existence and stability of the fully discrete approximations In this section we shall show that for each n the linear system represented by (1.10b) has a unique solution U n,m , 1 ≤ m ≤ q, and that the resulting overall fully discrete scheme (1.10) is stable (conservative) in the L2 norm. For this purpose we shall make use of the following well-known properties of the Gauss–Legendre methods, [6], [9]:
6
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
(2.1) For each q ≥ 1 there exists a diagonal q × q matrix D with positive
diagonal elements, such that the matrix DAD −1 is positive definite on Rq . (See e.g. [9, Theorem 5.5.6, Cor. 5.1.4 and (5.1.23)].)
(2.2)
bi > 0, bi aij + bj aji − bi bj = 0,
1 ≤ i, j ≤ q,
i.e. the Gauss–Legendre methods are conservative in the nonlinear context; cf. e.g. [7], [9, p.117]. We first examine the existence of solutions U n,m of the linear system (1.10b). On the product space (Sh )q we let Un denote the vector (U n,1 , . . . , U n,q )T and Lnh : (Sh )q → n,q T (Sh )q be the diagonal operator given, for V ∈ (Sh )q , by Lnh V = (Ln,1 h V1 , . . . , Lh Vq ) . We write then the equations (1.10a,b) respectively as (2.3a)
Un = U n e + ikALnh Un ,
(2.3b)
U n+1 = U n + ikbT Lnh Un .
In (2.3a,b) and in the sequel, abusing notation a bit to avoid tensor products, for V ∈ P P (Sh )q we let bT V = qi=1 bi Vi , AV ∈ (Sh )q : (AV)i = qj=1 aij Vj , e = (1, 1, . . . , 1)T ∈ Rq and for U ∈ Sh , Ue = (U, U, . . . , U)T ∈ (Sh )q . The existence of Un , solution of (2.3a), will follow from the following general lemma. (In the sequel we let ((·, ·)), resp. ||| · |||, denote the product inner product, resp. norm, in (L2 )q .) Lemma 2.1. Suppose that V = (Vi ) and W = (Wi ) in (S h )q satisfy the equation (2.4)
V = W + ikAF(V),
where F : (Sh )q → (Sh )q is a diagonal mapping such that (F(V ))i = Fi (Vi ), 1 ≤ i ≤ q, where Fi : Sh → Sh are given mappings with the property that Im(Fi (ϕ), ϕ) = 0, for ϕ ∈ Sh , 1 ≤ i ≤ q. Then, there exists a constant c, depending only on the constants of the Gauss–Legendre method, such that (2.5)
|||V||| ≤ c|||W|||.
Proof. Let D = diag (d1 , . . . , dq ), di > 0, be the diagonal matrix mentioned in (2.1). Multiplying (2.4) by D 2 A−1 on the left and taking the (L2 )q inner product with V we obtain (2.6)
((D 2 A−1 V, V)) = ((D 2 A−1 W, V)) + ik((D 2 F(V), V)).
By our hypotheses Re[ik((D 2 F(V), V))] = 0. Hence, taking real parts in (2.6) gives (2.7)
Re((D 2 A−1 V, V)) = Re((D 2 A−1 W, V)).
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
7
Denote DA−1 D −1 =: B = (bij ) ∈ Rq×q . By (2.1) there exists λ > 0 such that Pq Pq 2 q i=1 ξi for every ξ ∈ R . Hence, putting DV = Y we have i,j=1 bij ξi ξj ≥ λ Re((D 2 A−1 V, V)) = Re((DA−1 D −1 Y, Y)) = Re((BY, Y)) ≥ λ ||| Re Y|||2 + ||| Im Y|||2 = λ|||Y|||2 ≥ λ min d2i |||V|||2 . i
(2.5) follows then by (2.7) and the Cauchy–Schwarz inequality.
Given now U n ∈ Sh apply Lemma 2.1 to the linear system (2.3a) where F(V) = n,i Lnh V, Fi (Vi ) = Ln,i h Vi and Im(Lh ϕ, ϕ) = 0 for ϕ ∈ Sh . By (2.5) we see that the homogeneous system (U n = 0 in (2.3a)) has only the trivial solution. Hence, given U n ∈ Sh , (2.3a) has a unique solution Un = (U n,i ) in (Sh )q , which satisfies max kU n,i k ≤ ckU n k,
(2.8)
i
for some constant c that depends only on the Gauss–Legendre method. The stability (conservativeness) of the scheme (1.10) follows from the following result, stated again in slightly more general terms: Lemma 2.2. Suppose, given U ∈ Sh , that V ∈ (Sh )q and Y ∈ Sh satisfy the equations (2.9a)
V = Ue + ikAF(V),
(2.9b)
Y = U + ikbT F(V),
where F is a mapping that satisfies the hypotheses of Lemma 2.1. Then (2.10)
kY k = kUk.
Proof. (2.9b) gives kY k2 = kUk2 + ik(bT F(V), U) − ik(U, bT F(V)) + k 2 (bT F(V), bT F(V)) 2
= kUk − 2k Im
q X
bj (Wj , U) + k
j=1
2
q X
bj bℓ (Wj , Wℓ ),
j,ℓ=1
where Wj = Fj (Vj ). Let W = (W1 , . . . , Wq )T . Replacing U in the right-hand side of the above equation by its expression U = Vj − ik(AW)j that (2.9a) gives, we see, using the properties of F and (2.2) that 2
2
kY k = kUk − 2k 2
= kUk − k
2
2
q X
bi Re(Wi , (AW)i ) + k
q X
bi bj (Wi , Wj )
i,j=1
i=1
q hX
2
i
(bi aij + bj aji − bi bj )(Wi , Wj ) = kUk2 .
i,j=1
Applying this result to the scheme (1.10) we obtain (2.11)
kU n k = kU 0 k,
0 ≤ n ≤ M,
i.e. that our fully discrete method is conservative in the L2 sense.
8
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
3. Consistency and convergence In this section we shall study the consistency of the fully discrete scheme (1.10) and prove the error estimate (1.11). To this effect we shall first list some well-known algebraic properties of the Gauss–Legendre methods, [6], [9], that will be used in the sequel along with (2.1) and (2.2). (3.1) The q−stage Gauss–Legendre method is consistent of order 2q (i.e. has accuracy of order 2q when applied to an ordinary differential equation y ′ = f (t, y), where f and its partial derivatives of sufficiently high order with respect to y and t are smooth and bounded), and satisfies the following order (simplifying) conditions, [5]: q X
(3.1.1)
bi τiℓ−1 = ℓ−1 ,
i=1
q X
(3.1.2)
j=1
(3.1.3)
q X i=1
(3.2)
aij τjℓ−1 = τiℓ /ℓ,
1 ≤ ℓ ≤ 2q,
1 ≤ i, ℓ ≤ q,
bi τiℓ−1 aij = ℓ−1 bj (1 − τjℓ ),
1 ≤ j, ℓ ≤ q.
The q−stage Gauss–Legendre method corresponds to the q − th diagonal
Pad´e rational approximation to the exponential, i.e. if r(z) = 1 + zbT (I − zA)−1 e, e = (1, . . . , 1)T ∈ Rq ,
then r(z) is the q−th diagonal Pad´e approximant to exp(z). For the purposes of the proof of convergence we shall compare the solution U n of (1.10) to the elliptic projection W = W (t) ∈ Sh , 0 ≤ t ≤ T, of the solution u(t) of (1.1), defined as usual by (3.3)
(∆h W, ϕ) = −a(W, ϕ) = (∆u, ϕ) ∀ϕ ∈ Sh .
We shall denote the associated (time-independent) elliptic projection operator onto Sh (defined on H 2 ∩H01 ) by PI . In this notation, W (t) = PI u(t) and obviously W (j) (t) = PI u(j) (t); here and in the sequel v (j) (t) = (d/dt)j v(t). By our assumptions on Sh there follows that (3.4)
kv − PI vk + hkv − PI vk1 ≤ chr kvkr ,
∀v ∈ H r ∩ H01 .
Obviously (3.3) implies that kW (j) (t)k1 ≤ cku(j) (t)||1 ≤ cj , j ≥ 0. In addition, there exist constants cj such that kLh (t)W (j) (s)k ≤ cj , t, s ∈ [0, T ], j ≥ 0; this follows from
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
9
the observation that for any ϕ ∈ Sh Lh (t)W (j) (s), ϕ = α(∆h W (j) (s), ϕ) + (β(t)W (j) (s), ϕ) ≤ |α| k∆u(j)(s)k + |β(t)|∞ku(j) (s)k1 kϕk ≤ cku(j) (s)k2 kϕk.
(j)
In fact, if Lh (t) : Sh → Sh denotes the j−th time derivative of the operator Lh (t), (j) (j) (j) given by (Lh (t)ϕ, χ) = (∂tj β(·, t)ϕ, χ), j ≥ 1, for χ, ϕ ∈ Sh , i.e. by Lh = Bh , j ≥ 1, (i) we have kLh (t)W (j) (s)k ≤ ci ku(j)(s)k1 , i ≥ 1, j ≥ 0. Thus, we can generalize the previous estimate to (i)
kLh (t)W (j) (s)k ≤ cij ,
(3.5)
i, j ≥ 0,
t, s ∈ [0, T ],
and also note that (i)
(3.6)
kLh (t)k ≤ ci ,
i ≥ 1,
where k · k denotes here the L2 induced operator norm on Sh . We shall also make use of the following property of the elliptic projection, namely that for constants cij kLh (t)B (j) (s)W (i) (s)k ≤ cij ,
(3.7)
i, j ≥ 0,
t, s ∈ [0, T ],
which may be proved as follows. We have (3.8)
(j)
(j)
(j)
Lh (t)Bh (s)W (i) (s) = α∆h Bh (s)W (i) (s) + Bh (t)Bh (s)W (i) (s). (j)
(j)
Since for j ≥ 0, Bh (t)ϕ = P (∂tj β(t)ϕ), we have kBh ϕk ≤ cj kϕk for ϕ ∈ Sh . Hence in (3.8) (j)
kBh (t)Bh (s)W (i) (s)k ≤ cij .
(3.9)
Also for ϕ ∈ Sh we obtain, suppressing the dependence on s, (j)
(3.10)
(j)
−(∆h (t)Bh W (i) , ϕ) = a(Bh W (i) , ϕ) = a(P [β (j) W (i) ], ϕ)
= a(P [β (j)W (i) ] − β (j) u(i) , ϕ) + a(β (j) u(i) , ϕ).
We obviously have (3.11)
|a(β (j) u(i) , ϕ)| = |(∆(β (j) u(i) ), ϕ)| ≤ cij kϕk,
and by (1.3) (3.12)
|a(P [β (j) W (i) ] − β (j) u(i) , ϕ)| ≤ kP [β (j)W (i) ] − β (j) u(i) k1 kϕk1
≤ ch−1 kP [β (j) W (i) ] − β (j) u(i) k1 kϕk.
Now (3.13) kP [β (j)W (i) ] − β (j) u(i) k1 ≤ kP [β (j) (W (i) − u(i) )]k1 + kP [β (j)u(i) ] − β (j) u(i) k1 .
10
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
Since, as it may easily be seen from (1.3), (3.4), kP v −vk1 ≤ chr−1 kvkr for v ∈ H r ∩H01 , and since for β and u sufficiently smooth β (j) u(i) ∈ H r ∩ H01, there follows kP [β (j)(W (i) − u(i) )]k1 + kP [β (j) u(i) ] − β (j) u(i) k1
≤ ch−1 kP [β (j) (W (i) − u(i) )]k + cij hr−1 ≤ cij hr−1 .
These estimates, when substituted in (3.13), yield, in conjunction with (3.10)–(3.12), (j) that k∆h Bh W (i) k ≤ cij since r ≥ 2. Then (3.7) follows from (3.8) and (3.9). We now embark upon the proof of the main result of this section. For a function v defined on [0, T ] we generally denote v n = v(tn ). We first define, for the purposes of the proof of consistency, V n,m for 0 ≤ n ≤ M − 1, 1 ≤ m ≤ q and V n , 0 ≤ n ≤ M in Sh by (a) V 0 = W 0 , (b) V n,m = W n + ik (3.14)
q X
n,j amj Ln,j , h V
j=1
(c) V
n+1
n
= W + ik
q X
1 ≤ m ≤ q,
n,j bj Ln,j . h V
j=1
In Proposition 3.1 below we shall prove the consistency result, valid for u sufficiently smooth: (3.15) max kV n − W n k ≤ ck k min(2q,q+2) + hr . 0≤n≤M
If this holds, then a simple stability calculation, as the following theorem shows, gives the error bound (1.11): Theorem 3.1. Let u be sufficiently smooth and suppose that (3.15) holds. Then (3.16) max kU n − un k ≤ c k min(2q,q+2) + hr . 0≤n≤M
Proof. Let V n,m , V n be defined by (3.14) and let εn,m = U n,m −V n,m , εn = U n −V n , ζ n = U n − W n . Then (1.10) and (3.14) yield ε
n,m
n
= ζ + ik
q X
n,j amj Ln,j h ε ,
j=1
ε
n+1
n
= ζ + ik
q X
1 ≤ m ≤ q,
n,j bj Ln,j h ε .
j=1
The stability Lemma 2.2 gives then that kεn+1k = kζ n k. Hence, for 0 ≤ n ≤ M − 1, kζ n+1k ≤ kεn+1 k + kV n+1 − W n+1 k = kζ n k + kV n+1 − W n+1 k.
Therefore, by (3.15) kζ n k ≤ kζ 0 k + c(k min(2q,q+2) + hr ), 0 ≤ n ≤ M, and the result follows from (3.4), (1.5) and the triangle inequality. Hence our task is to prove consistency:
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
11
Proposition 3.1. If u is sufficiently smooth, (3.15) holds. Proof. We follow, up to a point, the technique of the consistency proof for Runge– Kutta discretizations of partial differential equations introduced in [11]. First define τij , 1 ≤ i ≤ q, j ≥ 0 by τi0 = 1, τij =
(3.17)
q X
aiℓ τℓ,j−1,
ℓ=1 j
⇐⇒ τij = (A e)i ,
1 ≤ i ≤ q, j ≥ 0,
1 ≤ i ≤ q, j ≥ 0.
Note that by (3.1.2) we may infer that τij = (τi )j /j!,
(3.18)
1 ≤ i ≤ q, 0 ≤ j ≤ q.
Also define, for 1 ≤ m ≤ q, 0 ≤ n ≤ M − 1 n
(3.19)
Λm W =
2q X
τmj k j W (j)n ,
j=0
en,m = V n,m − Λm W n .
(3.20)
We now make a preliminary observation. By (3.14) and (3.20) we have (3.21)
V
n+1
n
=W +
q X
m,j=1
bm (A−1 )mj (Λj W n − W n ) + bT A−1 en ,
where en = (en,1 , . . . , en,q )T ∈ (Sh )q . Using (3.19) and (3.17) we can write q X
m,j=1
bm (A−1 )mj (Λj W n − W n ) =
q X
bm (A−1 )mj
m,j=1
2q hX
(Aℓ e)j k ℓ W (ℓ)n
ℓ=1
i
2q X X T ℓ−1 ℓ (ℓ)n k ℓ W (ℓ)n /ℓ! (b A e)k W = = 2q
ℓ=1
ℓ=1
where in the last equality we used the identities (3.22)
bT Aℓ−1 e = 1/ℓ!,
1 ≤ ℓ ≤ 2q,
that follow from the fact that the rational approximation r(z), cf. (3.2), corresponding to the q−stage Gauss–Legendre method is an O(z 2q+1 ) approximation to exp(z) as z → 0. We conclude therefore by (3.21) and Taylor’s theorem that kV n+1 − W n+1 k ≤ ck 2q+1 + kbT A−1 en k. Hence, in order to prove (3.15) our preliminary observation is that it is sufficient to obtain an estimate of the form (3.23) kbT A−1 en k ≤ ck k min(2q,q+2) + hr .
12
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
Note that (3.20) when substituted in (3.14b) yields for 0 ≤ n ≤ M − 1 n,j
(3.24)
e
=E
n,j
+ ik
q X
n,ℓ ajℓ Ln,ℓ h e ,
ℓ=1
where we have put (3.25)
E
n,j
n
n
= −Λj W + W + ik
q X
1 ≤ j ≤ q,
n ajℓ Ln,ℓ h Λℓ W ,
ℓ=1
1 ≤ j ≤ q.
The first main step of the proof consists of some long intermediate calculations, inevitable in any order estimation of Runge–Kutta methods, aiming towards transforming the right-hand side of (3.25) to a form more suitable for our purposes. To simplify notation a bit, in the sequel by ϕ = O(k λ + hµ ) we shall mean that there exists a constant c > 0 independent of k and h such that kϕk ≤ c(k λ + hµ ) for k, h sufficiently small. First, using (3.25), (3.5), (3.17), (3.19) we have for 1 ≤ j ≤ q E
n,j
n
n
= −Λj W + W + ik + ik
q X
q X
ajℓ Ln,ℓ h
2q−1 X
τℓm k W
(m)n
2q−1 X
τℓm k m W (m)n
m
m=0
ℓ=1
(2q)n ajℓ τℓ,2q k 2q Ln,ℓ h W
ℓ=1
n
n
= −Λj W + W + ik q
+ ik
X
ajℓ (Ln,ℓ h
ℓ=1
(3.26)
−
q X
=−
Lnh )
m
τjm k W
(m)n
+
ajℓ (Ln,ℓ h
ℓ=1
−
Lnh )
m=1
+ ik
ℓ=1
=
I1n,j
+
2q−1 X
ℓ=1
n ajℓ (Ln,ℓ h − Lh )
I2n,j
+ O(k
2q+1
2q−1 X
),
m=0
+ O(k 2q+1 )
ajℓ τℓ,m−1 ik m Lnh W (m−1)n m
τℓm k W
τjm k m W (m)n − iLnh W (m−1)n
q X
(m)n
(m)n
m=0
2q
X
τℓm k W
q XX
m=1
X
m
m=0
q
=−
2q−1 X 2q
m=1
+ ik
m=0
ℓ=1
2q
X
ajℓ Ln,ℓ h
+ O(k 2q+1 )
τℓm k m W (m)n + O(k 2q+1 )
where (3.27)
I1n,j := −
2q X
m=1
τjm k m W (m)n − iLnh W (m−1)n ,
1 ≤ j ≤ q,
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
I2n,j := i
(3.28)
q X ℓ=1
In
I1n,j
n ajℓ (Ln,ℓ h − Lh )
2q−1 X m=0
τℓm k m+1 W (m)n ,
13
1 ≤ j ≤ q.
we have, using (1.8)
W (m)n − iLnh W (m−1)n = ∂tm−1 Wtn − iα∆h W n − iBhn W (m−1)n . Introducing now for 0 ≤ t ≤ T ψ(t) := Wt − ut − iβ(t) W (t) − u(t) , for which (3.4) shows that ψ (j) (t) = O(hr ), we can write using (3.3) and (1.1) that (3.29)
(3.30)
Wt − iα∆h W = P ψ + iBh (t)W.
Differentiating (3.30) with respect to t and using (3.29) we rewrite (3.27) as I1n,j (3.31)
2q X
= −i
m=1 2q X
τjm k m ∂tm−1 (Bhn W n ) − Bhn W (m−1)n + O(khr )
i h m−2 X m − 1 (m−1−ℓ)n Bh W (ℓ)n + O(khr ). τjm k = −i ℓ m=2 ℓ=0 m
Turning now to the term I2n,j we obtain by (3.6) and Taylor’s theorem (3.32)
I2n,j
=i
q X
2q 2q−1 X X µ (µ)n m (m−1)n ajℓ (τℓ k) Bh /µ! + O(k 2q+1 ). τℓ,m−1 k W µ=1
ℓ=1
m=1
Use of (3.5) yields for 1 ≤ ℓ ≤ q that 2q−1 2q
Zℓn
XX
:=
(j)n
(τℓ )j τℓ,m−1 k j+m Bh W (m−1)n /j!
j=1 m=1
2q
=
X
k
λ
λ−1 hX
(λ−m)n
(τℓ )λ−m τℓ,m−1 Bh
m=1
λ=2
i W (m−1)n /(λ − m)! + O(k 2q+1 ).
Therefore (3.32) gives for 1 ≤ j ≤ q (3.33)
I2n,j
=i
q X
ajℓ
ℓ=1
+ O(k
2q nX
k
λ
λ−1 hX
(λ−m)n
(τℓ )λ−m τℓ,m−1 Bh
m=1
λ=2 2q+1
).
io W (m−1)n /(λ − m)!
Summarizing, we obtain by (3.26), (3.31) and (3.33) for 1 ≤ j ≤ q (3.34)
E
n,j
=i
2q X
k
λ=2
+ O(k
λ
λ−1 hX
(λ−m)n
(δjm,λ − γjm,λ )Bh
m=1 2q+1
+ khr ),
i W (m−1)n /(λ − m)!
where, for 1 ≤ j ≤ q, λ = 2, . . . , 2q, m = 1, . . . , λ − 1: (3.35)
δjm,λ
=
q X ℓ=1
ajℓ (τℓ )λ−m τℓ,m−1 ,
γjm,λ = τjλ (λ − 1)!/(m − 1)!.
14
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
We immediately observe, using (3.18), that for λ = 1, . . . , q + 1, 1 ≤ m ≤ λ − 1 δjm,λ
=
q X
ajℓ (τℓ )
λ−1
ℓ=1
q X /(m − 1)! = (λ − 1)!/(m − 1)! ajℓ τℓ,λ−1 = γjm,λ . ℓ=1
Therefore, we finally conclude that if q ≥ 2, (3.34) yields for 1 ≤ j ≤ q (3.36)
E
n,j
=i
2q X
kλ
m=1
λ=q+2
+ O(k
λ−1 hX
2q+1
(λ−m)n
(δjm,λ − γjm,λ )Bh
+ khr ),
i W (m−1)n /(λ − m)!
while, simply, if q = 1 E n,1 = O k(k 2 + hr ) .
(3.37)
Note that (3.37) used in conjunction with
max ken,j k ≤ c max kE n,j k,
(3.38)
1≤j≤q
1≤j≤q
(which follows from the stability estimate of Lemma 2.1 applied to the equations (3.24)) gives the desired estimate (3.23) in the case q = 1. Therefore we henceforth concentrate on the cases q ≥ 2 for which E n,j is given by (3.36). To this end let for 1 ≤ j ≤ q, (3.39)
ϕ
n,j
=i
2q X
k
λ=q+2
λ
λ−1 hX
(δjm,λ
m=1
−
(λ−m)n γjm,λ )Bh W (m−1)n /(λ
i − m)! .
Then (3.36) is written as (3.40)
E n,j = ϕn,j + O(k 2q+1 + khr ),
1 ≤ j ≤ q.
Obviously, (3.5) gives that ϕn,j = O(k q+2) since δjm,λ − γjm,λ are not zero in general if λ ≥ q + 2. (3.40) then gives that E n,j = O(k q+2 + khr ) and (3.38) implies in turn that en,j = O(k q+2 + khr ) thus yielding the estimate bT A−1 en = O(k(k q+1 + hr )). (This concludes essentially the application of the idea of the consistency proof of [11] to the case at hand). The second step in the proof is the improvement of the q + 1 exponent of k to the better value q + 2. For this purpose we use the fact that we must actually estimate not the individual en,j but their particular linear combination bT A−1 en =
q X
bi (A−1 )ij en,j .
i,j=1
First note that defining (3.41)
E˜ n,j = E n,j − ϕn,j
e˜n,j = en,j − ϕn,j ,
we may rewrite (3.24) as q q X X n,ℓ n,ℓ n,j n,j ajℓ Ln,ℓ e˜n,ℓ . ajℓ L ϕ + ik (3.42) e˜ = E˜ + ik h
h
ℓ=1
ℓ=1
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
15
Recall from (3.40), (3.41) that E˜ n,j = O(k 2q+1 + khr ). Also, (3.7) and (3.39) yield that n,j q+2 kLn,j , 1 ≤ j ≤ q. Hence, the stability estimate of Lemma 2.1, applied to h ϕ k ≤ ck n,j (3.42), yields e˜ = O(k q+3 + khr ). Therefore, it follows from (3.41) that (3.43)
bT A−1 en = O(k(k q+2 + hr )) + bT A−1 ϕn ,
where ϕn = (ϕn,1 , . . . , ϕn,q )T ∈ (Sh )q . Hence, the desired estimate (3.23) will follow from the fact that actually bT A−1 ϕn = 0,
(3.44)
which is a consequence of (3.39) and the following cancellation property of the Gauss– Legendre methods that we state as a separate lemma: Lemma 3.1. Let δjm,λ , γjm,λ be defined for 1 ≤ j ≤ q, λ = 2, . . . , 2q, m = 1, . . . , λ − 1 by T T (3.35) and denote δ m,λ = δ1m,λ , . . . , δqm,λ and γ m,λ = γ1m,λ , . . . , γqm,λ ∈ Rq . Then
(3.45)
bT A−1 (δ m,λ − γ m,λ ) = 0.
Proof. For λ ≤ q + 1, 1 ≤ m ≤ λ − 1, we have already established by the simple calculation following (3.35) that δ m,λ = γ m,λ . Hence we restrict attention to the interesting case q + 2 ≤ λ ≤ 2q, 1 ≤ m ≤ λ − 1. Define for the purposes of this lemma T ∈ Rq×q as T = diag (τ1 , . . . , τq ) — no confusion with the T of (1.1) will arise — and note that (3.1.1), (3.1.2) can be written equivalently as (3.1.1′ )
bT T ℓ−1 e = 1/ℓ,
(3.1.2′ )
AT ℓ−1 e = T ℓ e/ℓ,
1 ≤ ℓ ≤ 2q, 1 ≤ ℓ ≤ q,
respectively; (3.1.2′ ) implies that (3.1.2′′ )
Aℓ e = T ℓ e/ℓ!,
0 ≤ ℓ ≤ q.
Now, (3.17) gives that δ m,λ = AT λ−m Am−1 e. On the other hand, using (3.22), since λ ≤ 2q, we have bT A−1 γ m,λ = (λ − 1)!bT Aλ−1 e/(m − 1)! = (λ(m − 1)!)−1 . Hence, to show (3.45) it suffices to establish (3.46)
bT T λ−m Am−1 e = (λ(m − 1)!)−1 ,
q + 2 ≤ λ ≤ 2q, 1 ≤ m ≤ λ − 1.
Obviously for each λ, q + 2 ≤ λ ≤ 2q, (3.46) holds for 1 ≤ m ≤ q + 1; this follows from (3.1.2′′ ) and (3.1.1′) that yield bT T λ−m Am−1 e = bT T λ−1 e/(m − 1)! = (λ(m − 1)!)−1 . Hence we focus on the case q + 2 ≤ λ ≤ 2q, q + 2 ≤ m ≤ λ − 1. Using (3.1.2′′ ) again gives, since m − 1 ≥ q + 1, (3.47)
bT T λ−m Am−1 e = bT T λ−m Am−1−q T q e/q!.
For integers k ≥ 1, 1 ≤ ℓ ≤ q define (3.48)
F (ℓ, k) = bT T ℓ−1 Ak T q e
16
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
and note that by (3.1.3) T
F (ℓ, k) = b T
ℓ−1
k−1
A(A
q
T e) =
q X
bi τiℓ−1
= ℓ−1
(3.49)
X j=1
=ℓ
−1
aij (Ak−1 T q e)j
j=1
i=1
q
q X
bj (1 − τjℓ )(Ak−1 T q e)j
bT Ak−1 T q e − bT T ℓ Ak−1 T q e
= ℓ−1 F (1, k − 1) − F (ℓ + 1, k − 1) ,
k ≥ 1, 1 ≤ ℓ ≤ q.
We can now calculate directly using (3.1.1), (3.1.3) for integer s, such that 1 ≤ s ≤ q−1, that F (s, 1) =
q X
bi τis−1
= s−1
aij τjq
q
hX j=1
bj τjq −
=s
−1
q X j=1
j=1
i=1
(3.50)
q X
q
X j=1
bj (1 − τjs )τjq
i −1 bj τjs+q = (q + 1)(s + q + 1) .
We claim now that (3.51)
−1 F (ℓ, k) = (q + 1) · · · (q + k)(q + ℓ + k)
for k ≥ 1, ℓ ≥ 1, ℓ + k ≤ q.
Indeed (3.50) shows that (3.51) holds for k = 1. For the inductive step assume that (3.51) holds for all k ≤ k ′ and 1 ≤ ℓ ≤ q − k. Then, for ℓ ≤ q − k ′ − 1 we have by (3.49) and the inductive hypothesis that F (ℓ, k ′ + 1) = ℓ−1 F (1, k ′) − F (ℓ + 1, k ′ ) −1 −1 = (q + 1) · · · (q + k ′ )(q + k ′ + 1) − (q + 1) · · · (q + k ′ )(q + k ′ + ℓ + 1) /ℓ −1 = (q + 1) · · · (q + k ′ )(q + k ′ + 1)(q + ℓ + k ′ + 1) . This completes the inductive step and shows the validity of (3.51). Using now (3.48), (3.51) and (3.47) we obtain with k := m − 1 − q ≥ 1, ℓ := λ − m + 1 ≥ 1, since k + ℓ = λ − q ≤ q in the region of interest, that for q + 2 ≤ λ ≤ 2q, q + 2 ≤ m ≤ λ − 1, −1 −1 bT T λ−m Am−1 e = q!(q + 1) · · · (m − 1)λ = λ(m − 1)! . This identity completes the proof of the validity of (3.46).
As a consequence of (3.45) and (3.39), (3.44) holds and the proof of Theorem 3.1 is now complete. 4. Practical implementation of the two-stage scheme In this section we shall study questions related to the efficient implementation of the fully discrete scheme generated by the two-stage Gauss–Legendre method which is
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
17
given of course by the tableau, [9], [6]:
(4.1)
A τ T
b
1 4 1 4
=
1 4
1 √ 2 3
1 4
1 √ 2 3
+
−
1 2
1 2 1 2
− +
1 √ 2 3 1 √ 2 3
.
1 2
With these values of aij , bi , τi the method is (a) U 0 = u0h , for n = 0, 1, . . . , M − 1 : (4.2)
(b) U
n,m
n
= U + ik
2 X
n,j amj Ln,j , h U
m = 1, 2,
j=1
(c) U
n+1
n
= U + ik
2 X
n,j bj Ln,j . h U
j=1
n,j Note that eliminating the Ln,j from (4.2b and c) and using the particular values h U of the constants aij and bi from (4.1), we may write (4.2c) simply as √ (4.2c′ ) U n+1 = U n + 3(U n,2 − U n,1 ).
In sections 2 and 3 we proved that for each 0 ≤ n ≤ M (4.2) has a unique solution, that the resulting scheme is L2 −conservative and that it satisfies the optimal in space and time estimate (4.3) max kU n − un k ≤ c k 4 + hr . 0≤n≤M
Let d = dim Sh . To determine U n,1 , U n,2 from (4.2b) one should solve, after choosing a basis for Sh , the 2d × 2d linear system J n Un = F n ,
(4.4) where (4.5)
Jn = Jn (tn,1 , tn,2 ) =
! n,2 I − ika11 Ln,1 −ika L 12 h h n,2 , −ika21 Ln,1 I − ika L 22 h h
Un = (U n,1 , U n,2 )T ,
Fn = (U n , U n )T .
With the aim of solving only (sparse) d × d systems of linear equations at each time step, we shall uncouple the two equations in (4.4) borrowing an idea from [3]. We first write (4.4) equivalently as (4.6)
J∗n Un = (J∗n − Jn )Un + Fn ,
where (4.7)
J∗n = Jn (t∗n , t∗n ),
t∗n = (tn,1 + tn,2 )/2 = tn + k/2,
18
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
the advantage being now that the operators Lh in the entries of J∗n are evaluated at the same point t∗n and that, in particular, they commute. This enables us to compute the solution Y = (Y1 , Y2)T ∈ (Sh )2 of systems of the form J∗n Y = Z,
(4.8)
given Z = (Z1 , Z2 )T ∈ (Sh )2 , “explicitly” as
∆∗n Y1 = Z1 + ikL∗n (a12 Z2 − a22 Z1 ),
(4.9)
∆∗n Y2 = Z2 + ikL∗n (a21 Z1 − a11 Z2 ),
where ∆∗n = I − k(iL∗n )/2 + k 2 (iL∗n )2 /12, L∗n = Lh (t∗n ). It is easily seen that ∆∗n is invertible since iL∗n is normal with purely imaginary eigenvalues and the polynomial 1 − z/2 + z 2 /12 has no roots on the imaginary axis. Systems with operator ∆∗n like the ones in (4.9) can then in turn be solved by the complex analog of the procedure proposed for Pad´e diagonal methods in [2]. Consider e.g. the first equation of (4.9) and, putting W1 = a12 Z2 − a22 Z1 , R1 = W1 + Z1 /2 − ikL∗n Z1 /12, Φ1 = Y1 − Z1 , rewrite it as ∆∗n Φ1 = ikL∗n R1 .
(4.10)
√ Since ∆∗n = (I − iµkL∗n )(I − i¯ µkL∗n ), where µ = 1/4 − i 3/12, we may rewrite (4.10), letting H n = I − iµkL∗n , K n = I − i¯ µkL∗n , as H n K n Φ1 = i(H n − K n )R1 /(2 Im µ), from which Φ1 = i[(K n )−1 − (H n )−1 ]R1/(2 Im µ).
(4.11)
To determine therefore Φ1 (and hence Y1 ) one must form R1 in the right-hand side of (4.11) and solve two complex linear systems of size d × d with operators K n and H n noting that at the level of matrix-vector operations the corresponding matrices will be sparse if a finite element basis consisting of functions of small support is chosen for Sh . Similarly for Φ2 with obvious parallelicity duly noted. In order to take advantage of the fact that systems of the form (4.8) can be “easily” solved in the manner outlined above, we shall solve the original system (4.6) by the simplest iterative method that its form immediately suggests. Let jn , n = 0, 1, . . . , M, be given positive integers — in practice jn = 1 or jn = 2 — representing the number of iterations that will be performed at step n to solve (4.6). For jn , n = 0, 1, . . . , M we shall compute Ujn approximating U n as follows: Let Uj00 = U 0 = u0h (e.g. take j0 = 0). For n = 0, 1, . . . , M: Compute suitable starting values U0n,1 , U0n,2 . Compute Ujn,m , m = 1, 2 by the iterative scheme: n+1 For j = 0, 1, . . . , jn+1 − 1 solve (4.12)
J∗n
n,1 Uj+1
Ujn,2 +1
!
= (J∗n − Jn )
Ujn,1 Ujn,2
!
+
Ujnn Ujn,2 n
!
.
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
19
√
− Ujn,1 ). Note that when j = jn+1 − 1 in the inner Define Ujn+1 = Ujnn + 3(Ujn,2 n+1 n+1 n+1 (j) loop in (4.12) there are important savings in the number of operations: since only the difference Ujn,2 − Ujn,1 is finally needed, we must solve directly only one system of n+1 n+1 ∗n the form ∆ (Φ2 − Φ1 ) = ikL∗n (R2 − R1 ). Hence it is important to try to get by with only jn = 1 iteration at each time step. We shall analyze the convergence of the iterative scheme in (4.12) and at the same time demonstrate the attendant stability of the resulting new fully discrete method. To this end assume for the time being that we have available approximations Ujnn ∈ Sh to U n for n = 0, 1, 2, 3, 4 such that (4.13) kU n − Ujnn k ≤ cn k 4 + hr , n = 0, 1, 2, 3, 4.
For n ≥ 3, given Ujnn approximating U n , we shall compute the required starting values U0n,m , m = 1, 2, for the inner (j) loop in (4.12) using extrapolation from previous values, i.e. as 3 X n,m µmλ Ujn−λ , m = 1, 2, 3 ≤ n ≤ M − 1, (4.14) U0 = n−λ λ=0
where the constants µmλ , m = 1, 2, λ = 0, 1, 2, 3, will be computed by letting p(t) be the cubic polynomial interpolating to the values y n−λ = y(tn−λ) of a smooth function y(t) at the points tn−λ , λ = 0, 1, 2, 3, and setting n
p(t + τm k) =
3 X
µmλ y n−λ. m = 1, 2.
λ=0
Proposition 4.1. Let u be sufficiently smooth, U n , 0 ≤ n ≤ M, be the solution of (4.2) and suppose that there exist Ujnn ∈ Sh , 0 ≤ n ≤ 4, such that (4.13) holds. For 4 ≤ n ≤ M −1 define Ujn+1 by the scheme (4.12) with starting values U0n,m as in (4.14). n+1 Then if jn+1 = 2 for all n, there exists k0 > 0 such that for k ≤ k0 (4.15) max kU n − Ujnn k ≤ c k 4 + hr . 0≤n≤M
If jn+1 = 1 for all n, then, given 0 < ε ≤ 1, there exists kε > 0 such that for k ≤ kε there holds (4.16) max kU n − Ujnn k ≤ c k 4−ε + hr . 0≤n≤M
In (4.15) and (4.16) c is a constant independent of k, h and ε.
Proof. Given Ujnn , let U˜ n,m , m = 1, 2, denote the exact solution of the system (4.17)
U˜ n,m = Ujnn + ik
2 X
˜ n,j , amj Ln,j h U
j=1
which we can write using (4.5) as (4.18)
Jn
U˜ n,1 U˜ n,2
!
=
! Ujnn . Ujnn
m = 1, 2,
20
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
First we prove a preliminary estimate that implies the convergence of the sequence ˜ n,m . From (4.12) and (4.18) suppressing n we let Yj,m = Ujn,m , j = 0, 1, 2, . . . to U Ujn,m − U˜ n,m , Yj = (Yj,1, Yj,2)T and obtain J∗n Yj+1 = (J∗n − Jn )Yj =: Zj = (Zj,1, Zj,2)T , i.e. Yj+1 = ikAL∗n Yj+1 + Zj ,
(4.19)
where L∗n = diag (L∗n , L∗n ) on (Sh )2 . There follows by Lemma 2.1 that |||Yj+1||| ≤ c|||Zj |||. On the other hand, using e.g. (3.6), we obtain kZj,mk ≤ ck 2 kYj,1k+kYj,2k), m = 1, 2 and conclude therefore that (4.20)
n,m ˜ n,m k ≤ γk 2 max kU n,m − U˜ n,m k, max kUj+1 −U j
m=1,2
m=1,2
j = 0, 1, . . . , jn+1 − 1,
where γ is a constant independent of h, k and the choice of jn . We next estimate the difference kU0n,m − U˜ n,m k. (In what follows we let n ≥ 4). Using (4.14) we have for m = 1, 2
(4.21)
3 X n,m n,m n,m n,m n,m n,m n,m ˜ ˜ − U ) + (U −u )+ u − µmλ un−λ U − U0 = (U λ=0
+
3 hX λ=0
µmλ u
n−λ
−U
n−λ
i
+
3 hX
µmλ U
n−λ
λ=0
−
Ujn−λ n−λ
i .
First, since U˜ n,m satisfies (4.17) and U n,m (4.2b) we obtain, again by stability (Lemma 2.1), that ˜ n,m − U n,m k ≤ ckU n − U n k. max kU jn
(4.22)
m=1,2
Next, recalling the definition of V n,m from (3.14), write, with W n,m = W (tn,m), (4.23)
U n,m − un,m = (U n,m − V n,m ) + (V n,m − W n,m) + (W n,m − un,m)
and observe that by (3.14b), (4.2b), Lemma 2.1 and (4.3) (4.24)
max kU n,m − V n,m k ≤ ckU n − W n k ≤ c(k 4 + hr ).
m=1,2
Now from (3.19), (3.20) V n,m − W n,m = en,m + (Λm W n − W n,m). By (3.38) and (3.36) we have ken,m k ≤ c(k 4 + khr ), m = 1, 2, whereas (3.18) and Taylor’s theorem yield for m = 1, 2 4
X
n n,m
τmj k j W (j)n − W n,m ≤ ck 3 . kΛm W − W k = j=0
We conclude that (4.25)
kV n,m − W n,m k ≤ c(k 3 + khr ),
m = 1, 2,
and, therefore, by (4.23)–(4.25) that for k ≤ 1, (4.26)
kU n,m − un,m k ≤ c(k 3 + hr ),
m = 1, 2.
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
21
Finally, the definition of the extrapolation procedure yields 3
n,m X
u µmλ un−λ ≤ ck 4 . −
(4.27)
λ=0
We conclude, by (4.21), (4.22), (4.26) and (4.27) that there exist positive constants η and θ, independent of k, h and the choice of jn such that (4.28)
˜ n,m − U0n,m k ≤ θ(k 3 + hr ) + ηkU n − U n k + max kU jn
m=1,2
3 X λ=0
k, µλ kU n−λ − Ujn−λ n−λ
where µλ = maxm=1,2 |µmλ |, 0 ≤ λ ≤ 3. We come now to the main part of the proof which is an induction √ 2 −1 step on n. First we treat the case jn = 2 for all n ≥ 4. We let k ≤ k0 = 2 3 γ where γ is as in (4.20). Assume that for 4 ≤ m ≤ n, (H1)
kU m − Ujmm k ≤ cm (k 4 + hr ),
where cm are positive constants satisfying (H2)
cm = [1 + k 3 (η + µ0 )]cm−1 + k 3 (µ1 cm−2 + µ2 cm−3 + µ3 cm−4 ) + θk 2 .
Clearly, it may be arranged, by taking c3 or c4 large enough in (4.13), that (H1)–(H2) hold for n = 4. Define U˜ n+1 ∈ Sh , conformal to the notation in (4.17) as √ ˜ n+1 = U n + 3(U˜ n,2 − U˜ n,1 ) (4.29) U jn and split (4.30)
U n+1 − Ujn+1 = (U n+1 − U˜ n+1 ) + (U˜ n+1 − Ujn+1 ). n+1 n+1
Since the time stepping procedure is conservative in the L2 sense, subtracting (4.17) from (4.2b) and (4.29) from (4.2c′ ) we obtain (4.31)
˜ n+1 k = kU n − U n k. kU n+1 − U jn
On the other hand, using (4.29) and (4.12) we have √ ˜ n+1 − U n+1 k = k 3(U˜ n,2 − U n,2 ) − (U˜ n,1 − U n,1 )k kU jn+1 jn+1 jn+1 √ n,m n,m ˜ − Ujn+1 k ≤ 2 3 max kU (4.32) m=1,2 √ ˜ n,m − U n,m k ≤ 2 3 γ 2 k 4 max kU 0 m=1,2
where in the last inequality we used the fact that jn+1 = 2 and (4.20). We conclude therefore by (4.30)–(4.32) that √ 2 4 n n ˜ n,m − U n,m k. 3 γ k max kU (4.33) kU n+1 − Ujn+1 k ≤ kU − U k + 2 0 jn n+1 m=1,2
√
Hence, (4.33) and (4.28) yield, if k ≤ k0 = 2 3 γ kU n+1 − Ujn+1 k ≤ (1 + ηk 3 )kU n − Ujnn k + k 3 n+1
3 X λ=0
2 −1
k + θk 2 (k 4 + khr ). µλ kU n−λ − Ujn−λ n−λ
22
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
Therefore, using the induction hypothesis (H1) we obtain kU n+1 − Ujn+1 k n+1 3 ≤ 1 + k (η + µ0 ) cn + k 3 (µ1 cn−1 + µ2 cn−2 + µ3 cn−3 ) + θk 2 (k 4 + khr ).
Hence, (H1)–(H2) hold for m = n+1 as well and the inductive step is complete. Clearly the constants cn are uniformly bounded for 0 ≤ n ≤ M by a constant independent of k and h; (4.15) follows. We examine now the consequences of taking jn = 1 for all n ≥ 4. Given 0 < ε ≤ 1 we √ −1/ε shall assume that k ≤ kε = 2 3 γ for reasons that will become apparent below. Our induction hypotheses that replace (H1), (H2) are now that for 4 ≤ m ≤ n, kU m − Ujmm k ≤ cm (k 4−ε + hr ),
(H1′ )
where cm are positive constants given by (H2′ )
cm = [1 + k 2−ε (η + µ0 )]cm−1 + k 2−ε (µ1 cm−2 + µ2 cm−3 + µ3 cm−4 ) + θk.
The verification of the inductive step follows the lines of the previous proof: (4.29) ˜ n+1 − U n+1 k ≤ to (4.31) still hold of course but since jn+1 = 1, (4.32) becomes kU jn+1 √ ˜ n,m − U0n,m k. Consequently we have kU n+1 − U n+1 k ≤ kU n − 2 3 γk 2 maxm=1,2 kU jn+1 √ −1/ε √ n,m 2 n,m n ˜ − U0 k, and therefore for k ≤ kε = 2 3 γ U k + 2 3 γk maxm=1,2 kU we jn
obtain by (H1′ )
h √ 2 n n 3 γk θ(k 3 + hr ) + ηkU n − Ujnn k k + 2 kU n+1 − Ujn+1 k ≤ kU − U jn n+1 +
3 X λ=0
≤
µλ kU n−λ − Ujn−λ k n−λ
i
1 + (η + µ0 )k 2−ε cn + k 2−ε (µ1 cn−1 + µ2 cn−2 + µ3 cn−3 ) (k 4−ε + hr )
+ θk(k 4−ε + hr ),
verifying that (H1′ ), (H2′ ) hold for m = n + 1 as well. Obviously the constants cn can be made uniformly bounded in n by a constant independent of ε, albeit larger than the corresponding constant in the jn = 2 case; we conclude that (4.16) holds. It is easy to construct Ujnn , 0 ≤ n ≤ 4, such that (4.13) is valid. For n = 0 we already have stipulated that Uj00 = U 0 . For n = 1 only the previous value U 0 is available. We set U00,1 = U00,2 = U 0 in (4.12) and generate the sequence Uj0,m taking j1 = 2. This will suffice since the analog of (4.27) is an O(k) bound implying that in the analog of (4.28) maxm=1,2 kU˜ 0,m − U00,m k = O(k + hr ). Similarly for n = 2, 3, 4 it suffices to generate Ujnn from (4.12) with jn = 1, computing for each n, U0n−1,m as the appropriate linear combination of the previous values Uj00 , Uj11 , . . . , Ujn−1 using the Lagrange interpolating n−1 polynomial of degree n − 1. In practice we noticed that taking cubic polynomial extrapolation to generate the starting values for n ≥ 4 and just jn = 1 was generally sufficient to preserve the overall order of accuracy and stability of the scheme. We report in [1] these and other relevant numerical experiments that we performed with the method, including experiments in
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
23
which the operator J∗n is not evaluated at every time step but rather every m∗ > 1 time steps. Our experiments indicate that it is possible to take in many interesting examples m∗ equal to, say, 20 and jn = 2 (the m∗ “large”, jn = 1 combination was unstable for some hard to integrate problems) and still preserve the overall order of accuracy and stability of the scheme. Note added in proof : For more recent work related to error estimates for the Nonlinear Schr¨odinger Equation we refer the reader to [24], [25]. References 1. G.D. Akrivis and V.A. Dougalis, On a conservative, high-order accurate finite element scheme for the ‘parabolic’ equation, to appear in the Proceedings for the 2nd IMACS symposium on computational acoustics, Princeton University, 15-17 March 1989. 2. G.A. Baker, J.H. Bramble and V. Thom´ee, Single step Galerkin approximations for parabolic problems, Math. Comp. 31 (1977) 818–847. 3. J.L. Bona, V.A. Dougalis, O.A. Karakashian and W. McKinney, Conservative high order schemes for the generalized Korteweg-de Vries equation, to appear. 4. A. Broc´ehn, Galerkin methods for approximation of solutions of second order partial differential equations of Schr¨ odinger type, Ph.D. Thesis, University of G¨oteborg, 1980. 5. J.C. Butcher, Implicit Runge-Kutta processes, Math. Comp. 18 (1964) 50–64. 6. J.C. Butcher, The numerical analysis of ordinary differential equations; Runge-Kutta methods and general linear methods, John Wiley, Chichester, 1987. 7. M. Crouzeix, Sur la B-stabilit´e des m´ethodes de Runge-Kutta, Numer. Math. 32 (1979) 75–82. 8. M. Crouzeix and V. Thom´ee, On the discretization in time of semilinear parabolic equations with nonsmooth initial data, Math. Comp. 49 (1987) 359–377. 9. K. Dekker and J.G. Verwer, Stability of Runge-Kutta methods for stiff nonlinear differential equations, North Holland, Amsterdam, 1984. 10. E. Dendy, Jr., An alternating direction method for Schr¨ odinger’s equation, SIAM J. Numer. Anal. 14 (1977) 1028–1032. 11. V.A. Dougalis and O.A. Karakashian, On some high order accurate fully discrete Galerkin methods for the Korteweg-de Vries equation, Math. Comp. 45 (1985) 329–345. 12. O.A. Karakashian and W. McKinney, On optimal high order in time approximations for the Korteweg-de Vries equation, Math. Comp. (to appear). 13. D. Lee and S.T. McDaniel, Ocean acoustic propagation by finite difference methods, Comput. Math. Appl. 14 No. 5. 14. D. Lee, R.L. Sternberg and M.H. Schultz, eds., Computational acoustics: wave propagation, Proceedings of the 1st IMACS symposium on computational acoustics, New Haven, 6–8 August 1986, vols.1,2, North Holland, Amsterdam, 1988. 15. J.L. Lions and E. Magenes, Probl`emes aux limites non homog`enes et applications, vol. 2, Dunod, Paris, 1968. 16. A. Quarteroni, Mixed approximations of evolution problems, Comput. Meths. Appl. Mech. Engrg. 24 (1980) 137–163. 17. J.M. Sanz-Serna, Methods for the numerical solution of the nonlinear Schr¨ odinger equation, Math. Comp. 43 (1984) 21–27. 18. J.M. Sanz-Serna and J.G. Verwer, Conservative and nonconservative schemes for the solution of the nonlinear Schr¨ odinger equation, IMA J. Numer. Anal. 6 (1986) 25–42. 19. M.H. Schultz and D. Lee, eds., Computational ocean acoustics, Invited lectures from the workshop held at Yale University, 1-3 August 1984, Comput. Math. Appl. 11 (1985) Nos 7-8.
24
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
20. F.D. Tappert, The parabolic approximation method. In: Wave propagation and underwater acoustics, J.B. Keller and J.S. Papadakis, eds., pp. 224–287, Lecture Notes in Physics v. 70, SpringerVerlag, Berlin-Heidelberg, 1977. 21. V. Thom´ee, Convergence estimates for semi-discrete Galerkin methods for initial-value problems. In: Numerische, insbesondere approximations-theoretische Behandlung von Funktionalgleichungen, R.Ansorge and W. Tornig, eds., pp. 243–262, Lecture Notes in Mathematics v. 333, SpringerVerlag, Berlin-Heidelberg, 1973. 22. V. Thom´ee, Galerkin finite element methods for parabolic problems, Lecture Notes in Mathematics v. 1054, Springer-Verlag, Berlin-Heidelberg, 1984. 23. L.B. Wahlbin, A dissipative Galerkin method for the numerical solution of first order hyperbolic equations. In: Mathematical aspects of finite elements in partial differential equations, C. de Boor, ed., pp. 147–169, Academic Press, New York, 1974. 24. G.D. Akrivis, V.A. Dougalis and O.A. Karakashian, On fully discrete Galerkin methods of secondorder temporal accuracy for the Nonlinear Schr¨ odinger Equation, to appear in Numer. Math. 25. O. Karakashian, G.D. Akrivis and V.A. Dougalis, On optimal-order error estimates for the Nonlinear Schr¨ odinger Equation, to appear. Mathematics Department, University of Crete, Iraklion, Greece Mathematics Department, University of Crete, Iraklion, Greece
Abstract. We construct and analyze efficient fully discrete Galerkin type methods that are of high order of accuracy and conservative in the L2 sense for approximating the solution of a form of the linear Schr¨odinger equation with a time-dependent coefficient, found e.g. in underwater acoustics. The time stepping procedures are based on the class of implicit Runge–Kutta methods known as the q−stage Gauss– Legendre schemes. L2 error estimates are proved that are of optimal order in space and of temporal order q + 2. An iterative procedure at each time step for the efficient implementation of the two-stage scheme is proposed and analyzed. ` ` . On construit et analyse des m´ RESUM E ethodes totalement discr`etes du type Galerkin, qui sont L2 −conservatives et d’order arbitraire, pour approcher la solution d’une forme de l’´equation lin´eaire de Schr¨odinger avec un coefficient qui d´epend du temps, trouv´ee par exemple dans l’acoustique sous-marine. La proc´edure de discr´etisation en temps est bas´ee sur la classe des m´ethodes implicites de Runge–Kutta connues comme les sch´emas de Gauss–Legendre `a q pas interm´ediaires. On obtient des estimations dans L2 pour les erreurs, qui sont d’order optimal en espace es d’order q + 2 en temps. On propose et analyse aussi une proc´educe it´erative pour r´esoudre les syst`emes lin´eaires ` a chaque pas de temps pour une application efficace de sch´ema Gauss– Legendre ` a deux pas interm´ediaires.
Dedicated to Professor Dr. G. H¨ ammerlin on the occasion of his 60th birthday, July 31, 1988.
1. Introduction In this paper we shall study conservative numerical methods of high order of accuracy for approximating the solution of the following initial- and boundary-value problem for a partial differential equation of the Schr¨odinger type. Let Ω be a bounded domain in RN with smooth boundary ∂Ω and let 0 < T < ∞ be given. We seek a complex-valued ¯ × [0, T ], satisfying: function u = u(x, t), (x, t) ∈ Ω ut = iL(t)u := i α∆u + β(x, t)u in Ω × [0, T ],
(1.1)
u=0
on ∂Ω × [0, T ], ¯ in Ω,
u(x, 0) = u0 (x)
where α is a given nonzero real number, β is a given smooth real-valued function on ¯ × [0, T ] and u0 is a given smooth complex-valued function on Ω. ¯ We shall assume Ω Work supported by the Institute of Applied and Computational Mathematics of the Research Center of Crete-FORTH. 1
2
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
that the data of (1.1) are smooth and compatible enough to ensure that the problem possesses a unique and smooth enough for our purposes solution; cf. e.g. [15, Chapter 5, Section 12] and [4] for relevant existence, uniqueness and regularity results. This form of the Schr¨odinger equation occurs, for example for N = 1, in underwater acoustics as “parabolic approximation” to the Helmholtz equation, cf. e.g. [20], posed with a variety of types of boundary conditions. For simplicity we consider here the case of homogeneous Dirichlet boundary conditions. In (1.1) the Laplacian ∆ could have as well been replaced by a second-order, symmetric, uniformly positive definite elliptic ¯ with space-dependent coefficients with no complications in the error operator on Ω estimates. We shall discretize (1.1) in space by a Galerkin-finite element type method and in time by a class of implicit Runge–Kutta schemes of arbitrary order, known as the Gauss–Legendre collocation type methods. We shall estimate the error of the fully discrete approximations in L2 and point out efficient ways for implementing the methods. Many finite difference and spectral schemes, usually of second-order temporal accuracy, have been proposed in the literature for problems such as (1.1); see e.g. the survey [13], the collections of papers in [19] and [14] and their references. Galerkin-finite element methods have also been considered. For the linear Schr¨odinger equation with time-independent coefficients early error estimates for semidiscrete approximations may be found in [21], [23]. Semidiscrete and fully discrete Galerkin approximations with Runge–Kutta time stepping have been analyzed in [4] in the general case where the Laplacian in the right-hand side of the p.d.e. is replaced by a second-order elliptic operator with space- and time-dependent coefficients. For alternative approaches based on separating real and imaginary parts, cf. [10], [16]; in this paper we shall discretize (1.1) directly using complex arithmetic. Among the growing literature on Galerkin type methods for nonlinear Schr¨odinger equations, some error estimates are shown in [17], where second-order time stepping procedures coupled with finite difference or Galerkin type space discretizations are analyzed; for computations with such schemes cf. e.g. [18]. We next introduce notation that will be used in the sequel. For integral s ≥ 0 let H s = H s (Ω) denote the usual, complex (Hilbert) Sobolev spaces with corresponding norm k · ks . For f, g ∈ L2 = H 0 , let (f, g) =
Z
f (x)g(x) dx Ω
be their L2 inner product; here the overbar denotes complex conjugation. Let k · k denote the associated L2 norm and | · |∞ be the norm of L∞ = L∞ (Ω). As usual, H01 will consist of the elements of H 1 that vanish on ∂Ω in the sence of trace. We shall discretize (1.1) in space by the standard Galerkin method, as follows. For 0 < h < 1, let Sh be a family of finite-dimensional subspaces of H01 in which approximations to the solution u(·, t) of (1.1) will be sought for given t ∈ [0, T ]. We assume that Sh satisfies the approximation property that there exists an integer r ≥ 2 and a constant c > 0
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
3
independent of h such that for v ∈ H s ∩ H01 (1.2)
inf kv − ϕk + hkv − ϕk1 ≤ chs kvks ,
ϕ∈Sh
1≤s≤r
and the inverse property that for some c > 0 independent of h, kϕk1 ≤ ch−1 kϕk,
(1.3)
∀ϕ ∈ Sh .
As in (1.2), (1.3), in the sequel the symbols c, C, ci etc. will denote generic constants independent of the discretization parameter h and the time step. Such constants may also depend on the solution and the data of (1.1). We define now the semidiscrete approximation of the solution u(t) = u(·, t) of (1.1) in Sh in the customary way as the map uh : [0, T ] → Sh satisfying (with β(t) = β(·, t)): (1.4)
(uht, ϕ) = −iαa(uh , ϕ) + i(β(t)uh , ϕ),
uh (0) = u0h ,
where, for ϕ, χ ∈ H 1 , a(ϕ, χ) =
N X
∀ϕ ∈ Sh , 0 ≤ t ≤ T,
(∂j ϕ, ∂j χ).
j=1
In addition, we shall henceforth assume that u0 ∈ H r ∩ H01 and that u0h is an element of Sh such that (1.5)
ku0 − u0h k ≤ chr ku0 kr ;
for example u0h = P u0, where P : L2 → Sh is the L2 −projection operator onto Sh . If we introduce the linear operators ∆h : Sh → Sh and Bh (t) : Sh → Sh , Lh (t) : Sh → Sh , 0 ≤ t ≤ T, by (1.6) (1.7) (1.8)
(∆h ϕ, χ) = −a(ϕ, χ), (Bh (t)ϕ, χ) = (β(t)ϕ, χ),
∀ϕ, χ ∈ Sh ,
∀ϕ, χ ∈ Sh ,
Lh (t) = α∆h + Bh (t),
i.e.,
Bh (t)ϕ = P [β(·, t)ϕ],
0 ≤ t ≤ T,
and note that, since α and β are real, ∆h , Bh (t) and Lh (t), 0 ≤ t ≤ T, are Hermitian operators on {Sh , (·, ·)}, we may write (1.4) as (1.9)
uht = iLh (t)uh ,
0 ≤ t ≤ T,
uh (0) = u0h .
The equations (1.4) or (1.9) represent an initial-value problem for a system of ordinary differential equations that obviously has a unique solution uh (t) for 0 ≤ t ≤ T ; they will only be used in the sequel in order to motivate the fully discrete approximations. Let us only remark here that it is not hard to show, by comparing e.g. uh to the elliptic projection of u in the standard way, cf. [22], [4] that if u0h satisfies (1.5), then for t ∈ [0, T ] there holds Z t n o 0 r k(u − uh )(t)k ≤ ch ku kr + ku(s)kr + kut (s)kr ds , 0
4
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
i.e. that uh (t) satisfies an optimal-order L2 error estimate if u is smooth enough. We shall discretize (1.9) in time using the well-known class of implicit Runge–Kutta procedures of collocation type known as the q−stage (q ≥ 1 integer) Gauss–Legendre methods, [5], [6], [9]. The methods are defined by constants A = (aij ) ∈ Rq×q , b = (b1 , . . . , bq )T ∈ Rq , τ = (τ1 , . . . , τq )T ∈ Rq that are constructed in the standard way. Specifically, the τi are the —distinct, in (0, 1)— zeros of the shifted Legendre polynomials (d/dx)q xq (1 − x)q , the weights bi are defined as the solution of the q × q linear system of equations represented by (3.1.1) below for 1 ≤ ℓ ≤ q, while the coefficients aij are defined for each i, 1 ≤ i ≤ q, as solutions of the q × q linear system of equations represented by (3.1.2). Let k > 0 be the (constant) time step, let tn = nk, n = 0, 1, . . . , M, where T = Mk, and tn,i = tn + τi k, 1 ≤ i ≤ q. Then, the q−stage Gauss–Legendre methods applied to the system of ordinary differential equations represented by (1.9) yield the following full discretization of (1.1): For 0 ≤ n ≤ M, we seek U n ∈ Sh , approximating u(tn ), and U n,m ∈ Sh , 1 ≤ m ≤ q, that satisfy (a) U 0 = u0h ,
for n = 0, 1, . . . , M − 1 : (1.10)
(b) U
n,m
n
= U + ik
q X
n,j amj Ln,j , h U
j=1
(c) U
n+1
n
= U + ik
q X
1 ≤ m ≤ q,
n,j bj Ln,j , h U
j=1
Lh (tn ), Ln,j h
where Lnh = = Lh (tn,j ). In Section 2 we shall show that, for each n, the system (1.10b) has a unique solution {U n,m }, 1 ≤ m ≤ q, in (Sh )q and therefore that U n+1 is defined uniquely in Sh for 0 ≤ n ≤ M − 1, by (1.10c). We shall also verify that (1.10) is unconditionally stable in L2 and, in particular, conservative, i.e. that it satisfies kU n k = kU 0 k, 0 ≤ n ≤ M; thus it mimicks the behavior of (1.1) and (1.9) for which ku(t)k = ku0k and kuh (t)k = ku0h k hold, respectively, for 0 ≤ t ≤ T. In Section 3 we shall estimate the error of the approximation U n in the L2 norm; specifically, we shall show in Theorem 3.1, which is the main result of the paper, that (1.11) max kU n − u(tn )k ≤ c k min(2q,q+2) + hr . 0≤n≤M
This error bound is of optimal order in space. As far as the temporal rate of convergence is concerned, it is well-known that the q−stage Gauss–Legendre methods have (classical) order 2q when applied to nonstiff ordinary differential equations. Therefore our proof certainly gives optimal-order temporal convergence, resp. 2, 4, for the one- and two-stage, resp., methods that seem to enjoy current practical importance. For q > 2 our result — O(k q+2) in time — shows the effect of “reduction of order due to stiffness”. This result is no worse than analogous estimates proven in the literature for Runge– Kutta full discretizations of initial- and boundary-value problems for p.d.e.’s with timedependent coefficients or nonlinear terms, posed with Dirichlet boundary conditions: In
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
5
his thesis Broc´ehn, [4], considers full discretizations of the Schr¨odinger equation with a general second-order elliptic operator with time-dependent coefficients using some semi-implicit Runge–Kutta methods. (The class of Gauss–Legendre schemes considered here is not semi-implicit with the exception of the one-stage scheme.) For his schemes, using different estimation techniques, he proves error bounds with temporal order of accuracy equal to min(p, q + 1), where, in his notation, q is the order of the quadrature rule associated with the intermediate stages of the Runge–Kutta method and plays an analogous role to the q used here, and p is the classical order (equal to 2q for the Gauss–Legendre methods). For the special elliptic operator of the right-hand side of (1.1) he obtains a temporal rate of q + 2, if p = q − 2, and under the mesh condition that kh−2 remain bounded as k, h → 0. In Theorem 2 of [8], the temporal discretization, by a class of Runge–Kutta schemes (disjoint from the Gauss–Legendre methods but suitable for parabolic problems) of an abstract semilinear parabolic equation is shown to have a rate of convergence exhibiting an analogous limitation to our q + 2 result. In [3] the Gauss–Legendre methods are applied to a nonlinear p.d.e., the generalized Korteweg–de Vries equation, posed in one space dimension with periodic boundary conditions and discretized in space with smooth periodic splines on a uniform mesh. The exact analog to (1.11) is proved then by a different technique from the one at hand, with the details of the space discretization and the periodicity of the exact and discrete solutions playing a crucial role. (After the completion of the original version of the paper at hand, we learnt that Karakashian and McKinney, [12], proved the optimal order 2q for the Korteweg–de Vries equation; their remarkable proof again relies heavily on the periodic boundary conditions.) Finally, in Section 4 we confine attention to the 2−stage Gauss–Legendre method and devise a scheme that avoids solving the 2 dim Sh × 2 dim Sh linear system represented by (1.10b) for q = 2. A suitable decoupling strategy and an iteration scheme enables us to produce stable and optimal-order accurate approximations to U n by solving a number of linear systems of size dim Sh × dim Sh at each time step; these systems will have sparse matrices if Sh is furnished with a finite element basis with elements of small support. Acknowledgement. The authors record their thanks to a referee for pointing out reference [4] to them.
2. Existence and stability of the fully discrete approximations In this section we shall show that for each n the linear system represented by (1.10b) has a unique solution U n,m , 1 ≤ m ≤ q, and that the resulting overall fully discrete scheme (1.10) is stable (conservative) in the L2 norm. For this purpose we shall make use of the following well-known properties of the Gauss–Legendre methods, [6], [9]:
6
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
(2.1) For each q ≥ 1 there exists a diagonal q × q matrix D with positive
diagonal elements, such that the matrix DAD −1 is positive definite on Rq . (See e.g. [9, Theorem 5.5.6, Cor. 5.1.4 and (5.1.23)].)
(2.2)
bi > 0, bi aij + bj aji − bi bj = 0,
1 ≤ i, j ≤ q,
i.e. the Gauss–Legendre methods are conservative in the nonlinear context; cf. e.g. [7], [9, p.117]. We first examine the existence of solutions U n,m of the linear system (1.10b). On the product space (Sh )q we let Un denote the vector (U n,1 , . . . , U n,q )T and Lnh : (Sh )q → n,q T (Sh )q be the diagonal operator given, for V ∈ (Sh )q , by Lnh V = (Ln,1 h V1 , . . . , Lh Vq ) . We write then the equations (1.10a,b) respectively as (2.3a)
Un = U n e + ikALnh Un ,
(2.3b)
U n+1 = U n + ikbT Lnh Un .
In (2.3a,b) and in the sequel, abusing notation a bit to avoid tensor products, for V ∈ P P (Sh )q we let bT V = qi=1 bi Vi , AV ∈ (Sh )q : (AV)i = qj=1 aij Vj , e = (1, 1, . . . , 1)T ∈ Rq and for U ∈ Sh , Ue = (U, U, . . . , U)T ∈ (Sh )q . The existence of Un , solution of (2.3a), will follow from the following general lemma. (In the sequel we let ((·, ·)), resp. ||| · |||, denote the product inner product, resp. norm, in (L2 )q .) Lemma 2.1. Suppose that V = (Vi ) and W = (Wi ) in (S h )q satisfy the equation (2.4)
V = W + ikAF(V),
where F : (Sh )q → (Sh )q is a diagonal mapping such that (F(V ))i = Fi (Vi ), 1 ≤ i ≤ q, where Fi : Sh → Sh are given mappings with the property that Im(Fi (ϕ), ϕ) = 0, for ϕ ∈ Sh , 1 ≤ i ≤ q. Then, there exists a constant c, depending only on the constants of the Gauss–Legendre method, such that (2.5)
|||V||| ≤ c|||W|||.
Proof. Let D = diag (d1 , . . . , dq ), di > 0, be the diagonal matrix mentioned in (2.1). Multiplying (2.4) by D 2 A−1 on the left and taking the (L2 )q inner product with V we obtain (2.6)
((D 2 A−1 V, V)) = ((D 2 A−1 W, V)) + ik((D 2 F(V), V)).
By our hypotheses Re[ik((D 2 F(V), V))] = 0. Hence, taking real parts in (2.6) gives (2.7)
Re((D 2 A−1 V, V)) = Re((D 2 A−1 W, V)).
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
7
Denote DA−1 D −1 =: B = (bij ) ∈ Rq×q . By (2.1) there exists λ > 0 such that Pq Pq 2 q i=1 ξi for every ξ ∈ R . Hence, putting DV = Y we have i,j=1 bij ξi ξj ≥ λ Re((D 2 A−1 V, V)) = Re((DA−1 D −1 Y, Y)) = Re((BY, Y)) ≥ λ ||| Re Y|||2 + ||| Im Y|||2 = λ|||Y|||2 ≥ λ min d2i |||V|||2 . i
(2.5) follows then by (2.7) and the Cauchy–Schwarz inequality.
Given now U n ∈ Sh apply Lemma 2.1 to the linear system (2.3a) where F(V) = n,i Lnh V, Fi (Vi ) = Ln,i h Vi and Im(Lh ϕ, ϕ) = 0 for ϕ ∈ Sh . By (2.5) we see that the homogeneous system (U n = 0 in (2.3a)) has only the trivial solution. Hence, given U n ∈ Sh , (2.3a) has a unique solution Un = (U n,i ) in (Sh )q , which satisfies max kU n,i k ≤ ckU n k,
(2.8)
i
for some constant c that depends only on the Gauss–Legendre method. The stability (conservativeness) of the scheme (1.10) follows from the following result, stated again in slightly more general terms: Lemma 2.2. Suppose, given U ∈ Sh , that V ∈ (Sh )q and Y ∈ Sh satisfy the equations (2.9a)
V = Ue + ikAF(V),
(2.9b)
Y = U + ikbT F(V),
where F is a mapping that satisfies the hypotheses of Lemma 2.1. Then (2.10)
kY k = kUk.
Proof. (2.9b) gives kY k2 = kUk2 + ik(bT F(V), U) − ik(U, bT F(V)) + k 2 (bT F(V), bT F(V)) 2
= kUk − 2k Im
q X
bj (Wj , U) + k
j=1
2
q X
bj bℓ (Wj , Wℓ ),
j,ℓ=1
where Wj = Fj (Vj ). Let W = (W1 , . . . , Wq )T . Replacing U in the right-hand side of the above equation by its expression U = Vj − ik(AW)j that (2.9a) gives, we see, using the properties of F and (2.2) that 2
2
kY k = kUk − 2k 2
= kUk − k
2
2
q X
bi Re(Wi , (AW)i ) + k
q X
bi bj (Wi , Wj )
i,j=1
i=1
q hX
2
i
(bi aij + bj aji − bi bj )(Wi , Wj ) = kUk2 .
i,j=1
Applying this result to the scheme (1.10) we obtain (2.11)
kU n k = kU 0 k,
0 ≤ n ≤ M,
i.e. that our fully discrete method is conservative in the L2 sense.
8
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
3. Consistency and convergence In this section we shall study the consistency of the fully discrete scheme (1.10) and prove the error estimate (1.11). To this effect we shall first list some well-known algebraic properties of the Gauss–Legendre methods, [6], [9], that will be used in the sequel along with (2.1) and (2.2). (3.1) The q−stage Gauss–Legendre method is consistent of order 2q (i.e. has accuracy of order 2q when applied to an ordinary differential equation y ′ = f (t, y), where f and its partial derivatives of sufficiently high order with respect to y and t are smooth and bounded), and satisfies the following order (simplifying) conditions, [5]: q X
(3.1.1)
bi τiℓ−1 = ℓ−1 ,
i=1
q X
(3.1.2)
j=1
(3.1.3)
q X i=1
(3.2)
aij τjℓ−1 = τiℓ /ℓ,
1 ≤ ℓ ≤ 2q,
1 ≤ i, ℓ ≤ q,
bi τiℓ−1 aij = ℓ−1 bj (1 − τjℓ ),
1 ≤ j, ℓ ≤ q.
The q−stage Gauss–Legendre method corresponds to the q − th diagonal
Pad´e rational approximation to the exponential, i.e. if r(z) = 1 + zbT (I − zA)−1 e, e = (1, . . . , 1)T ∈ Rq ,
then r(z) is the q−th diagonal Pad´e approximant to exp(z). For the purposes of the proof of convergence we shall compare the solution U n of (1.10) to the elliptic projection W = W (t) ∈ Sh , 0 ≤ t ≤ T, of the solution u(t) of (1.1), defined as usual by (3.3)
(∆h W, ϕ) = −a(W, ϕ) = (∆u, ϕ) ∀ϕ ∈ Sh .
We shall denote the associated (time-independent) elliptic projection operator onto Sh (defined on H 2 ∩H01 ) by PI . In this notation, W (t) = PI u(t) and obviously W (j) (t) = PI u(j) (t); here and in the sequel v (j) (t) = (d/dt)j v(t). By our assumptions on Sh there follows that (3.4)
kv − PI vk + hkv − PI vk1 ≤ chr kvkr ,
∀v ∈ H r ∩ H01 .
Obviously (3.3) implies that kW (j) (t)k1 ≤ cku(j) (t)||1 ≤ cj , j ≥ 0. In addition, there exist constants cj such that kLh (t)W (j) (s)k ≤ cj , t, s ∈ [0, T ], j ≥ 0; this follows from
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
9
the observation that for any ϕ ∈ Sh Lh (t)W (j) (s), ϕ = α(∆h W (j) (s), ϕ) + (β(t)W (j) (s), ϕ) ≤ |α| k∆u(j)(s)k + |β(t)|∞ku(j) (s)k1 kϕk ≤ cku(j) (s)k2 kϕk.
(j)
In fact, if Lh (t) : Sh → Sh denotes the j−th time derivative of the operator Lh (t), (j) (j) (j) given by (Lh (t)ϕ, χ) = (∂tj β(·, t)ϕ, χ), j ≥ 1, for χ, ϕ ∈ Sh , i.e. by Lh = Bh , j ≥ 1, (i) we have kLh (t)W (j) (s)k ≤ ci ku(j)(s)k1 , i ≥ 1, j ≥ 0. Thus, we can generalize the previous estimate to (i)
kLh (t)W (j) (s)k ≤ cij ,
(3.5)
i, j ≥ 0,
t, s ∈ [0, T ],
and also note that (i)
(3.6)
kLh (t)k ≤ ci ,
i ≥ 1,
where k · k denotes here the L2 induced operator norm on Sh . We shall also make use of the following property of the elliptic projection, namely that for constants cij kLh (t)B (j) (s)W (i) (s)k ≤ cij ,
(3.7)
i, j ≥ 0,
t, s ∈ [0, T ],
which may be proved as follows. We have (3.8)
(j)
(j)
(j)
Lh (t)Bh (s)W (i) (s) = α∆h Bh (s)W (i) (s) + Bh (t)Bh (s)W (i) (s). (j)
(j)
Since for j ≥ 0, Bh (t)ϕ = P (∂tj β(t)ϕ), we have kBh ϕk ≤ cj kϕk for ϕ ∈ Sh . Hence in (3.8) (j)
kBh (t)Bh (s)W (i) (s)k ≤ cij .
(3.9)
Also for ϕ ∈ Sh we obtain, suppressing the dependence on s, (j)
(3.10)
(j)
−(∆h (t)Bh W (i) , ϕ) = a(Bh W (i) , ϕ) = a(P [β (j) W (i) ], ϕ)
= a(P [β (j)W (i) ] − β (j) u(i) , ϕ) + a(β (j) u(i) , ϕ).
We obviously have (3.11)
|a(β (j) u(i) , ϕ)| = |(∆(β (j) u(i) ), ϕ)| ≤ cij kϕk,
and by (1.3) (3.12)
|a(P [β (j) W (i) ] − β (j) u(i) , ϕ)| ≤ kP [β (j)W (i) ] − β (j) u(i) k1 kϕk1
≤ ch−1 kP [β (j) W (i) ] − β (j) u(i) k1 kϕk.
Now (3.13) kP [β (j)W (i) ] − β (j) u(i) k1 ≤ kP [β (j) (W (i) − u(i) )]k1 + kP [β (j)u(i) ] − β (j) u(i) k1 .
10
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
Since, as it may easily be seen from (1.3), (3.4), kP v −vk1 ≤ chr−1 kvkr for v ∈ H r ∩H01 , and since for β and u sufficiently smooth β (j) u(i) ∈ H r ∩ H01, there follows kP [β (j)(W (i) − u(i) )]k1 + kP [β (j) u(i) ] − β (j) u(i) k1
≤ ch−1 kP [β (j) (W (i) − u(i) )]k + cij hr−1 ≤ cij hr−1 .
These estimates, when substituted in (3.13), yield, in conjunction with (3.10)–(3.12), (j) that k∆h Bh W (i) k ≤ cij since r ≥ 2. Then (3.7) follows from (3.8) and (3.9). We now embark upon the proof of the main result of this section. For a function v defined on [0, T ] we generally denote v n = v(tn ). We first define, for the purposes of the proof of consistency, V n,m for 0 ≤ n ≤ M − 1, 1 ≤ m ≤ q and V n , 0 ≤ n ≤ M in Sh by (a) V 0 = W 0 , (b) V n,m = W n + ik (3.14)
q X
n,j amj Ln,j , h V
j=1
(c) V
n+1
n
= W + ik
q X
1 ≤ m ≤ q,
n,j bj Ln,j . h V
j=1
In Proposition 3.1 below we shall prove the consistency result, valid for u sufficiently smooth: (3.15) max kV n − W n k ≤ ck k min(2q,q+2) + hr . 0≤n≤M
If this holds, then a simple stability calculation, as the following theorem shows, gives the error bound (1.11): Theorem 3.1. Let u be sufficiently smooth and suppose that (3.15) holds. Then (3.16) max kU n − un k ≤ c k min(2q,q+2) + hr . 0≤n≤M
Proof. Let V n,m , V n be defined by (3.14) and let εn,m = U n,m −V n,m , εn = U n −V n , ζ n = U n − W n . Then (1.10) and (3.14) yield ε
n,m
n
= ζ + ik
q X
n,j amj Ln,j h ε ,
j=1
ε
n+1
n
= ζ + ik
q X
1 ≤ m ≤ q,
n,j bj Ln,j h ε .
j=1
The stability Lemma 2.2 gives then that kεn+1k = kζ n k. Hence, for 0 ≤ n ≤ M − 1, kζ n+1k ≤ kεn+1 k + kV n+1 − W n+1 k = kζ n k + kV n+1 − W n+1 k.
Therefore, by (3.15) kζ n k ≤ kζ 0 k + c(k min(2q,q+2) + hr ), 0 ≤ n ≤ M, and the result follows from (3.4), (1.5) and the triangle inequality. Hence our task is to prove consistency:
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
11
Proposition 3.1. If u is sufficiently smooth, (3.15) holds. Proof. We follow, up to a point, the technique of the consistency proof for Runge– Kutta discretizations of partial differential equations introduced in [11]. First define τij , 1 ≤ i ≤ q, j ≥ 0 by τi0 = 1, τij =
(3.17)
q X
aiℓ τℓ,j−1,
ℓ=1 j
⇐⇒ τij = (A e)i ,
1 ≤ i ≤ q, j ≥ 0,
1 ≤ i ≤ q, j ≥ 0.
Note that by (3.1.2) we may infer that τij = (τi )j /j!,
(3.18)
1 ≤ i ≤ q, 0 ≤ j ≤ q.
Also define, for 1 ≤ m ≤ q, 0 ≤ n ≤ M − 1 n
(3.19)
Λm W =
2q X
τmj k j W (j)n ,
j=0
en,m = V n,m − Λm W n .
(3.20)
We now make a preliminary observation. By (3.14) and (3.20) we have (3.21)
V
n+1
n
=W +
q X
m,j=1
bm (A−1 )mj (Λj W n − W n ) + bT A−1 en ,
where en = (en,1 , . . . , en,q )T ∈ (Sh )q . Using (3.19) and (3.17) we can write q X
m,j=1
bm (A−1 )mj (Λj W n − W n ) =
q X
bm (A−1 )mj
m,j=1
2q hX
(Aℓ e)j k ℓ W (ℓ)n
ℓ=1
i
2q X X T ℓ−1 ℓ (ℓ)n k ℓ W (ℓ)n /ℓ! (b A e)k W = = 2q
ℓ=1
ℓ=1
where in the last equality we used the identities (3.22)
bT Aℓ−1 e = 1/ℓ!,
1 ≤ ℓ ≤ 2q,
that follow from the fact that the rational approximation r(z), cf. (3.2), corresponding to the q−stage Gauss–Legendre method is an O(z 2q+1 ) approximation to exp(z) as z → 0. We conclude therefore by (3.21) and Taylor’s theorem that kV n+1 − W n+1 k ≤ ck 2q+1 + kbT A−1 en k. Hence, in order to prove (3.15) our preliminary observation is that it is sufficient to obtain an estimate of the form (3.23) kbT A−1 en k ≤ ck k min(2q,q+2) + hr .
12
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
Note that (3.20) when substituted in (3.14b) yields for 0 ≤ n ≤ M − 1 n,j
(3.24)
e
=E
n,j
+ ik
q X
n,ℓ ajℓ Ln,ℓ h e ,
ℓ=1
where we have put (3.25)
E
n,j
n
n
= −Λj W + W + ik
q X
1 ≤ j ≤ q,
n ajℓ Ln,ℓ h Λℓ W ,
ℓ=1
1 ≤ j ≤ q.
The first main step of the proof consists of some long intermediate calculations, inevitable in any order estimation of Runge–Kutta methods, aiming towards transforming the right-hand side of (3.25) to a form more suitable for our purposes. To simplify notation a bit, in the sequel by ϕ = O(k λ + hµ ) we shall mean that there exists a constant c > 0 independent of k and h such that kϕk ≤ c(k λ + hµ ) for k, h sufficiently small. First, using (3.25), (3.5), (3.17), (3.19) we have for 1 ≤ j ≤ q E
n,j
n
n
= −Λj W + W + ik + ik
q X
q X
ajℓ Ln,ℓ h
2q−1 X
τℓm k W
(m)n
2q−1 X
τℓm k m W (m)n
m
m=0
ℓ=1
(2q)n ajℓ τℓ,2q k 2q Ln,ℓ h W
ℓ=1
n
n
= −Λj W + W + ik q
+ ik
X
ajℓ (Ln,ℓ h
ℓ=1
(3.26)
−
q X
=−
Lnh )
m
τjm k W
(m)n
+
ajℓ (Ln,ℓ h
ℓ=1
−
Lnh )
m=1
+ ik
ℓ=1
=
I1n,j
+
2q−1 X
ℓ=1
n ajℓ (Ln,ℓ h − Lh )
I2n,j
+ O(k
2q+1
2q−1 X
),
m=0
+ O(k 2q+1 )
ajℓ τℓ,m−1 ik m Lnh W (m−1)n m
τℓm k W
τjm k m W (m)n − iLnh W (m−1)n
q X
(m)n
(m)n
m=0
2q
X
τℓm k W
q XX
m=1
X
m
m=0
q
=−
2q−1 X 2q
m=1
+ ik
m=0
ℓ=1
2q
X
ajℓ Ln,ℓ h
+ O(k 2q+1 )
τℓm k m W (m)n + O(k 2q+1 )
where (3.27)
I1n,j := −
2q X
m=1
τjm k m W (m)n − iLnh W (m−1)n ,
1 ≤ j ≤ q,
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
I2n,j := i
(3.28)
q X ℓ=1
In
I1n,j
n ajℓ (Ln,ℓ h − Lh )
2q−1 X m=0
τℓm k m+1 W (m)n ,
13
1 ≤ j ≤ q.
we have, using (1.8)
W (m)n − iLnh W (m−1)n = ∂tm−1 Wtn − iα∆h W n − iBhn W (m−1)n . Introducing now for 0 ≤ t ≤ T ψ(t) := Wt − ut − iβ(t) W (t) − u(t) , for which (3.4) shows that ψ (j) (t) = O(hr ), we can write using (3.3) and (1.1) that (3.29)
(3.30)
Wt − iα∆h W = P ψ + iBh (t)W.
Differentiating (3.30) with respect to t and using (3.29) we rewrite (3.27) as I1n,j (3.31)
2q X
= −i
m=1 2q X
τjm k m ∂tm−1 (Bhn W n ) − Bhn W (m−1)n + O(khr )
i h m−2 X m − 1 (m−1−ℓ)n Bh W (ℓ)n + O(khr ). τjm k = −i ℓ m=2 ℓ=0 m
Turning now to the term I2n,j we obtain by (3.6) and Taylor’s theorem (3.32)
I2n,j
=i
q X
2q 2q−1 X X µ (µ)n m (m−1)n ajℓ (τℓ k) Bh /µ! + O(k 2q+1 ). τℓ,m−1 k W µ=1
ℓ=1
m=1
Use of (3.5) yields for 1 ≤ ℓ ≤ q that 2q−1 2q
Zℓn
XX
:=
(j)n
(τℓ )j τℓ,m−1 k j+m Bh W (m−1)n /j!
j=1 m=1
2q
=
X
k
λ
λ−1 hX
(λ−m)n
(τℓ )λ−m τℓ,m−1 Bh
m=1
λ=2
i W (m−1)n /(λ − m)! + O(k 2q+1 ).
Therefore (3.32) gives for 1 ≤ j ≤ q (3.33)
I2n,j
=i
q X
ajℓ
ℓ=1
+ O(k
2q nX
k
λ
λ−1 hX
(λ−m)n
(τℓ )λ−m τℓ,m−1 Bh
m=1
λ=2 2q+1
).
io W (m−1)n /(λ − m)!
Summarizing, we obtain by (3.26), (3.31) and (3.33) for 1 ≤ j ≤ q (3.34)
E
n,j
=i
2q X
k
λ=2
+ O(k
λ
λ−1 hX
(λ−m)n
(δjm,λ − γjm,λ )Bh
m=1 2q+1
+ khr ),
i W (m−1)n /(λ − m)!
where, for 1 ≤ j ≤ q, λ = 2, . . . , 2q, m = 1, . . . , λ − 1: (3.35)
δjm,λ
=
q X ℓ=1
ajℓ (τℓ )λ−m τℓ,m−1 ,
γjm,λ = τjλ (λ − 1)!/(m − 1)!.
14
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
We immediately observe, using (3.18), that for λ = 1, . . . , q + 1, 1 ≤ m ≤ λ − 1 δjm,λ
=
q X
ajℓ (τℓ )
λ−1
ℓ=1
q X /(m − 1)! = (λ − 1)!/(m − 1)! ajℓ τℓ,λ−1 = γjm,λ . ℓ=1
Therefore, we finally conclude that if q ≥ 2, (3.34) yields for 1 ≤ j ≤ q (3.36)
E
n,j
=i
2q X
kλ
m=1
λ=q+2
+ O(k
λ−1 hX
2q+1
(λ−m)n
(δjm,λ − γjm,λ )Bh
+ khr ),
i W (m−1)n /(λ − m)!
while, simply, if q = 1 E n,1 = O k(k 2 + hr ) .
(3.37)
Note that (3.37) used in conjunction with
max ken,j k ≤ c max kE n,j k,
(3.38)
1≤j≤q
1≤j≤q
(which follows from the stability estimate of Lemma 2.1 applied to the equations (3.24)) gives the desired estimate (3.23) in the case q = 1. Therefore we henceforth concentrate on the cases q ≥ 2 for which E n,j is given by (3.36). To this end let for 1 ≤ j ≤ q, (3.39)
ϕ
n,j
=i
2q X
k
λ=q+2
λ
λ−1 hX
(δjm,λ
m=1
−
(λ−m)n γjm,λ )Bh W (m−1)n /(λ
i − m)! .
Then (3.36) is written as (3.40)
E n,j = ϕn,j + O(k 2q+1 + khr ),
1 ≤ j ≤ q.
Obviously, (3.5) gives that ϕn,j = O(k q+2) since δjm,λ − γjm,λ are not zero in general if λ ≥ q + 2. (3.40) then gives that E n,j = O(k q+2 + khr ) and (3.38) implies in turn that en,j = O(k q+2 + khr ) thus yielding the estimate bT A−1 en = O(k(k q+1 + hr )). (This concludes essentially the application of the idea of the consistency proof of [11] to the case at hand). The second step in the proof is the improvement of the q + 1 exponent of k to the better value q + 2. For this purpose we use the fact that we must actually estimate not the individual en,j but their particular linear combination bT A−1 en =
q X
bi (A−1 )ij en,j .
i,j=1
First note that defining (3.41)
E˜ n,j = E n,j − ϕn,j
e˜n,j = en,j − ϕn,j ,
we may rewrite (3.24) as q q X X n,ℓ n,ℓ n,j n,j ajℓ Ln,ℓ e˜n,ℓ . ajℓ L ϕ + ik (3.42) e˜ = E˜ + ik h
h
ℓ=1
ℓ=1
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
15
Recall from (3.40), (3.41) that E˜ n,j = O(k 2q+1 + khr ). Also, (3.7) and (3.39) yield that n,j q+2 kLn,j , 1 ≤ j ≤ q. Hence, the stability estimate of Lemma 2.1, applied to h ϕ k ≤ ck n,j (3.42), yields e˜ = O(k q+3 + khr ). Therefore, it follows from (3.41) that (3.43)
bT A−1 en = O(k(k q+2 + hr )) + bT A−1 ϕn ,
where ϕn = (ϕn,1 , . . . , ϕn,q )T ∈ (Sh )q . Hence, the desired estimate (3.23) will follow from the fact that actually bT A−1 ϕn = 0,
(3.44)
which is a consequence of (3.39) and the following cancellation property of the Gauss– Legendre methods that we state as a separate lemma: Lemma 3.1. Let δjm,λ , γjm,λ be defined for 1 ≤ j ≤ q, λ = 2, . . . , 2q, m = 1, . . . , λ − 1 by T T (3.35) and denote δ m,λ = δ1m,λ , . . . , δqm,λ and γ m,λ = γ1m,λ , . . . , γqm,λ ∈ Rq . Then
(3.45)
bT A−1 (δ m,λ − γ m,λ ) = 0.
Proof. For λ ≤ q + 1, 1 ≤ m ≤ λ − 1, we have already established by the simple calculation following (3.35) that δ m,λ = γ m,λ . Hence we restrict attention to the interesting case q + 2 ≤ λ ≤ 2q, 1 ≤ m ≤ λ − 1. Define for the purposes of this lemma T ∈ Rq×q as T = diag (τ1 , . . . , τq ) — no confusion with the T of (1.1) will arise — and note that (3.1.1), (3.1.2) can be written equivalently as (3.1.1′ )
bT T ℓ−1 e = 1/ℓ,
(3.1.2′ )
AT ℓ−1 e = T ℓ e/ℓ,
1 ≤ ℓ ≤ 2q, 1 ≤ ℓ ≤ q,
respectively; (3.1.2′ ) implies that (3.1.2′′ )
Aℓ e = T ℓ e/ℓ!,
0 ≤ ℓ ≤ q.
Now, (3.17) gives that δ m,λ = AT λ−m Am−1 e. On the other hand, using (3.22), since λ ≤ 2q, we have bT A−1 γ m,λ = (λ − 1)!bT Aλ−1 e/(m − 1)! = (λ(m − 1)!)−1 . Hence, to show (3.45) it suffices to establish (3.46)
bT T λ−m Am−1 e = (λ(m − 1)!)−1 ,
q + 2 ≤ λ ≤ 2q, 1 ≤ m ≤ λ − 1.
Obviously for each λ, q + 2 ≤ λ ≤ 2q, (3.46) holds for 1 ≤ m ≤ q + 1; this follows from (3.1.2′′ ) and (3.1.1′) that yield bT T λ−m Am−1 e = bT T λ−1 e/(m − 1)! = (λ(m − 1)!)−1 . Hence we focus on the case q + 2 ≤ λ ≤ 2q, q + 2 ≤ m ≤ λ − 1. Using (3.1.2′′ ) again gives, since m − 1 ≥ q + 1, (3.47)
bT T λ−m Am−1 e = bT T λ−m Am−1−q T q e/q!.
For integers k ≥ 1, 1 ≤ ℓ ≤ q define (3.48)
F (ℓ, k) = bT T ℓ−1 Ak T q e
16
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
and note that by (3.1.3) T
F (ℓ, k) = b T
ℓ−1
k−1
A(A
q
T e) =
q X
bi τiℓ−1
= ℓ−1
(3.49)
X j=1
=ℓ
−1
aij (Ak−1 T q e)j
j=1
i=1
q
q X
bj (1 − τjℓ )(Ak−1 T q e)j
bT Ak−1 T q e − bT T ℓ Ak−1 T q e
= ℓ−1 F (1, k − 1) − F (ℓ + 1, k − 1) ,
k ≥ 1, 1 ≤ ℓ ≤ q.
We can now calculate directly using (3.1.1), (3.1.3) for integer s, such that 1 ≤ s ≤ q−1, that F (s, 1) =
q X
bi τis−1
= s−1
aij τjq
q
hX j=1
bj τjq −
=s
−1
q X j=1
j=1
i=1
(3.50)
q X
q
X j=1
bj (1 − τjs )τjq
i −1 bj τjs+q = (q + 1)(s + q + 1) .
We claim now that (3.51)
−1 F (ℓ, k) = (q + 1) · · · (q + k)(q + ℓ + k)
for k ≥ 1, ℓ ≥ 1, ℓ + k ≤ q.
Indeed (3.50) shows that (3.51) holds for k = 1. For the inductive step assume that (3.51) holds for all k ≤ k ′ and 1 ≤ ℓ ≤ q − k. Then, for ℓ ≤ q − k ′ − 1 we have by (3.49) and the inductive hypothesis that F (ℓ, k ′ + 1) = ℓ−1 F (1, k ′) − F (ℓ + 1, k ′ ) −1 −1 = (q + 1) · · · (q + k ′ )(q + k ′ + 1) − (q + 1) · · · (q + k ′ )(q + k ′ + ℓ + 1) /ℓ −1 = (q + 1) · · · (q + k ′ )(q + k ′ + 1)(q + ℓ + k ′ + 1) . This completes the inductive step and shows the validity of (3.51). Using now (3.48), (3.51) and (3.47) we obtain with k := m − 1 − q ≥ 1, ℓ := λ − m + 1 ≥ 1, since k + ℓ = λ − q ≤ q in the region of interest, that for q + 2 ≤ λ ≤ 2q, q + 2 ≤ m ≤ λ − 1, −1 −1 bT T λ−m Am−1 e = q!(q + 1) · · · (m − 1)λ = λ(m − 1)! . This identity completes the proof of the validity of (3.46).
As a consequence of (3.45) and (3.39), (3.44) holds and the proof of Theorem 3.1 is now complete. 4. Practical implementation of the two-stage scheme In this section we shall study questions related to the efficient implementation of the fully discrete scheme generated by the two-stage Gauss–Legendre method which is
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
17
given of course by the tableau, [9], [6]:
(4.1)
A τ T
b
1 4 1 4
=
1 4
1 √ 2 3
1 4
1 √ 2 3
+
−
1 2
1 2 1 2
− +
1 √ 2 3 1 √ 2 3
.
1 2
With these values of aij , bi , τi the method is (a) U 0 = u0h , for n = 0, 1, . . . , M − 1 : (4.2)
(b) U
n,m
n
= U + ik
2 X
n,j amj Ln,j , h U
m = 1, 2,
j=1
(c) U
n+1
n
= U + ik
2 X
n,j bj Ln,j . h U
j=1
n,j Note that eliminating the Ln,j from (4.2b and c) and using the particular values h U of the constants aij and bi from (4.1), we may write (4.2c) simply as √ (4.2c′ ) U n+1 = U n + 3(U n,2 − U n,1 ).
In sections 2 and 3 we proved that for each 0 ≤ n ≤ M (4.2) has a unique solution, that the resulting scheme is L2 −conservative and that it satisfies the optimal in space and time estimate (4.3) max kU n − un k ≤ c k 4 + hr . 0≤n≤M
Let d = dim Sh . To determine U n,1 , U n,2 from (4.2b) one should solve, after choosing a basis for Sh , the 2d × 2d linear system J n Un = F n ,
(4.4) where (4.5)
Jn = Jn (tn,1 , tn,2 ) =
! n,2 I − ika11 Ln,1 −ika L 12 h h n,2 , −ika21 Ln,1 I − ika L 22 h h
Un = (U n,1 , U n,2 )T ,
Fn = (U n , U n )T .
With the aim of solving only (sparse) d × d systems of linear equations at each time step, we shall uncouple the two equations in (4.4) borrowing an idea from [3]. We first write (4.4) equivalently as (4.6)
J∗n Un = (J∗n − Jn )Un + Fn ,
where (4.7)
J∗n = Jn (t∗n , t∗n ),
t∗n = (tn,1 + tn,2 )/2 = tn + k/2,
18
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
the advantage being now that the operators Lh in the entries of J∗n are evaluated at the same point t∗n and that, in particular, they commute. This enables us to compute the solution Y = (Y1 , Y2)T ∈ (Sh )2 of systems of the form J∗n Y = Z,
(4.8)
given Z = (Z1 , Z2 )T ∈ (Sh )2 , “explicitly” as
∆∗n Y1 = Z1 + ikL∗n (a12 Z2 − a22 Z1 ),
(4.9)
∆∗n Y2 = Z2 + ikL∗n (a21 Z1 − a11 Z2 ),
where ∆∗n = I − k(iL∗n )/2 + k 2 (iL∗n )2 /12, L∗n = Lh (t∗n ). It is easily seen that ∆∗n is invertible since iL∗n is normal with purely imaginary eigenvalues and the polynomial 1 − z/2 + z 2 /12 has no roots on the imaginary axis. Systems with operator ∆∗n like the ones in (4.9) can then in turn be solved by the complex analog of the procedure proposed for Pad´e diagonal methods in [2]. Consider e.g. the first equation of (4.9) and, putting W1 = a12 Z2 − a22 Z1 , R1 = W1 + Z1 /2 − ikL∗n Z1 /12, Φ1 = Y1 − Z1 , rewrite it as ∆∗n Φ1 = ikL∗n R1 .
(4.10)
√ Since ∆∗n = (I − iµkL∗n )(I − i¯ µkL∗n ), where µ = 1/4 − i 3/12, we may rewrite (4.10), letting H n = I − iµkL∗n , K n = I − i¯ µkL∗n , as H n K n Φ1 = i(H n − K n )R1 /(2 Im µ), from which Φ1 = i[(K n )−1 − (H n )−1 ]R1/(2 Im µ).
(4.11)
To determine therefore Φ1 (and hence Y1 ) one must form R1 in the right-hand side of (4.11) and solve two complex linear systems of size d × d with operators K n and H n noting that at the level of matrix-vector operations the corresponding matrices will be sparse if a finite element basis consisting of functions of small support is chosen for Sh . Similarly for Φ2 with obvious parallelicity duly noted. In order to take advantage of the fact that systems of the form (4.8) can be “easily” solved in the manner outlined above, we shall solve the original system (4.6) by the simplest iterative method that its form immediately suggests. Let jn , n = 0, 1, . . . , M, be given positive integers — in practice jn = 1 or jn = 2 — representing the number of iterations that will be performed at step n to solve (4.6). For jn , n = 0, 1, . . . , M we shall compute Ujn approximating U n as follows: Let Uj00 = U 0 = u0h (e.g. take j0 = 0). For n = 0, 1, . . . , M: Compute suitable starting values U0n,1 , U0n,2 . Compute Ujn,m , m = 1, 2 by the iterative scheme: n+1 For j = 0, 1, . . . , jn+1 − 1 solve (4.12)
J∗n
n,1 Uj+1
Ujn,2 +1
!
= (J∗n − Jn )
Ujn,1 Ujn,2
!
+
Ujnn Ujn,2 n
!
.
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
19
√
− Ujn,1 ). Note that when j = jn+1 − 1 in the inner Define Ujn+1 = Ujnn + 3(Ujn,2 n+1 n+1 n+1 (j) loop in (4.12) there are important savings in the number of operations: since only the difference Ujn,2 − Ujn,1 is finally needed, we must solve directly only one system of n+1 n+1 ∗n the form ∆ (Φ2 − Φ1 ) = ikL∗n (R2 − R1 ). Hence it is important to try to get by with only jn = 1 iteration at each time step. We shall analyze the convergence of the iterative scheme in (4.12) and at the same time demonstrate the attendant stability of the resulting new fully discrete method. To this end assume for the time being that we have available approximations Ujnn ∈ Sh to U n for n = 0, 1, 2, 3, 4 such that (4.13) kU n − Ujnn k ≤ cn k 4 + hr , n = 0, 1, 2, 3, 4.
For n ≥ 3, given Ujnn approximating U n , we shall compute the required starting values U0n,m , m = 1, 2, for the inner (j) loop in (4.12) using extrapolation from previous values, i.e. as 3 X n,m µmλ Ujn−λ , m = 1, 2, 3 ≤ n ≤ M − 1, (4.14) U0 = n−λ λ=0
where the constants µmλ , m = 1, 2, λ = 0, 1, 2, 3, will be computed by letting p(t) be the cubic polynomial interpolating to the values y n−λ = y(tn−λ) of a smooth function y(t) at the points tn−λ , λ = 0, 1, 2, 3, and setting n
p(t + τm k) =
3 X
µmλ y n−λ. m = 1, 2.
λ=0
Proposition 4.1. Let u be sufficiently smooth, U n , 0 ≤ n ≤ M, be the solution of (4.2) and suppose that there exist Ujnn ∈ Sh , 0 ≤ n ≤ 4, such that (4.13) holds. For 4 ≤ n ≤ M −1 define Ujn+1 by the scheme (4.12) with starting values U0n,m as in (4.14). n+1 Then if jn+1 = 2 for all n, there exists k0 > 0 such that for k ≤ k0 (4.15) max kU n − Ujnn k ≤ c k 4 + hr . 0≤n≤M
If jn+1 = 1 for all n, then, given 0 < ε ≤ 1, there exists kε > 0 such that for k ≤ kε there holds (4.16) max kU n − Ujnn k ≤ c k 4−ε + hr . 0≤n≤M
In (4.15) and (4.16) c is a constant independent of k, h and ε.
Proof. Given Ujnn , let U˜ n,m , m = 1, 2, denote the exact solution of the system (4.17)
U˜ n,m = Ujnn + ik
2 X
˜ n,j , amj Ln,j h U
j=1
which we can write using (4.5) as (4.18)
Jn
U˜ n,1 U˜ n,2
!
=
! Ujnn . Ujnn
m = 1, 2,
20
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
First we prove a preliminary estimate that implies the convergence of the sequence ˜ n,m . From (4.12) and (4.18) suppressing n we let Yj,m = Ujn,m , j = 0, 1, 2, . . . to U Ujn,m − U˜ n,m , Yj = (Yj,1, Yj,2)T and obtain J∗n Yj+1 = (J∗n − Jn )Yj =: Zj = (Zj,1, Zj,2)T , i.e. Yj+1 = ikAL∗n Yj+1 + Zj ,
(4.19)
where L∗n = diag (L∗n , L∗n ) on (Sh )2 . There follows by Lemma 2.1 that |||Yj+1||| ≤ c|||Zj |||. On the other hand, using e.g. (3.6), we obtain kZj,mk ≤ ck 2 kYj,1k+kYj,2k), m = 1, 2 and conclude therefore that (4.20)
n,m ˜ n,m k ≤ γk 2 max kU n,m − U˜ n,m k, max kUj+1 −U j
m=1,2
m=1,2
j = 0, 1, . . . , jn+1 − 1,
where γ is a constant independent of h, k and the choice of jn . We next estimate the difference kU0n,m − U˜ n,m k. (In what follows we let n ≥ 4). Using (4.14) we have for m = 1, 2
(4.21)
3 X n,m n,m n,m n,m n,m n,m n,m ˜ ˜ − U ) + (U −u )+ u − µmλ un−λ U − U0 = (U λ=0
+
3 hX λ=0
µmλ u
n−λ
−U
n−λ
i
+
3 hX
µmλ U
n−λ
λ=0
−
Ujn−λ n−λ
i .
First, since U˜ n,m satisfies (4.17) and U n,m (4.2b) we obtain, again by stability (Lemma 2.1), that ˜ n,m − U n,m k ≤ ckU n − U n k. max kU jn
(4.22)
m=1,2
Next, recalling the definition of V n,m from (3.14), write, with W n,m = W (tn,m), (4.23)
U n,m − un,m = (U n,m − V n,m ) + (V n,m − W n,m) + (W n,m − un,m)
and observe that by (3.14b), (4.2b), Lemma 2.1 and (4.3) (4.24)
max kU n,m − V n,m k ≤ ckU n − W n k ≤ c(k 4 + hr ).
m=1,2
Now from (3.19), (3.20) V n,m − W n,m = en,m + (Λm W n − W n,m). By (3.38) and (3.36) we have ken,m k ≤ c(k 4 + khr ), m = 1, 2, whereas (3.18) and Taylor’s theorem yield for m = 1, 2 4
X
n n,m
τmj k j W (j)n − W n,m ≤ ck 3 . kΛm W − W k = j=0
We conclude that (4.25)
kV n,m − W n,m k ≤ c(k 3 + khr ),
m = 1, 2,
and, therefore, by (4.23)–(4.25) that for k ≤ 1, (4.26)
kU n,m − un,m k ≤ c(k 3 + hr ),
m = 1, 2.
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
21
Finally, the definition of the extrapolation procedure yields 3
n,m X
u µmλ un−λ ≤ ck 4 . −
(4.27)
λ=0
We conclude, by (4.21), (4.22), (4.26) and (4.27) that there exist positive constants η and θ, independent of k, h and the choice of jn such that (4.28)
˜ n,m − U0n,m k ≤ θ(k 3 + hr ) + ηkU n − U n k + max kU jn
m=1,2
3 X λ=0
k, µλ kU n−λ − Ujn−λ n−λ
where µλ = maxm=1,2 |µmλ |, 0 ≤ λ ≤ 3. We come now to the main part of the proof which is an induction √ 2 −1 step on n. First we treat the case jn = 2 for all n ≥ 4. We let k ≤ k0 = 2 3 γ where γ is as in (4.20). Assume that for 4 ≤ m ≤ n, (H1)
kU m − Ujmm k ≤ cm (k 4 + hr ),
where cm are positive constants satisfying (H2)
cm = [1 + k 3 (η + µ0 )]cm−1 + k 3 (µ1 cm−2 + µ2 cm−3 + µ3 cm−4 ) + θk 2 .
Clearly, it may be arranged, by taking c3 or c4 large enough in (4.13), that (H1)–(H2) hold for n = 4. Define U˜ n+1 ∈ Sh , conformal to the notation in (4.17) as √ ˜ n+1 = U n + 3(U˜ n,2 − U˜ n,1 ) (4.29) U jn and split (4.30)
U n+1 − Ujn+1 = (U n+1 − U˜ n+1 ) + (U˜ n+1 − Ujn+1 ). n+1 n+1
Since the time stepping procedure is conservative in the L2 sense, subtracting (4.17) from (4.2b) and (4.29) from (4.2c′ ) we obtain (4.31)
˜ n+1 k = kU n − U n k. kU n+1 − U jn
On the other hand, using (4.29) and (4.12) we have √ ˜ n+1 − U n+1 k = k 3(U˜ n,2 − U n,2 ) − (U˜ n,1 − U n,1 )k kU jn+1 jn+1 jn+1 √ n,m n,m ˜ − Ujn+1 k ≤ 2 3 max kU (4.32) m=1,2 √ ˜ n,m − U n,m k ≤ 2 3 γ 2 k 4 max kU 0 m=1,2
where in the last inequality we used the fact that jn+1 = 2 and (4.20). We conclude therefore by (4.30)–(4.32) that √ 2 4 n n ˜ n,m − U n,m k. 3 γ k max kU (4.33) kU n+1 − Ujn+1 k ≤ kU − U k + 2 0 jn n+1 m=1,2
√
Hence, (4.33) and (4.28) yield, if k ≤ k0 = 2 3 γ kU n+1 − Ujn+1 k ≤ (1 + ηk 3 )kU n − Ujnn k + k 3 n+1
3 X λ=0
2 −1
k + θk 2 (k 4 + khr ). µλ kU n−λ − Ujn−λ n−λ
22
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
Therefore, using the induction hypothesis (H1) we obtain kU n+1 − Ujn+1 k n+1 3 ≤ 1 + k (η + µ0 ) cn + k 3 (µ1 cn−1 + µ2 cn−2 + µ3 cn−3 ) + θk 2 (k 4 + khr ).
Hence, (H1)–(H2) hold for m = n+1 as well and the inductive step is complete. Clearly the constants cn are uniformly bounded for 0 ≤ n ≤ M by a constant independent of k and h; (4.15) follows. We examine now the consequences of taking jn = 1 for all n ≥ 4. Given 0 < ε ≤ 1 we √ −1/ε shall assume that k ≤ kε = 2 3 γ for reasons that will become apparent below. Our induction hypotheses that replace (H1), (H2) are now that for 4 ≤ m ≤ n, kU m − Ujmm k ≤ cm (k 4−ε + hr ),
(H1′ )
where cm are positive constants given by (H2′ )
cm = [1 + k 2−ε (η + µ0 )]cm−1 + k 2−ε (µ1 cm−2 + µ2 cm−3 + µ3 cm−4 ) + θk.
The verification of the inductive step follows the lines of the previous proof: (4.29) ˜ n+1 − U n+1 k ≤ to (4.31) still hold of course but since jn+1 = 1, (4.32) becomes kU jn+1 √ ˜ n,m − U0n,m k. Consequently we have kU n+1 − U n+1 k ≤ kU n − 2 3 γk 2 maxm=1,2 kU jn+1 √ −1/ε √ n,m 2 n,m n ˜ − U0 k, and therefore for k ≤ kε = 2 3 γ U k + 2 3 γk maxm=1,2 kU we jn
obtain by (H1′ )
h √ 2 n n 3 γk θ(k 3 + hr ) + ηkU n − Ujnn k k + 2 kU n+1 − Ujn+1 k ≤ kU − U jn n+1 +
3 X λ=0
≤
µλ kU n−λ − Ujn−λ k n−λ
i
1 + (η + µ0 )k 2−ε cn + k 2−ε (µ1 cn−1 + µ2 cn−2 + µ3 cn−3 ) (k 4−ε + hr )
+ θk(k 4−ε + hr ),
verifying that (H1′ ), (H2′ ) hold for m = n + 1 as well. Obviously the constants cn can be made uniformly bounded in n by a constant independent of ε, albeit larger than the corresponding constant in the jn = 2 case; we conclude that (4.16) holds. It is easy to construct Ujnn , 0 ≤ n ≤ 4, such that (4.13) is valid. For n = 0 we already have stipulated that Uj00 = U 0 . For n = 1 only the previous value U 0 is available. We set U00,1 = U00,2 = U 0 in (4.12) and generate the sequence Uj0,m taking j1 = 2. This will suffice since the analog of (4.27) is an O(k) bound implying that in the analog of (4.28) maxm=1,2 kU˜ 0,m − U00,m k = O(k + hr ). Similarly for n = 2, 3, 4 it suffices to generate Ujnn from (4.12) with jn = 1, computing for each n, U0n−1,m as the appropriate linear combination of the previous values Uj00 , Uj11 , . . . , Ujn−1 using the Lagrange interpolating n−1 polynomial of degree n − 1. In practice we noticed that taking cubic polynomial extrapolation to generate the starting values for n ≥ 4 and just jn = 1 was generally sufficient to preserve the overall order of accuracy and stability of the scheme. We report in [1] these and other relevant numerical experiments that we performed with the method, including experiments in
¨ CONSERVATIVE GALERKIN METHODS FOR THE SCHRODINGER EQUATION
23
which the operator J∗n is not evaluated at every time step but rather every m∗ > 1 time steps. Our experiments indicate that it is possible to take in many interesting examples m∗ equal to, say, 20 and jn = 2 (the m∗ “large”, jn = 1 combination was unstable for some hard to integrate problems) and still preserve the overall order of accuracy and stability of the scheme. Note added in proof : For more recent work related to error estimates for the Nonlinear Schr¨odinger Equation we refer the reader to [24], [25]. References 1. G.D. Akrivis and V.A. Dougalis, On a conservative, high-order accurate finite element scheme for the ‘parabolic’ equation, to appear in the Proceedings for the 2nd IMACS symposium on computational acoustics, Princeton University, 15-17 March 1989. 2. G.A. Baker, J.H. Bramble and V. Thom´ee, Single step Galerkin approximations for parabolic problems, Math. Comp. 31 (1977) 818–847. 3. J.L. Bona, V.A. Dougalis, O.A. Karakashian and W. McKinney, Conservative high order schemes for the generalized Korteweg-de Vries equation, to appear. 4. A. Broc´ehn, Galerkin methods for approximation of solutions of second order partial differential equations of Schr¨ odinger type, Ph.D. Thesis, University of G¨oteborg, 1980. 5. J.C. Butcher, Implicit Runge-Kutta processes, Math. Comp. 18 (1964) 50–64. 6. J.C. Butcher, The numerical analysis of ordinary differential equations; Runge-Kutta methods and general linear methods, John Wiley, Chichester, 1987. 7. M. Crouzeix, Sur la B-stabilit´e des m´ethodes de Runge-Kutta, Numer. Math. 32 (1979) 75–82. 8. M. Crouzeix and V. Thom´ee, On the discretization in time of semilinear parabolic equations with nonsmooth initial data, Math. Comp. 49 (1987) 359–377. 9. K. Dekker and J.G. Verwer, Stability of Runge-Kutta methods for stiff nonlinear differential equations, North Holland, Amsterdam, 1984. 10. E. Dendy, Jr., An alternating direction method for Schr¨ odinger’s equation, SIAM J. Numer. Anal. 14 (1977) 1028–1032. 11. V.A. Dougalis and O.A. Karakashian, On some high order accurate fully discrete Galerkin methods for the Korteweg-de Vries equation, Math. Comp. 45 (1985) 329–345. 12. O.A. Karakashian and W. McKinney, On optimal high order in time approximations for the Korteweg-de Vries equation, Math. Comp. (to appear). 13. D. Lee and S.T. McDaniel, Ocean acoustic propagation by finite difference methods, Comput. Math. Appl. 14 No. 5. 14. D. Lee, R.L. Sternberg and M.H. Schultz, eds., Computational acoustics: wave propagation, Proceedings of the 1st IMACS symposium on computational acoustics, New Haven, 6–8 August 1986, vols.1,2, North Holland, Amsterdam, 1988. 15. J.L. Lions and E. Magenes, Probl`emes aux limites non homog`enes et applications, vol. 2, Dunod, Paris, 1968. 16. A. Quarteroni, Mixed approximations of evolution problems, Comput. Meths. Appl. Mech. Engrg. 24 (1980) 137–163. 17. J.M. Sanz-Serna, Methods for the numerical solution of the nonlinear Schr¨ odinger equation, Math. Comp. 43 (1984) 21–27. 18. J.M. Sanz-Serna and J.G. Verwer, Conservative and nonconservative schemes for the solution of the nonlinear Schr¨ odinger equation, IMA J. Numer. Anal. 6 (1986) 25–42. 19. M.H. Schultz and D. Lee, eds., Computational ocean acoustics, Invited lectures from the workshop held at Yale University, 1-3 August 1984, Comput. Math. Appl. 11 (1985) Nos 7-8.
24
GEORGIOS D. AKRIVIS AND VASSILIOS A. DOUGALIS
20. F.D. Tappert, The parabolic approximation method. In: Wave propagation and underwater acoustics, J.B. Keller and J.S. Papadakis, eds., pp. 224–287, Lecture Notes in Physics v. 70, SpringerVerlag, Berlin-Heidelberg, 1977. 21. V. Thom´ee, Convergence estimates for semi-discrete Galerkin methods for initial-value problems. In: Numerische, insbesondere approximations-theoretische Behandlung von Funktionalgleichungen, R.Ansorge and W. Tornig, eds., pp. 243–262, Lecture Notes in Mathematics v. 333, SpringerVerlag, Berlin-Heidelberg, 1973. 22. V. Thom´ee, Galerkin finite element methods for parabolic problems, Lecture Notes in Mathematics v. 1054, Springer-Verlag, Berlin-Heidelberg, 1984. 23. L.B. Wahlbin, A dissipative Galerkin method for the numerical solution of first order hyperbolic equations. In: Mathematical aspects of finite elements in partial differential equations, C. de Boor, ed., pp. 147–169, Academic Press, New York, 1974. 24. G.D. Akrivis, V.A. Dougalis and O.A. Karakashian, On fully discrete Galerkin methods of secondorder temporal accuracy for the Nonlinear Schr¨ odinger Equation, to appear in Numer. Math. 25. O. Karakashian, G.D. Akrivis and V.A. Dougalis, On optimal-order error estimates for the Nonlinear Schr¨ odinger Equation, to appear. Mathematics Department, University of Crete, Iraklion, Greece Mathematics Department, University of Crete, Iraklion, Greece
