Stiff convergence of force-gradient operator splitting methods
Stiff convergence of force-gradient operator splitting methods Emil Kieri∗ February 28, 2014
Abstract We consider force-gradient, also called modified potential, operator splitting methods for problems with unbounded operators. We prove that force-gradient operator splitting schemes retain their classical orders of accuracy for time-dependent partial differential equations of parabolic or Schr¨ odinger type, provided that the solution is sufficiently regular.
1
Introduction
We consider the numerical time-integration of differential equations on the form ut = Au + Bu,
u(0) = u0 ,
(1)
evolving on a Banach space X, with the operators A : D(A) ⊂ X → X, B : D(B) ⊂ X → X. Splitting methods are advantageous if the subproblems with only one of the operators can be solved easily. An example of such a problem, where splitting methods often are used, is the time-dependent Schr¨odinger equation, iut = −∆u + V (x)u, u(x, 0) = u0 ∈ L2 (Rd ). If one of the terms is neglected we get the subproblems ivt = −∆v, iwt = V (x)w,
⇒
v(x, 0) = v0 ,
⇒
w(x, 0) = w0 ,
v(x, t) = eit∆ v0 , w(x, t) = e
−itV (x)
and w0 .
The operator exponentials are readily computed since the potential V (x) and the Laplacian ∆ are diagonal in coordinate and frequency space, respectively. The combined operator (−∆ + V (x)), however, is not as easily exponentiated. Splitting methods consist in taking small time steps, alternatingly applying the action of the two operators. An s-stage standard splitting method for (1) can be written as un+1 = Sun ,
S=
s Y
ebj hB eaj hA ,
(2)
j=1 ∗ Division of Scientific Computing, Department of Information Technology, Uppsala University, Box 337, 751 05 Uppsala, Sweden ([email protected]).
1
where h is the time step and the coefficients aj , bj define the splitting scheme. Throughout this paper products are defined in falling order, i.e., s Y
Aj = As As−1 . . . A1 ,
j=1
Ql and j=k aj = 1 if k > l. We let C denote a generic constant which may take different values at different occurrences. We say that a splitting method is of classical order p if the local error estimate k(S − eh(A+B) )uk ≤ Chp+1 holds for bounded operators A and B, and sufficiently small time steps h. Such estimates are obtained using Taylor expansion or the Baker–Campbell– Hausdorff formula [11, Ch. III.5]. They are adequate, for instance, for ordinary differential equations, with u ∈ RN , and A, B ∈ RN ×N . In infinite-dimensional settings, classical error analysis is no longer rigorous since the possibly unbounded operators A and B appear in the remainder. Still, splitting methods often retain their classical order also for such problems. Error estimates which hold independently of the norms of the operators are called stiff estimates. In an early work, Jahnke and Lubich proved that the Strang splitting scheme keeps its second order accuracy for Schr¨odinger equations [15]. Thalhammer [22], and Hansen and Ostermann [13], proved that splitting methods keep their classical order for certain classes of problems with unbounded operators, given that the solution is sufficiently regular. See also [5, 14], which enables splitting methods of order higher than two for parabolic problems using the theory of analytic semigroups. Using the calculus of Lie derivatives, the theory was extended to non-linear problems in [16, 23]. If the problem has a particular structure, this can often be used to construct more efficient methods. One example is Runge–Kutta–Nystr¨om methods [3, 11], which apply when commutators of the operators satisfy [B, [B, [A, B]]] = 0. This reduces the number of order conditions, and gives more freedom to the choice of method coefficients. In this paper we will study force-gradient methods, which were proposed by Suzuki [20, 21]. They rely on the same commutator assumption as Runge–Kutta–Nystr¨om methods, and in addition that [B, [A, B]] can be conveniently evaluated. A general force-gradient method for the problem (1) reads s Y 3 un+1 = Sun , S = ebj hB+cj h [B,[A,B]] eaj hA . (3) j=1
The perhaps most used force-gradient scheme is the fourth order method of Chin [6, 8, 9], with s = 3. An attractive property of the Chin scheme is that all its coefficients are real and non-negative. The method can therefore be used both for parabolic and Schr¨ odinger-type problems. Force-gradient methods of order higher than four have been derived, e.g., by Omelyan et al. [18]. These methods do however have some negative coefficients. This makes them unsuitable for parabolic problems, which cannot be solved backwards in time. Chin has proven that a 6th order force-gradient method necessarily has negative coefficients, unless higher order commutators are included in the scheme [7]. As for standard splitting methods, this order barrier can be circumvented if the coefficients are 2
allowed to be complex-valued. In [1], Bader et al. derive several high order forcegradient schemes with complex coefficients with positive real parts, suitable for parabolic problems. In this paper we derive stiff order conditions for force-gradient methods. That is, conditions under which the methods are of a certain order of accuracy also for unbounded operators. We then prove that the stiff order conditions are equivalent to the classical ones, as long as the exact solution is sufficiently regular. The result holds for parabolic and Schr¨odinger-type partial differential equations. Problems where force-gradient methods are applicable can be found in those classes, since [V, [∆, V ]] = −2|∇V |2 . The remainder of this paper is organised as follows. In Section 2 we specify the assumptions on the operators and the regularity requirements, and state local and global error estimates. The global estimate follows directly from the local estimate and the assumptions. The local error estimate, which is the main result of this paper, is proven in Section 3. In Section 4 we use the calculus of Lie derivatives to show how force-gradient methods can be used for nonautonomous problems. We conclude the paper with numerical experiments in Section 5, which corroborate the theoretical results.
2
Statement of the error estimate
In this section we list the assumptions we rely on and prove the global error estimate, Theorem 2. The global estimate relies on a local error estimate, Theorem 1, which we will prove in Section 3. The first assumption specifies the classes of problems under consideration. In particular, linear autonomous problems of parabolic and Schr¨ odinger type are included. Assumption 1. The problem (1) is assumed to be well-posed, and the commutators of A and B are assumed to satisfy [B, [B, [A, B]]] = 0. We further assume either of the following two conditions. 1. (Schr¨ odinger-type problem) The operators (A + B), A, B and [B, [A, B]] are the infinitesimal generators of strongly continuous groups of operators {etF }t∈R . Furthermore, there exists a constant ω ≥ 0 such that ketF kX←X ≤ eω|t|
∀ t ∈ R,
F = (A + B), A, B, and [B, [A, B]]. 2. (Parabolic problem) The operators (A + B), A, B, and [B, [A, B]] are the infinitesimal generators of analytic semigroups {etF }t∈Σθ , defined in the sectors Σθ = {t ∈ C : −θ ≤ arg t ≤ θ} for some 0 < θ ≤ π/2. Furthermore, there exists a constant ω ≥ 0 such that ketF kX←X ≤ eω|t|
∀ t ∈ Σθ ,
F = (A + B), A, B, and [B, [A, B]]. Next, we assume that the solution and the operators A and B are sufficiently regular. This assumption was justified for Schr¨odinger equations with sufficiently regular, bounded potential B in [22, Sec. 2]. The arguments transfer directly to parabolic problems on the form ut = ∆u + V (x)u. We here use the 3
notation |µ| = |µ1 | + · · · + |µk | for the modulus of a multi-index, and iterated commutators are denoted p adp+1 A (B) = [A, adA (B)],
ad0A (B) = B.
Assumption 2. For a force-gradient splitting method of classical order p, assume that k
Y
µ adAj (Bj )eζj A
≤ C,
X←Yp
j=1
µ ∈ Nk , |µ| = p + 1 − k,
Bj = B and [B, [A, B]],
j = 1, . . . , k,
(4) (5)
holds for k = 1, . . . , p + 1, with ζj ∈ R (Schr¨ odinger) or ζj ∈ Σθ (parabolic), |ζj | ≤ C, and some suitably chosen normed space Yp ⊂ X. Assume further that keζF kYp ←Yp ≤ C,
F = (A + B), A, B, and [B, [A, B]],
(6)
and that the exact solution to (1) resides in the space Yp . For a problem with smooth, bounded potential we get the standard Sobolev spaces Yp = H p . Under these assumptions we can prove the following local and global error estimates. Theorem 1. Assume that Assumptions 1 and 2 hold, and that the coefficients aj , bj , cj define a force-gradient operator splitting scheme of classical order p. Assume further that the coefficients reside in R if (1) is of Schr¨ odinger type, or in Σθ if parabolic. Then, the scheme satisfies the local error estimate kSu(t) − u(t + h)kX ≤ Chp+1 , also for unbounded operators A and B. Proof. We will prove this result in Section 3. Theorem 2. Assume the conditions of Theorem 1. Then, the force-gradient operator splitting scheme (3) satisfies the global error estimate kS n u0 − u(nh)kX ≤ Chp ,
0 ≤ nh ≤ T.
Proof. By Assumption 1, kSkX←X ≤ ech for some c. Using this in connection with the local error estimate, the result follows through a standard Lady Windermere’s fan argument [12, Ch. II.3]. See, e.g., [13, Thm. 2.3] for a proof of a much similar result.
3
Proof of the local error estimate
This section is devoted to the proof of Theorem 1. The line of attack is as follows. First, we summarise the proof of a similar result [22], a local error estimate for standard splitting methods on the form (2). This introduces notation and techniques which we will use when we extend the proof to force-gradient methods. Second, we take the commutators [B, [A, B]] into account, and manipulate 4
the expressions so that they fit in the framework of the stiff order conditions for standard methods. Finally, we observe that the stiff order conditions hold if and only if the classical order conditions do. Since regularity is assumed uniformly in time, we can without restriction consider the error in the first time step, e0 = (S − eh(A+B) )u0 . As a first step, we express the exact and approximate solutions in forms that facilitate comparison. Repeating calculations from [22], we express the exact solution after one time step as u(h) = ehA u0 +
p X
Z Ik u0 + Rp+1 ,
Ik =
gk (τ ) dτ ,
(7)
∆k
k=1
with τ = (τ1 , . . . , τk ), τ0 = h, and k Y
fk (τ ) =
e(τk−l −τk−l+1 )A B,
gk (τ ) = fk (τ )eτk A .
l=1
The integrals are evaluated over the domains ∆k = {τ ∈ Rk : 0 ≤ τk ≤ τk−1 ≤ · · · ≤ τ1 ≤ h}, and the remainder term reads Z Rp+1 =
fp+1 (τ )u(τp+1 ) dτ . ∆p+1
Throughout this section, we let R0p+1 and Rp+1 denote arbitrary linear combinations of terms on the form R0p+1 ∼ hp+1
k Y
µ adAj (Bj )eζj A χ,
Rp+1 = R0p+1 u,
|µ| ≤ p − k + 1,
j=1
with the linear operator χ bounded Yp → Yp . By Assumption 2, this means that kRp+1 kX ≤ Chp+1 . Note that Rp+1 = Rp+1 . Now, if we define Aj = aj hA,
Bj = bj hB,
Cj = cj h3 [B, [A, B]],
the approximate solution can be expressed similarly, S=
s Y j=1
e
Aj
+
p X
Qk + rp+1 .
k=1
Here, Qk =
X
κλ G(λ),
Λk = {λ ∈ Nk : 1 ≤ λk ≤ λk−1 ≤ · · · ≤ λ1 ≤ s},
λ∈Λk
F (λ) =
k Y
λk−l
Y
eAj (Bλk−l+1 + Cλk−l+1 ),
G(λ) = F (λ)
λk Y j=1
l=1 j=λk−l+1 +1
5
eAj ,
with λ = (λ1 , . . . , λk ) ∈ Nk and λ0 = s. The remainder term reads rp+1 =
p+1 X X
κ0j,λ F (λ)ϕj (Bλk + Cλk )eAλk
λY k −1
eBj +Cj eAj .
j=1
j=1 λ∈Λp+1
The function ϕj (z), which appear in the remainder, is defined recursively as ϕj (z) =
1 + zϕj+1 (z), j!
j ≥ 0,
ϕ0 (z) = ez ,
(8)
and κλ and κ0j,λ are well-behaved expansion coefficients. The expression ϕj (Bλk + Cλk )e
A λk
λY k −1
eBj +Cj eAj
j=1
is a bounded operator Yp → Yp due to (6), and therefore rp+1 = R0p+1 . The values of κ0λ are of lesser importance, but κλ appears in the order conditions. If σ0 = 0,
λσ0 +···+σl−1 +1 = · · · = λσ0 +···+σl ,
l = 1, . . . , m,
and σ0 + · · · + σm = k, then κλ =
1 . σ0 ! · · · σm !
The local error can now be expressed as p s Y X e0 = ehA − eAj u0 + (Ik − Qk )u0 + Rp+1 . j=1
(9)
k=1
In [22] it was shown how each of the terms in (9) could be bounded separately, yielding the desired result for a standard operator splitting scheme. For forcegradient methods this is not possible since some of the order conditions for standard splitting schemes are violated. Error cancellation between the terms will instead recover the desired order of accuracy. We start off ignoring this, and re-derive the order conditions for a standard method. After that, we investigate how corrections from the [B, [A, B]]-terms kick in. In what is coming, we need some notation for multi-indices. For µ ∈ Nk and x ∈ Ck , define µ! = µ1 !µ2 ! · · · µk ! and xµ = xµ1 1 xµ2 2 · · · xµk k . Furthermore, if 1 ≤ µj ≤ l and x ∈ Cl , let xµ = (xµ1 , . . . , xµk ). The first error term, s Y (0) e0 = ehA − eaj hA u0 , j=1
Ps is zero as long as the consistency condition j=1 aj = 1 is fulfilled. Taylor expansion of the integrands in (7) gives Z X 1 Z Ik = τ µ dτ ∂τµ gk (0) + %k,p−k+1 (τ ) dτ , (10) µ! ∆k ∆k |µ|≤p−k
6
with the remainder X
%k,n+1 =
ν∈Nk |ν|=n+1
n+1 ν!
Z
1
(1 − s)n ∂τν gk (sτ )τ ν ds,
0
and %k,n+1 = R0n+1 . The expansion coefficients in (10) can be evaluated as Z
τ µ dτ = hk+|µ|
k Y
∆k
l=1
1 , µl + · · · + µk + k − l + 1
and the partial derivatives of gk are ∂τµ gk (τ ) = (−1)|µ|
k Y
µ e(τk−l −τk−l+1 )A adAk−l+1 (B) eτk A .
(11)
l=1
Next, we let (0)
Qk = Qk |cj ≡0 , Pj and Taylor expand that. Defining αj = i=1 ai , we get (0) Qk
k
=h
X
κλ
k Y
λ∈Λk
=
X
l=1
h
k+|µ|
|µ|≤p−k
+ hk
bλl gk (αλ h)
X
k Y 1 X µ µ κλ bλl αλ ∂τ gk (0)+ µ! λ∈Λk
κλ
λ∈Λk
k Y
l=1
bλl %k,p−k+1 (αλ h).
l=1
Defining k
Πµ =
1 Y 1 , µ! µl + · · · + µk + k − l + 1 l=1
Σ(0) µ =
k Y 1 X µ κλ bλl αλ , µ! λ∈Λk
l=1
the kth order error term reads X (0) Ik − Qk = hk+|µ| Πµ − Σ(0) ∂τµ gk (0) + R0p+1 . µ |µ|≤p−k
We thereby achieve pth order of accuracy if αs = 1 and X µ Πµ − Σ(0) k = 1, . . . , p, σ = 0, . . . , p − k. µ ∂τ gk (0) = 0,
(12)
k
µ∈N |µ|=σ
These are the stiff order conditions for standard operator splitting methods derived in [22]. For force-gradient methods these conditions are typically violated, since the contributions from Cj are not yet taken into account. We will next extend the proof to address that. The key observation is that changing a B 7
(0)
for [B, [A, B]] in the expansion of Qk transforms gk into two first order partial derivatives of gk+1 . Take, for example, 3 Y
g2 (α(2,1) h) =
eaj hA B
j=3
2 Y
1 Y
eaj hA B
j=2
eaj hA ,
j=1
and compare it to 3 Y
=
eaj hA B
2 Y
j=3
j=2
3 Y
2 Y
eaj hA B
j=3
−
eaj hA B
1 Y
eaj hA
j=1 1 Y
eaj hA B[A, B]
j=2
3 Y j=3
=
eaj hA [B, [A, B]]
eaj hA +
j=1
2 Y
eaj hA [A, B]B
j=2
1 Y
eaj hA
j=1
−∂τ(0,0,1) g3 (α(2,1,1) h)
+
∂τ(0,1,0) g3 (α(2,1,1) h),
cf. (11). We introduce the operators Dl± as µ
Dl± ∂τµ gk (τ ) = ∂τ ±,l gk (τ 0 ), with µ ∈ Nk , τ ∈ Rk , and µ+,l = (µ1 , . . . , µl−1 , µl , µl + 1, µl+1 , . . . , µk ) ∈ Nk+1 , µ−,l = (µ1 , . . . , µl−1 , µl + 1, µl , µl+1 , . . . , µk ) ∈ Nk+1 , τ 0 = (τ1 , . . . , τl−1 , τl , τl , τl+1 , . . . , τk ) ∈ Rk+1 . With this notation we can write k
Qk = h
X λ∈Λk
cλ
k Y
+ − bλk−l+1 − h2 cλk−l+1 Dk−l+1 + h2 cλk−l+1 Dk−l+1 gk (αλ h).
l=1
(13) We expand the products in (13), and order the terms with respect to the number of D-operators in them. We let Qk =
k X
(m)
Qk ,
m=0 (0)
with Qk as before and (m)
Qk
X
= hk+2m
λ∈Λk
×
m Y
κλ
k X l1 =m
···
lm−1 −1 k X Y lm =1
l=1 l6=lj
− + − Dk−l cλk−lj +1 Dk−l j +1 j +1
j=1
8
bλk−l+1 ×
gk (αλ h),
m = 1, . . . k.
The point with this notation is that the error now can be expressed as
e0 = e
hA
−
s Y
k
e
Aj
b2c p X X (m) Ik − u0 + Qk−m u0 + Rp+1 .
j=1
m=0
k=1
The disregarded terms p X
k X
(m)
Qk u0
k=0 m=p−k+1
are of the form Rq , with q ≥ p + 1, and can thereby be controlled. With this formulation, we require no cross-cancellation between the terms k
b2c X
Ik −
(m) Qk−m .
m=0 (m)
In order to match terms with the same derivative of gk we Taylor expand Qk−m . (1)
For illustration, we first expand Qk−1 . (1)
X
Qk−1 = h(k−1)+2
k−1 X
k−1 Y
l1 =1
l=1 l6=l1
κλ
λ∈Λk−1
bλ(k−1)−l+1 cλ(k−1)−l1 +1 ×
− + gk−1 (αλ h) × D(k−1)−l − D(k−1)−l 1 +1 1 +1 X
= hk+1
κλ
λ∈Λk−1 e
k−1 X
k−1 Y
l1 =1
l=1 l6=l1
ek−l1 +1
× ∂τ k−l1 − ∂τ X
=
hk+1+|µ|
|µ|≤p−(k+1) µ+ek−l1
µ × αλ (l) ∂τ
bλk−l cλk−l1 ×
gk (αλ(l) h)
1 µ!
X
κλ
λ∈Λk−1
µ+ek−l1 +1
− ∂τ
k−1 X
k−1 Y
l1 =1
l=1 l6=l1
bλk−l cλk−l1 ×
gk (0) + R0p+1 ,
with λ(l) = (λ1 , . . . , λl−1 , λl , λl , λl+1 , . . . , λk−1 ) ∈ Nk . We rewrite this as (1)
Qk−1 =
X
µ 0 hk+|µ| Σ(1) µ ∂τ gk (0) + Rp+1 ,
with
1≤|µ|≤p−k
Σ(1) µ =
k−1 X
1 (µ − eν )! ν=1 k X
X
κλ
λ∈Λk−1
1 − (µ − eν )! ν=2
k−1 Y
µ−eν + bλk−l cλν αλ (l)
l=1 l6=k−ν
X
κλ
λ∈Λk−1
k−1 Y
µ−eν bλk−l cλν−1 αλ . (l)
l=1 l6=k−ν+1
Terms with µ − eν ∈ / Nk are excluded from the summation. Similarly in the 9
general case, (m)
X
Qk−m = h(k−m)+2m
κλ
λ∈Λk−m
×
m Y
k−m X
lm−1 −1 k−m
···
l1 =m
X
Y
lm =1
l=1 l6=lj
bλ(k−m)−l+1 ×
− + cλ(k−m)−lj +1 D(k−m)−l − D(k−m)−l j +1 j +1
gk−m (αλ h)
j=1
X
= hk+m
λ∈Λk−m
×
m Y
k−m X
κλ
···
lm−1 −1 k−m
l1 =m
X
Y
lm =1
l=1 l6=lj
e(k−m)−lj +j
cλ(k−m)−lj +1 ∂τ
bλ(k−m)−l+1 ×
e(k−m)−lj +j+1
− ∂τ
gk (αλ(l) h)
j=1
X
=
k+m+|µ|
h
|µ|≤p−(k+m)
×
m Y
1 µ!
X
κλ
λ∈Λk−m
e(k−m)−lj +j
cλ(k−m)−lj +1 ∂τ
k−m X
···
lm−1 −1 k−m
l1 =m
X
Y
lm =1
l=1 l6=lj
e(k−m)−lj +j+1
− ∂τ
bλ(k−m)−l+1 ×
µ µ αλ (l) ∂τ gk (0)+
j=1
+ R0p+1 . λ(l) is defined analogously to λ(l) , with all the elements λlj , j = 1, . . . , m, doubled. In order to turn this expression inside out we need to introduce some notation. Let B = {0, 1}, and define the function ν : Nm × Bm → Nm as νj (l, v) = k − m − lj + j + vj ,
j = 1, . . . , m,
l = (l1 , . . . , lm ).
Then the product of binomials above can be expanded as m Y
e(k−m)−lj +j
∂τ
e(k−m)−lj +j+1
− ∂τ
=
X
(−1)|v| ∂τν(l,v) .
v∈Bm
j=1 (m) Qk−m
We can now rewrite the expression for as X (m) µ 0 Qk−m = hk+|µ| Σ(m) µ ∂τ gk (0) + Rp+1 ,
with
µ∈Nk m≤|µ|≤p−k
Σ(m) µ
=
k−m X
lm−1 −1
X
···
lm =1 v∈Bm
l1 =m
×
X
k−m Y
bλ(k−m)−l+1
(−1)|v| (µ − ν(l, v))!
m Y
X
κλ ×
λ∈Λk−m
µ−ν(l,v) cλ(k−m)−lj +1 αλ(l) .
j=1
l=1 l6=lj
As before, terms with µ − ν(l, v) ∈ / Nk are omitted from the summation. Assuming the consistency condition αs = 1, we can write the local error as e0 =
p X k=1
k
X k
µ∈N |µ|≤p−k
k+|µ|
h
Πµ −
b2c X m=0
10
Σ(m) ∂τµ gk (0)u0 + Rp+1 . µ
(14)
Noting that ∂τµ gk (0) = (−1)|µ| ehA
k Y
µ
adAk−l+1 (B),
l=1
we can write the stiff order conditions for force-gradient methods as1 X k
µ∈N |µ|=σ
Eµ
k Y
k
µ adAk−l+1 (B)
= 0,
Eµ = Πµ −
b2c X
Σ(m) µ ,
(15)
m=0
l=1
for all k = 1, . . . , p, σ = 1, . . . , p − k. Note the similarity with the order conditions (12) for standard splitting methods. Note also that terms in (15) can be matched using the properties of commutators, in particular the Jacobi identity and the assumed relation [B, [B, [A, B]]] = 0 of Runge–Kutta–Nystr¨om methods. Finally, by repeating the arguments from [22], one can show that the conditions (15) are equivalent to the classical order conditions. We assume that A and B are bounded operators, e.g., square matrices of fixed size, and claim that (15) are sufficient and necessary conditions for pth order convergence. They are sufficient conditions for bounded operators, since they are for unbounded. By Taylor expanding ehA in the partial derivatives of gk (0) in (14) and collecting terms of similar order in h, we see that the conditions are necessary also for non-stiff problems.
4
Generalisation to non-autonomous problems
Using the calculus of Lie derivatives [11, Ch. III.5], the results presented in this paper can be generalised to non-linear problems. That is, the proofs of the theorems are readily repeated, but computing the flows of the subproblems might not always be practical. The requirement that the commutator [B, [A, B]] should be easy to handle is also restrictive. The practical use of force-gradient methods for non-linear problems is therefore limited. The non-linear framework does, however, provide a convenient way of generalising the results to nonautonomous linear problems. That is the goal of this section. Given two non-linear operators A, B : D ⊂ X → X, the Lie derivative DA of A and its exponential are defined by DA Bu = B 0 (u)A(u),
t1 ,t2 e(t2 −t1 )DA Bu = B(ΨA u).
ΨtA1 ,t2 is here the exact flow of the differential equation ut = A(u) from t = t1 to t = t2 , and B 0 (u) denotes the Fr´echet derivative of B(u). Taking B as the identity mapping yields DA u = A(u),
t1 ,t2 e(t2 −t1 )DA u = ΨA u.
The Lie commutator [DA , DB ] will be an essential tool in the following. It is defined through [DA , DB ]u = DA DB u − DB DA u = B 0 (u)A(u) − A0 (u)B(u). 1A
Matlab script for evaluating the expressions Eµ is available upon request.
11
To generalise Theorems 1 and 2 to non-linear problems on the form ut = A(u) + B(u),
u(0) = u0 ,
we only need to assert that the problem satisfies a modified Assumption 2. In the modified assumption, (4)–(5) are replaced by k
Y
µ adDjA (DBj )eζj DA
≤ C,
X←Yp
j=1
µ ∈ Nk , |µ| = p + 1 − k,
DBj = DB or [DB , [DA , DB ]],
j = 1, . . . , k,
for k = 1, . . . , p + 1. Provided that this holds, the proofs can be carried out as previously, only with the order of the exponentials reversed. We will now show that this assumption holds for non-autonomous linear problems ut = A(t)u + B(t)u, u(0) = u0 , satisfying [A(s), A(t)] = [B(s), B(t)] = 0 for all s, t. We will furthermore show that the Lie commutator [DB , [DA , DB ]] has a simple form and can be used as before. We introduce time as an auxiliary variable and write the problem as u Ut = F (U ) + G(U ), with U = , s and
F (U ) =
A(s)u , 1
G(U ) =
B(s)u . 0
(16)
The Fr´echet derivatives of F and G are then A(s) A0 (s)u B(s) B 0 (s)u F 0 (U ) = , G0 (U ) = , 0 0 0 0 and the iterated Lie commutators can be determined as ! p f ad (B(s))u p A(s) adDF (DG )U = , 0 with p+1 p d fp f f ad (adA(s) (B(s))), A(s) (B(s)) = [adA(s) (B(s)), A(s)] + ds 0 f ad (B(s)) = B(s). A(s)
We can conclude that in order to fulfil the modified Assumption 2 we, in addition to the conditions from the autonomous case, must require the operators A and B to be temporally smooth. This is a reasonable requirement. The force-gradient correction similarly evaluates to dB(s) −[B(s), [B(s), A(s)]]u + dB(s) ds B(s)u − B(s) ds u [DG , [DF , DG ]]U = 0 [B(s), [A(s), B(s)]]u = . (17) 0 12
Since B(s) commutes with itself evaluated at any other time, it also commutes with its time derivative. The commutator in (17) is the same commutator as in the autonomous problem, evaluated at a single point in time. Furthermore, since the operator A evaluated at different times commutes, its exponentiation is given by the first term of its Magnus series, Z t2 t1 ,t2 ΨA = exp A(σ) dσ . t1
In a computational context, it is sufficient to approximate the integral in the exponent with a pth order quadrature rule. Using the recurrence relation (8) of the functions ϕj the quadrature remainder can be isolated from the scheme. As an example, consider Strang splitting with the integral of A(t) approximated with the midpoint rule. Writing the quadrature rule as the exact integral and a truncation error, we get Z tn +h ehB/2 exp A(σ) dσ + Ch3 A00 (ξ) ehB/2 = tn
e
hB/2
exp
Z
tn +h
A(σ) dσ ehB/2 +
tn
+ Ch3 ehB/2 exp
Z
tn +h
A(σ) dσ A00 (ξ)ϕ1 Ch3 A00 (ξ) ehB/2 ,
tn
with ξ ∈ [tn , tn + h]. Note that it is possible to introduce time as an auxiliary variable in several ways [2]. One may even prefer to use more than one auxiliary variable. In connection with force-gradient methods, however, the arrangement (16) seems the most practical. The reason for this is that we want time to “stand still” while G(U ) is integrated, otherwise the commutator (17) will become complicated. When using schemes with complex coefficients for non-autonomous problems, one may have to evaluate the operators at complex times. This can lead to poor conditioning. In order to circumvent this difficulty, Blanes and Seydao˘ glu suggest methods with aj ∈ R+ , bj ∈ Σθ [4]. No force-gradient method with this property (and cj ∈ Σθ ) of order higher than four has to our knowledge been proposed, but could in principle be derived.
5
Numerical experiments
Preservation of the classical order for force-gradient schemes applied to partial differential equations have been observed previously. The Chin scheme, with s = 3 and 1 2 1 1 , a1 = 0, a2 = a3 = , b1 = b3 = , b2 = , c1 = c3 = 0, c2 = 2 6 3 72 (18) was studied for the Schr¨ odinger equation in [8, 9]. In [19], fourth order accuracy was observed for a non-autonomous Schr¨odinger-like problem in four dimensions. For completeness, we present convergence studies for the Chin scheme applied to the two-dimensional Schr¨ odinger-type and parabolic problems 1 iut = − ∇ · A(t)∇u + V (x, t)u, 2 13
ut = ∇ · A(t)∇u + V (x, t)u,
on x ∈ [−π, π]2 . We use periodic boundary conditions and the Gaussian initial condition T u(x, 0) = e−5 x x . The matrix A and the potential V are given by t e 0 A(t) = , V (x, t) = 2 − cos(x1 ) − sin(t) cos(x2 ). 0 e−t We discretise in space with the Fourier pseudospectral method [10, 17] using 256 grid points per dimension. We compare to a reference solution, which is computed with ∆t = 2−10 . The results, which confirm the 4th order convergence, are displayed in Table 1. Table 1: Verification of the convergence rate for the Chin scheme (18). The method is confirmed to be 4th order accurate in L2 -norm for non-autonomous Schr¨ odinger-type and parabolic problems. Schr¨odinger Parabolic error rate error rate ∆t 2−1 1.364e-2 3.295e-4 2−2 7.855e-4 4.12 2.299e-5 3.84 2−3 3.287e-5 4.58 1.568e-6 3.87 2−4 1.853e-6 4.15 1.015e-7 3.95 2−5 1.130e-7 4.03 6.411e-9 3.99 2−6 7.022e-9 4.01 4.018e-10 4.00 2−7 4.381e-10 4.00 2.512e-11 4.00 2−8 2.728e-11 4.01 1.565e-12 4.00 We emphasise the independence of the norms of the operators with another numerical study. We consider the harmonic oscillator in two dimensions, 1 1 iut = − ∆u + (x21 + x22 )u, 2 2
(19)
and truncate the spatial domain to [−3π, 3π]. We start in the first eigenstate, 1
T
u(x, 0) = π −1/2 e− 2 x
x
,
and solve over one period, up to the time t = 2π, using 50 time steps. We use the Fourier pseudospectral method in space, and time step with the Chin scheme. The behaviour of the L2 -error, compared to the analytical solution, as the spatial step length ∆x1 = ∆x2 = 6π/N is decreased is shown in Table 2. When the solution is well resolved and the temporal error dominates, the error stays constant as the spatial grid is refined further. The quadratic growth with N of the discrete Laplace operator does not affect the error.
Acknowledgements Discussions with Vasile Gradinaru and Sverker Holmgren are gratefully acknowledged. 14
Table 2: Behaviour of the error in L2 for the harmonic oscillator (19) when the spatial grid is refined. We use ∆x1 = ∆x2 = 6π/N , and the constant time step ∆t = π/25. When the temporal error is dominant the error stagnates, the growing norm of the spatial operator does not influence the error. N error 8 2.481 16 1.261e-1 32 1.872e-5 64 3.630e-7 128 3.630e-7 256 3.630e-7 512 3.630e-7 1024 3.630e-7
References [1] P. Bader, S. Blanes, and F. Casas. Solving the Schr¨odinger eigenvalue problem by the imaginary time propagation technique using splitting methods with complex coefficients. J. Chem. Phys., 139:124117, 2013. doi:10.1063/1.4821126. [2] S. Blanes, F. Diele, C. Marangi, and S. Ragni. Splitting and composition methods for explicit time dependence in separable dynamical systems. J. Comput. Appl. Math., 235:646–659, 2010. doi:10.1016/j.cam.2010.06. 018. [3] S. Blanes and P. C. Moan. Practical symplectic partitioned Runge–Kutta and Runge–Kutta–Nystr¨ om methods. J. Comput. Appl. Math., 142:313– 330, 2002. doi:10.1016/S0377-0427(01)00492-7. [4] S. Blanes and M. Seydao˘ glu. High-order splitting methods for separable non-autonomous parabolic equations. Preprint, 2013. [5] F. Castella, P. Chartier, S. Descombes, and G. Vilmart. Splitting methods with complex times for parabolic equations. BIT, 49:487–508, 2009. doi: 10.1007/s10543-009-0235-y. [6] S. A. Chin. Symplectic integrators from composite operator factorizations. Phys. Lett. A, 226:344–348, 1997. doi:10.1016/S0375-9601(97)00003-0. [7] S. A. Chin. Structure of positive decompositions of exponential operators. Phys. Rev. E, 71:016703, 2005. doi:10.1103/PhysRevE.71.016703. [8] S. A. Chin and C. R. Chen. Fourth order gradient symplectic integrator methods for solving the time-dependent Schr¨odinger equation. J. Chem. Phys., 114:7338–7341, 2001. doi:10.1063/1.1362288. [9] S. A. Chin and C. R. Chen. Gradient symplectic algorithms for solving the Schr¨ odinger equation with time-dependent potentials. J. Chem. Phys., 117:1409–1415, 2002. doi:10.1063/1.1485725. 15
[10] B. Fornberg. A Practical Guide to Pseudospectral Methods. Cambridge University Press, Cambridge, 1996. [11] E. Hairer, C. Lubich, and G. Wanner. Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations. Springer, Berlin, 2002. [12] E. Hairer, S. P. Nørsett, and G. Wanner. Solving Ordinary Differential Equations I. Nonstiff Problems. Springer, Berlin, 2nd edition, 1993. [13] E. Hansen and A. Ostermann. Exponential splitting for unbounded operators. Math. Comput., 78:1485–1496, 2009. doi:10.1090/ S0025-5718-09-02213-3. [14] E. Hansen and A. Ostermann. High order splitting methods for analytic semigroups exist. BIT, 49:527–542, 2009. doi:10.1007/ s10543-009-0236-x. [15] T. Jahnke and C. Lubich. Error bounds for exponential operator splittings. BIT, 40:735–744, 2000. doi:10.1023/A:1022396519656. [16] O. Koch, C. Neuhauser, and M. Thalhammer. Error analysis of high-order splitting methods for nonlinear evolutionary Schr¨odinger equations and application to the MCTDHF equations in electron dynamics. M2AN Math. Model. Numer. Anal., 47:1265–1286, 2013. doi:10.1051/m2an/2013067. [17] H.-O. Kreiss and J. Oliger. Comparison of accurate methods for the integration of hyperbolic equations. Tellus, 24:199–215, 1972. doi: 10.1111/j.2153-3490.1972.tb01547.x. [18] I. P. Omelyan, I. M. Mryglod, and R. Folk. Construction of high-order forcegradient algorithms for integration of motion in classical and quantum systems. Phys. Rev. E, 66:026701, 2002. doi:10.1103/PhysRevE.66.026701. [19] G. Russo and P. Smereka. The Gaussian wave packet transform: Efficient computation of the semi-classical limit of the Schr¨odinger equation. Part 2 – multidimensional case. J. Comput. Phys., 257:1022–1038, 2014. doi: 10.1016/j.jcp.2013.09.023. [20] M. Suzuki. Hybrid exponential product formulas for unbounded operators with possible applications to Monte Carlo simulations. Phys. Lett. A, 201:425–428, 1995. doi:10.1016/0375-9601(95)00266-6. [21] M. Suzuki. New scheme of hybrid exponential product formulas with applications to quantum Monte-Carlo simulations. In D. P. Landau, K. K. Mon, and H.-B. Sch¨ uttler, editors, Computer Simulation Studies in CondensedMatter Physics VIII, volume 80 of Springer Proceedings in Physics, pages 169–174. Springer, Berlin, 1995. doi:10.1007/978-3-642-79991-4_21. [22] M. Thalhammer. High-order exponential operator splitting methods for time-dependent Schr¨odinger equations. SIAM J. Numer. Anal., 46:2022– 2038, 2008. doi:10.1137/060674636. 16
[23] M. Thalhammer. Convergence analysis of high-order time-splitting pseudospectral methods for nonlinear Schr¨odinger equations. SIAM J. Numer. Anal., 50:3231–3258, 2012. doi:10.1137/120866373.
17
Abstract We consider force-gradient, also called modified potential, operator splitting methods for problems with unbounded operators. We prove that force-gradient operator splitting schemes retain their classical orders of accuracy for time-dependent partial differential equations of parabolic or Schr¨ odinger type, provided that the solution is sufficiently regular.
1
Introduction
We consider the numerical time-integration of differential equations on the form ut = Au + Bu,
u(0) = u0 ,
(1)
evolving on a Banach space X, with the operators A : D(A) ⊂ X → X, B : D(B) ⊂ X → X. Splitting methods are advantageous if the subproblems with only one of the operators can be solved easily. An example of such a problem, where splitting methods often are used, is the time-dependent Schr¨odinger equation, iut = −∆u + V (x)u, u(x, 0) = u0 ∈ L2 (Rd ). If one of the terms is neglected we get the subproblems ivt = −∆v, iwt = V (x)w,
⇒
v(x, 0) = v0 ,
⇒
w(x, 0) = w0 ,
v(x, t) = eit∆ v0 , w(x, t) = e
−itV (x)
and w0 .
The operator exponentials are readily computed since the potential V (x) and the Laplacian ∆ are diagonal in coordinate and frequency space, respectively. The combined operator (−∆ + V (x)), however, is not as easily exponentiated. Splitting methods consist in taking small time steps, alternatingly applying the action of the two operators. An s-stage standard splitting method for (1) can be written as un+1 = Sun ,
S=
s Y
ebj hB eaj hA ,
(2)
j=1 ∗ Division of Scientific Computing, Department of Information Technology, Uppsala University, Box 337, 751 05 Uppsala, Sweden ([email protected]).
1
where h is the time step and the coefficients aj , bj define the splitting scheme. Throughout this paper products are defined in falling order, i.e., s Y
Aj = As As−1 . . . A1 ,
j=1
Ql and j=k aj = 1 if k > l. We let C denote a generic constant which may take different values at different occurrences. We say that a splitting method is of classical order p if the local error estimate k(S − eh(A+B) )uk ≤ Chp+1 holds for bounded operators A and B, and sufficiently small time steps h. Such estimates are obtained using Taylor expansion or the Baker–Campbell– Hausdorff formula [11, Ch. III.5]. They are adequate, for instance, for ordinary differential equations, with u ∈ RN , and A, B ∈ RN ×N . In infinite-dimensional settings, classical error analysis is no longer rigorous since the possibly unbounded operators A and B appear in the remainder. Still, splitting methods often retain their classical order also for such problems. Error estimates which hold independently of the norms of the operators are called stiff estimates. In an early work, Jahnke and Lubich proved that the Strang splitting scheme keeps its second order accuracy for Schr¨odinger equations [15]. Thalhammer [22], and Hansen and Ostermann [13], proved that splitting methods keep their classical order for certain classes of problems with unbounded operators, given that the solution is sufficiently regular. See also [5, 14], which enables splitting methods of order higher than two for parabolic problems using the theory of analytic semigroups. Using the calculus of Lie derivatives, the theory was extended to non-linear problems in [16, 23]. If the problem has a particular structure, this can often be used to construct more efficient methods. One example is Runge–Kutta–Nystr¨om methods [3, 11], which apply when commutators of the operators satisfy [B, [B, [A, B]]] = 0. This reduces the number of order conditions, and gives more freedom to the choice of method coefficients. In this paper we will study force-gradient methods, which were proposed by Suzuki [20, 21]. They rely on the same commutator assumption as Runge–Kutta–Nystr¨om methods, and in addition that [B, [A, B]] can be conveniently evaluated. A general force-gradient method for the problem (1) reads s Y 3 un+1 = Sun , S = ebj hB+cj h [B,[A,B]] eaj hA . (3) j=1
The perhaps most used force-gradient scheme is the fourth order method of Chin [6, 8, 9], with s = 3. An attractive property of the Chin scheme is that all its coefficients are real and non-negative. The method can therefore be used both for parabolic and Schr¨ odinger-type problems. Force-gradient methods of order higher than four have been derived, e.g., by Omelyan et al. [18]. These methods do however have some negative coefficients. This makes them unsuitable for parabolic problems, which cannot be solved backwards in time. Chin has proven that a 6th order force-gradient method necessarily has negative coefficients, unless higher order commutators are included in the scheme [7]. As for standard splitting methods, this order barrier can be circumvented if the coefficients are 2
allowed to be complex-valued. In [1], Bader et al. derive several high order forcegradient schemes with complex coefficients with positive real parts, suitable for parabolic problems. In this paper we derive stiff order conditions for force-gradient methods. That is, conditions under which the methods are of a certain order of accuracy also for unbounded operators. We then prove that the stiff order conditions are equivalent to the classical ones, as long as the exact solution is sufficiently regular. The result holds for parabolic and Schr¨odinger-type partial differential equations. Problems where force-gradient methods are applicable can be found in those classes, since [V, [∆, V ]] = −2|∇V |2 . The remainder of this paper is organised as follows. In Section 2 we specify the assumptions on the operators and the regularity requirements, and state local and global error estimates. The global estimate follows directly from the local estimate and the assumptions. The local error estimate, which is the main result of this paper, is proven in Section 3. In Section 4 we use the calculus of Lie derivatives to show how force-gradient methods can be used for nonautonomous problems. We conclude the paper with numerical experiments in Section 5, which corroborate the theoretical results.
2
Statement of the error estimate
In this section we list the assumptions we rely on and prove the global error estimate, Theorem 2. The global estimate relies on a local error estimate, Theorem 1, which we will prove in Section 3. The first assumption specifies the classes of problems under consideration. In particular, linear autonomous problems of parabolic and Schr¨ odinger type are included. Assumption 1. The problem (1) is assumed to be well-posed, and the commutators of A and B are assumed to satisfy [B, [B, [A, B]]] = 0. We further assume either of the following two conditions. 1. (Schr¨ odinger-type problem) The operators (A + B), A, B and [B, [A, B]] are the infinitesimal generators of strongly continuous groups of operators {etF }t∈R . Furthermore, there exists a constant ω ≥ 0 such that ketF kX←X ≤ eω|t|
∀ t ∈ R,
F = (A + B), A, B, and [B, [A, B]]. 2. (Parabolic problem) The operators (A + B), A, B, and [B, [A, B]] are the infinitesimal generators of analytic semigroups {etF }t∈Σθ , defined in the sectors Σθ = {t ∈ C : −θ ≤ arg t ≤ θ} for some 0 < θ ≤ π/2. Furthermore, there exists a constant ω ≥ 0 such that ketF kX←X ≤ eω|t|
∀ t ∈ Σθ ,
F = (A + B), A, B, and [B, [A, B]]. Next, we assume that the solution and the operators A and B are sufficiently regular. This assumption was justified for Schr¨odinger equations with sufficiently regular, bounded potential B in [22, Sec. 2]. The arguments transfer directly to parabolic problems on the form ut = ∆u + V (x)u. We here use the 3
notation |µ| = |µ1 | + · · · + |µk | for the modulus of a multi-index, and iterated commutators are denoted p adp+1 A (B) = [A, adA (B)],
ad0A (B) = B.
Assumption 2. For a force-gradient splitting method of classical order p, assume that k
Y
µ adAj (Bj )eζj A
≤ C,
X←Yp
j=1
µ ∈ Nk , |µ| = p + 1 − k,
Bj = B and [B, [A, B]],
j = 1, . . . , k,
(4) (5)
holds for k = 1, . . . , p + 1, with ζj ∈ R (Schr¨ odinger) or ζj ∈ Σθ (parabolic), |ζj | ≤ C, and some suitably chosen normed space Yp ⊂ X. Assume further that keζF kYp ←Yp ≤ C,
F = (A + B), A, B, and [B, [A, B]],
(6)
and that the exact solution to (1) resides in the space Yp . For a problem with smooth, bounded potential we get the standard Sobolev spaces Yp = H p . Under these assumptions we can prove the following local and global error estimates. Theorem 1. Assume that Assumptions 1 and 2 hold, and that the coefficients aj , bj , cj define a force-gradient operator splitting scheme of classical order p. Assume further that the coefficients reside in R if (1) is of Schr¨ odinger type, or in Σθ if parabolic. Then, the scheme satisfies the local error estimate kSu(t) − u(t + h)kX ≤ Chp+1 , also for unbounded operators A and B. Proof. We will prove this result in Section 3. Theorem 2. Assume the conditions of Theorem 1. Then, the force-gradient operator splitting scheme (3) satisfies the global error estimate kS n u0 − u(nh)kX ≤ Chp ,
0 ≤ nh ≤ T.
Proof. By Assumption 1, kSkX←X ≤ ech for some c. Using this in connection with the local error estimate, the result follows through a standard Lady Windermere’s fan argument [12, Ch. II.3]. See, e.g., [13, Thm. 2.3] for a proof of a much similar result.
3
Proof of the local error estimate
This section is devoted to the proof of Theorem 1. The line of attack is as follows. First, we summarise the proof of a similar result [22], a local error estimate for standard splitting methods on the form (2). This introduces notation and techniques which we will use when we extend the proof to force-gradient methods. Second, we take the commutators [B, [A, B]] into account, and manipulate 4
the expressions so that they fit in the framework of the stiff order conditions for standard methods. Finally, we observe that the stiff order conditions hold if and only if the classical order conditions do. Since regularity is assumed uniformly in time, we can without restriction consider the error in the first time step, e0 = (S − eh(A+B) )u0 . As a first step, we express the exact and approximate solutions in forms that facilitate comparison. Repeating calculations from [22], we express the exact solution after one time step as u(h) = ehA u0 +
p X
Z Ik u0 + Rp+1 ,
Ik =
gk (τ ) dτ ,
(7)
∆k
k=1
with τ = (τ1 , . . . , τk ), τ0 = h, and k Y
fk (τ ) =
e(τk−l −τk−l+1 )A B,
gk (τ ) = fk (τ )eτk A .
l=1
The integrals are evaluated over the domains ∆k = {τ ∈ Rk : 0 ≤ τk ≤ τk−1 ≤ · · · ≤ τ1 ≤ h}, and the remainder term reads Z Rp+1 =
fp+1 (τ )u(τp+1 ) dτ . ∆p+1
Throughout this section, we let R0p+1 and Rp+1 denote arbitrary linear combinations of terms on the form R0p+1 ∼ hp+1
k Y
µ adAj (Bj )eζj A χ,
Rp+1 = R0p+1 u,
|µ| ≤ p − k + 1,
j=1
with the linear operator χ bounded Yp → Yp . By Assumption 2, this means that kRp+1 kX ≤ Chp+1 . Note that Rp+1 = Rp+1 . Now, if we define Aj = aj hA,
Bj = bj hB,
Cj = cj h3 [B, [A, B]],
the approximate solution can be expressed similarly, S=
s Y j=1
e
Aj
+
p X
Qk + rp+1 .
k=1
Here, Qk =
X
κλ G(λ),
Λk = {λ ∈ Nk : 1 ≤ λk ≤ λk−1 ≤ · · · ≤ λ1 ≤ s},
λ∈Λk
F (λ) =
k Y
λk−l
Y
eAj (Bλk−l+1 + Cλk−l+1 ),
G(λ) = F (λ)
λk Y j=1
l=1 j=λk−l+1 +1
5
eAj ,
with λ = (λ1 , . . . , λk ) ∈ Nk and λ0 = s. The remainder term reads rp+1 =
p+1 X X
κ0j,λ F (λ)ϕj (Bλk + Cλk )eAλk
λY k −1
eBj +Cj eAj .
j=1
j=1 λ∈Λp+1
The function ϕj (z), which appear in the remainder, is defined recursively as ϕj (z) =
1 + zϕj+1 (z), j!
j ≥ 0,
ϕ0 (z) = ez ,
(8)
and κλ and κ0j,λ are well-behaved expansion coefficients. The expression ϕj (Bλk + Cλk )e
A λk
λY k −1
eBj +Cj eAj
j=1
is a bounded operator Yp → Yp due to (6), and therefore rp+1 = R0p+1 . The values of κ0λ are of lesser importance, but κλ appears in the order conditions. If σ0 = 0,
λσ0 +···+σl−1 +1 = · · · = λσ0 +···+σl ,
l = 1, . . . , m,
and σ0 + · · · + σm = k, then κλ =
1 . σ0 ! · · · σm !
The local error can now be expressed as p s Y X e0 = ehA − eAj u0 + (Ik − Qk )u0 + Rp+1 . j=1
(9)
k=1
In [22] it was shown how each of the terms in (9) could be bounded separately, yielding the desired result for a standard operator splitting scheme. For forcegradient methods this is not possible since some of the order conditions for standard splitting schemes are violated. Error cancellation between the terms will instead recover the desired order of accuracy. We start off ignoring this, and re-derive the order conditions for a standard method. After that, we investigate how corrections from the [B, [A, B]]-terms kick in. In what is coming, we need some notation for multi-indices. For µ ∈ Nk and x ∈ Ck , define µ! = µ1 !µ2 ! · · · µk ! and xµ = xµ1 1 xµ2 2 · · · xµk k . Furthermore, if 1 ≤ µj ≤ l and x ∈ Cl , let xµ = (xµ1 , . . . , xµk ). The first error term, s Y (0) e0 = ehA − eaj hA u0 , j=1
Ps is zero as long as the consistency condition j=1 aj = 1 is fulfilled. Taylor expansion of the integrands in (7) gives Z X 1 Z Ik = τ µ dτ ∂τµ gk (0) + %k,p−k+1 (τ ) dτ , (10) µ! ∆k ∆k |µ|≤p−k
6
with the remainder X
%k,n+1 =
ν∈Nk |ν|=n+1
n+1 ν!
Z
1
(1 − s)n ∂τν gk (sτ )τ ν ds,
0
and %k,n+1 = R0n+1 . The expansion coefficients in (10) can be evaluated as Z
τ µ dτ = hk+|µ|
k Y
∆k
l=1
1 , µl + · · · + µk + k − l + 1
and the partial derivatives of gk are ∂τµ gk (τ ) = (−1)|µ|
k Y
µ e(τk−l −τk−l+1 )A adAk−l+1 (B) eτk A .
(11)
l=1
Next, we let (0)
Qk = Qk |cj ≡0 , Pj and Taylor expand that. Defining αj = i=1 ai , we get (0) Qk
k
=h
X
κλ
k Y
λ∈Λk
=
X
l=1
h
k+|µ|
|µ|≤p−k
+ hk
bλl gk (αλ h)
X
k Y 1 X µ µ κλ bλl αλ ∂τ gk (0)+ µ! λ∈Λk
κλ
λ∈Λk
k Y
l=1
bλl %k,p−k+1 (αλ h).
l=1
Defining k
Πµ =
1 Y 1 , µ! µl + · · · + µk + k − l + 1 l=1
Σ(0) µ =
k Y 1 X µ κλ bλl αλ , µ! λ∈Λk
l=1
the kth order error term reads X (0) Ik − Qk = hk+|µ| Πµ − Σ(0) ∂τµ gk (0) + R0p+1 . µ |µ|≤p−k
We thereby achieve pth order of accuracy if αs = 1 and X µ Πµ − Σ(0) k = 1, . . . , p, σ = 0, . . . , p − k. µ ∂τ gk (0) = 0,
(12)
k
µ∈N |µ|=σ
These are the stiff order conditions for standard operator splitting methods derived in [22]. For force-gradient methods these conditions are typically violated, since the contributions from Cj are not yet taken into account. We will next extend the proof to address that. The key observation is that changing a B 7
(0)
for [B, [A, B]] in the expansion of Qk transforms gk into two first order partial derivatives of gk+1 . Take, for example, 3 Y
g2 (α(2,1) h) =
eaj hA B
j=3
2 Y
1 Y
eaj hA B
j=2
eaj hA ,
j=1
and compare it to 3 Y
=
eaj hA B
2 Y
j=3
j=2
3 Y
2 Y
eaj hA B
j=3
−
eaj hA B
1 Y
eaj hA
j=1 1 Y
eaj hA B[A, B]
j=2
3 Y j=3
=
eaj hA [B, [A, B]]
eaj hA +
j=1
2 Y
eaj hA [A, B]B
j=2
1 Y
eaj hA
j=1
−∂τ(0,0,1) g3 (α(2,1,1) h)
+
∂τ(0,1,0) g3 (α(2,1,1) h),
cf. (11). We introduce the operators Dl± as µ
Dl± ∂τµ gk (τ ) = ∂τ ±,l gk (τ 0 ), with µ ∈ Nk , τ ∈ Rk , and µ+,l = (µ1 , . . . , µl−1 , µl , µl + 1, µl+1 , . . . , µk ) ∈ Nk+1 , µ−,l = (µ1 , . . . , µl−1 , µl + 1, µl , µl+1 , . . . , µk ) ∈ Nk+1 , τ 0 = (τ1 , . . . , τl−1 , τl , τl , τl+1 , . . . , τk ) ∈ Rk+1 . With this notation we can write k
Qk = h
X λ∈Λk
cλ
k Y
+ − bλk−l+1 − h2 cλk−l+1 Dk−l+1 + h2 cλk−l+1 Dk−l+1 gk (αλ h).
l=1
(13) We expand the products in (13), and order the terms with respect to the number of D-operators in them. We let Qk =
k X
(m)
Qk ,
m=0 (0)
with Qk as before and (m)
Qk
X
= hk+2m
λ∈Λk
×
m Y
κλ
k X l1 =m
···
lm−1 −1 k X Y lm =1
l=1 l6=lj
− + − Dk−l cλk−lj +1 Dk−l j +1 j +1
j=1
8
bλk−l+1 ×
gk (αλ h),
m = 1, . . . k.
The point with this notation is that the error now can be expressed as
e0 = e
hA
−
s Y
k
e
Aj
b2c p X X (m) Ik − u0 + Qk−m u0 + Rp+1 .
j=1
m=0
k=1
The disregarded terms p X
k X
(m)
Qk u0
k=0 m=p−k+1
are of the form Rq , with q ≥ p + 1, and can thereby be controlled. With this formulation, we require no cross-cancellation between the terms k
b2c X
Ik −
(m) Qk−m .
m=0 (m)
In order to match terms with the same derivative of gk we Taylor expand Qk−m . (1)
For illustration, we first expand Qk−1 . (1)
X
Qk−1 = h(k−1)+2
k−1 X
k−1 Y
l1 =1
l=1 l6=l1
κλ
λ∈Λk−1
bλ(k−1)−l+1 cλ(k−1)−l1 +1 ×
− + gk−1 (αλ h) × D(k−1)−l − D(k−1)−l 1 +1 1 +1 X
= hk+1
κλ
λ∈Λk−1 e
k−1 X
k−1 Y
l1 =1
l=1 l6=l1
ek−l1 +1
× ∂τ k−l1 − ∂τ X
=
hk+1+|µ|
|µ|≤p−(k+1) µ+ek−l1
µ × αλ (l) ∂τ
bλk−l cλk−l1 ×
gk (αλ(l) h)
1 µ!
X
κλ
λ∈Λk−1
µ+ek−l1 +1
− ∂τ
k−1 X
k−1 Y
l1 =1
l=1 l6=l1
bλk−l cλk−l1 ×
gk (0) + R0p+1 ,
with λ(l) = (λ1 , . . . , λl−1 , λl , λl , λl+1 , . . . , λk−1 ) ∈ Nk . We rewrite this as (1)
Qk−1 =
X
µ 0 hk+|µ| Σ(1) µ ∂τ gk (0) + Rp+1 ,
with
1≤|µ|≤p−k
Σ(1) µ =
k−1 X
1 (µ − eν )! ν=1 k X
X
κλ
λ∈Λk−1
1 − (µ − eν )! ν=2
k−1 Y
µ−eν + bλk−l cλν αλ (l)
l=1 l6=k−ν
X
κλ
λ∈Λk−1
k−1 Y
µ−eν bλk−l cλν−1 αλ . (l)
l=1 l6=k−ν+1
Terms with µ − eν ∈ / Nk are excluded from the summation. Similarly in the 9
general case, (m)
X
Qk−m = h(k−m)+2m
κλ
λ∈Λk−m
×
m Y
k−m X
lm−1 −1 k−m
···
l1 =m
X
Y
lm =1
l=1 l6=lj
bλ(k−m)−l+1 ×
− + cλ(k−m)−lj +1 D(k−m)−l − D(k−m)−l j +1 j +1
gk−m (αλ h)
j=1
X
= hk+m
λ∈Λk−m
×
m Y
k−m X
κλ
···
lm−1 −1 k−m
l1 =m
X
Y
lm =1
l=1 l6=lj
e(k−m)−lj +j
cλ(k−m)−lj +1 ∂τ
bλ(k−m)−l+1 ×
e(k−m)−lj +j+1
− ∂τ
gk (αλ(l) h)
j=1
X
=
k+m+|µ|
h
|µ|≤p−(k+m)
×
m Y
1 µ!
X
κλ
λ∈Λk−m
e(k−m)−lj +j
cλ(k−m)−lj +1 ∂τ
k−m X
···
lm−1 −1 k−m
l1 =m
X
Y
lm =1
l=1 l6=lj
e(k−m)−lj +j+1
− ∂τ
bλ(k−m)−l+1 ×
µ µ αλ (l) ∂τ gk (0)+
j=1
+ R0p+1 . λ(l) is defined analogously to λ(l) , with all the elements λlj , j = 1, . . . , m, doubled. In order to turn this expression inside out we need to introduce some notation. Let B = {0, 1}, and define the function ν : Nm × Bm → Nm as νj (l, v) = k − m − lj + j + vj ,
j = 1, . . . , m,
l = (l1 , . . . , lm ).
Then the product of binomials above can be expanded as m Y
e(k−m)−lj +j
∂τ
e(k−m)−lj +j+1
− ∂τ
=
X
(−1)|v| ∂τν(l,v) .
v∈Bm
j=1 (m) Qk−m
We can now rewrite the expression for as X (m) µ 0 Qk−m = hk+|µ| Σ(m) µ ∂τ gk (0) + Rp+1 ,
with
µ∈Nk m≤|µ|≤p−k
Σ(m) µ
=
k−m X
lm−1 −1
X
···
lm =1 v∈Bm
l1 =m
×
X
k−m Y
bλ(k−m)−l+1
(−1)|v| (µ − ν(l, v))!
m Y
X
κλ ×
λ∈Λk−m
µ−ν(l,v) cλ(k−m)−lj +1 αλ(l) .
j=1
l=1 l6=lj
As before, terms with µ − ν(l, v) ∈ / Nk are omitted from the summation. Assuming the consistency condition αs = 1, we can write the local error as e0 =
p X k=1
k
X k
µ∈N |µ|≤p−k
k+|µ|
h
Πµ −
b2c X m=0
10
Σ(m) ∂τµ gk (0)u0 + Rp+1 . µ
(14)
Noting that ∂τµ gk (0) = (−1)|µ| ehA
k Y
µ
adAk−l+1 (B),
l=1
we can write the stiff order conditions for force-gradient methods as1 X k
µ∈N |µ|=σ
Eµ
k Y
k
µ adAk−l+1 (B)
= 0,
Eµ = Πµ −
b2c X
Σ(m) µ ,
(15)
m=0
l=1
for all k = 1, . . . , p, σ = 1, . . . , p − k. Note the similarity with the order conditions (12) for standard splitting methods. Note also that terms in (15) can be matched using the properties of commutators, in particular the Jacobi identity and the assumed relation [B, [B, [A, B]]] = 0 of Runge–Kutta–Nystr¨om methods. Finally, by repeating the arguments from [22], one can show that the conditions (15) are equivalent to the classical order conditions. We assume that A and B are bounded operators, e.g., square matrices of fixed size, and claim that (15) are sufficient and necessary conditions for pth order convergence. They are sufficient conditions for bounded operators, since they are for unbounded. By Taylor expanding ehA in the partial derivatives of gk (0) in (14) and collecting terms of similar order in h, we see that the conditions are necessary also for non-stiff problems.
4
Generalisation to non-autonomous problems
Using the calculus of Lie derivatives [11, Ch. III.5], the results presented in this paper can be generalised to non-linear problems. That is, the proofs of the theorems are readily repeated, but computing the flows of the subproblems might not always be practical. The requirement that the commutator [B, [A, B]] should be easy to handle is also restrictive. The practical use of force-gradient methods for non-linear problems is therefore limited. The non-linear framework does, however, provide a convenient way of generalising the results to nonautonomous linear problems. That is the goal of this section. Given two non-linear operators A, B : D ⊂ X → X, the Lie derivative DA of A and its exponential are defined by DA Bu = B 0 (u)A(u),
t1 ,t2 e(t2 −t1 )DA Bu = B(ΨA u).
ΨtA1 ,t2 is here the exact flow of the differential equation ut = A(u) from t = t1 to t = t2 , and B 0 (u) denotes the Fr´echet derivative of B(u). Taking B as the identity mapping yields DA u = A(u),
t1 ,t2 e(t2 −t1 )DA u = ΨA u.
The Lie commutator [DA , DB ] will be an essential tool in the following. It is defined through [DA , DB ]u = DA DB u − DB DA u = B 0 (u)A(u) − A0 (u)B(u). 1A
Matlab script for evaluating the expressions Eµ is available upon request.
11
To generalise Theorems 1 and 2 to non-linear problems on the form ut = A(u) + B(u),
u(0) = u0 ,
we only need to assert that the problem satisfies a modified Assumption 2. In the modified assumption, (4)–(5) are replaced by k
Y
µ adDjA (DBj )eζj DA
≤ C,
X←Yp
j=1
µ ∈ Nk , |µ| = p + 1 − k,
DBj = DB or [DB , [DA , DB ]],
j = 1, . . . , k,
for k = 1, . . . , p + 1. Provided that this holds, the proofs can be carried out as previously, only with the order of the exponentials reversed. We will now show that this assumption holds for non-autonomous linear problems ut = A(t)u + B(t)u, u(0) = u0 , satisfying [A(s), A(t)] = [B(s), B(t)] = 0 for all s, t. We will furthermore show that the Lie commutator [DB , [DA , DB ]] has a simple form and can be used as before. We introduce time as an auxiliary variable and write the problem as u Ut = F (U ) + G(U ), with U = , s and
F (U ) =
A(s)u , 1
G(U ) =
B(s)u . 0
(16)
The Fr´echet derivatives of F and G are then A(s) A0 (s)u B(s) B 0 (s)u F 0 (U ) = , G0 (U ) = , 0 0 0 0 and the iterated Lie commutators can be determined as ! p f ad (B(s))u p A(s) adDF (DG )U = , 0 with p+1 p d fp f f ad (adA(s) (B(s))), A(s) (B(s)) = [adA(s) (B(s)), A(s)] + ds 0 f ad (B(s)) = B(s). A(s)
We can conclude that in order to fulfil the modified Assumption 2 we, in addition to the conditions from the autonomous case, must require the operators A and B to be temporally smooth. This is a reasonable requirement. The force-gradient correction similarly evaluates to dB(s) −[B(s), [B(s), A(s)]]u + dB(s) ds B(s)u − B(s) ds u [DG , [DF , DG ]]U = 0 [B(s), [A(s), B(s)]]u = . (17) 0 12
Since B(s) commutes with itself evaluated at any other time, it also commutes with its time derivative. The commutator in (17) is the same commutator as in the autonomous problem, evaluated at a single point in time. Furthermore, since the operator A evaluated at different times commutes, its exponentiation is given by the first term of its Magnus series, Z t2 t1 ,t2 ΨA = exp A(σ) dσ . t1
In a computational context, it is sufficient to approximate the integral in the exponent with a pth order quadrature rule. Using the recurrence relation (8) of the functions ϕj the quadrature remainder can be isolated from the scheme. As an example, consider Strang splitting with the integral of A(t) approximated with the midpoint rule. Writing the quadrature rule as the exact integral and a truncation error, we get Z tn +h ehB/2 exp A(σ) dσ + Ch3 A00 (ξ) ehB/2 = tn
e
hB/2
exp
Z
tn +h
A(σ) dσ ehB/2 +
tn
+ Ch3 ehB/2 exp
Z
tn +h
A(σ) dσ A00 (ξ)ϕ1 Ch3 A00 (ξ) ehB/2 ,
tn
with ξ ∈ [tn , tn + h]. Note that it is possible to introduce time as an auxiliary variable in several ways [2]. One may even prefer to use more than one auxiliary variable. In connection with force-gradient methods, however, the arrangement (16) seems the most practical. The reason for this is that we want time to “stand still” while G(U ) is integrated, otherwise the commutator (17) will become complicated. When using schemes with complex coefficients for non-autonomous problems, one may have to evaluate the operators at complex times. This can lead to poor conditioning. In order to circumvent this difficulty, Blanes and Seydao˘ glu suggest methods with aj ∈ R+ , bj ∈ Σθ [4]. No force-gradient method with this property (and cj ∈ Σθ ) of order higher than four has to our knowledge been proposed, but could in principle be derived.
5
Numerical experiments
Preservation of the classical order for force-gradient schemes applied to partial differential equations have been observed previously. The Chin scheme, with s = 3 and 1 2 1 1 , a1 = 0, a2 = a3 = , b1 = b3 = , b2 = , c1 = c3 = 0, c2 = 2 6 3 72 (18) was studied for the Schr¨ odinger equation in [8, 9]. In [19], fourth order accuracy was observed for a non-autonomous Schr¨odinger-like problem in four dimensions. For completeness, we present convergence studies for the Chin scheme applied to the two-dimensional Schr¨ odinger-type and parabolic problems 1 iut = − ∇ · A(t)∇u + V (x, t)u, 2 13
ut = ∇ · A(t)∇u + V (x, t)u,
on x ∈ [−π, π]2 . We use periodic boundary conditions and the Gaussian initial condition T u(x, 0) = e−5 x x . The matrix A and the potential V are given by t e 0 A(t) = , V (x, t) = 2 − cos(x1 ) − sin(t) cos(x2 ). 0 e−t We discretise in space with the Fourier pseudospectral method [10, 17] using 256 grid points per dimension. We compare to a reference solution, which is computed with ∆t = 2−10 . The results, which confirm the 4th order convergence, are displayed in Table 1. Table 1: Verification of the convergence rate for the Chin scheme (18). The method is confirmed to be 4th order accurate in L2 -norm for non-autonomous Schr¨ odinger-type and parabolic problems. Schr¨odinger Parabolic error rate error rate ∆t 2−1 1.364e-2 3.295e-4 2−2 7.855e-4 4.12 2.299e-5 3.84 2−3 3.287e-5 4.58 1.568e-6 3.87 2−4 1.853e-6 4.15 1.015e-7 3.95 2−5 1.130e-7 4.03 6.411e-9 3.99 2−6 7.022e-9 4.01 4.018e-10 4.00 2−7 4.381e-10 4.00 2.512e-11 4.00 2−8 2.728e-11 4.01 1.565e-12 4.00 We emphasise the independence of the norms of the operators with another numerical study. We consider the harmonic oscillator in two dimensions, 1 1 iut = − ∆u + (x21 + x22 )u, 2 2
(19)
and truncate the spatial domain to [−3π, 3π]. We start in the first eigenstate, 1
T
u(x, 0) = π −1/2 e− 2 x
x
,
and solve over one period, up to the time t = 2π, using 50 time steps. We use the Fourier pseudospectral method in space, and time step with the Chin scheme. The behaviour of the L2 -error, compared to the analytical solution, as the spatial step length ∆x1 = ∆x2 = 6π/N is decreased is shown in Table 2. When the solution is well resolved and the temporal error dominates, the error stays constant as the spatial grid is refined further. The quadratic growth with N of the discrete Laplace operator does not affect the error.
Acknowledgements Discussions with Vasile Gradinaru and Sverker Holmgren are gratefully acknowledged. 14
Table 2: Behaviour of the error in L2 for the harmonic oscillator (19) when the spatial grid is refined. We use ∆x1 = ∆x2 = 6π/N , and the constant time step ∆t = π/25. When the temporal error is dominant the error stagnates, the growing norm of the spatial operator does not influence the error. N error 8 2.481 16 1.261e-1 32 1.872e-5 64 3.630e-7 128 3.630e-7 256 3.630e-7 512 3.630e-7 1024 3.630e-7
References [1] P. Bader, S. Blanes, and F. Casas. Solving the Schr¨odinger eigenvalue problem by the imaginary time propagation technique using splitting methods with complex coefficients. J. Chem. Phys., 139:124117, 2013. doi:10.1063/1.4821126. [2] S. Blanes, F. Diele, C. Marangi, and S. Ragni. Splitting and composition methods for explicit time dependence in separable dynamical systems. J. Comput. Appl. Math., 235:646–659, 2010. doi:10.1016/j.cam.2010.06. 018. [3] S. Blanes and P. C. Moan. Practical symplectic partitioned Runge–Kutta and Runge–Kutta–Nystr¨ om methods. J. Comput. Appl. Math., 142:313– 330, 2002. doi:10.1016/S0377-0427(01)00492-7. [4] S. Blanes and M. Seydao˘ glu. High-order splitting methods for separable non-autonomous parabolic equations. Preprint, 2013. [5] F. Castella, P. Chartier, S. Descombes, and G. Vilmart. Splitting methods with complex times for parabolic equations. BIT, 49:487–508, 2009. doi: 10.1007/s10543-009-0235-y. [6] S. A. Chin. Symplectic integrators from composite operator factorizations. Phys. Lett. A, 226:344–348, 1997. doi:10.1016/S0375-9601(97)00003-0. [7] S. A. Chin. Structure of positive decompositions of exponential operators. Phys. Rev. E, 71:016703, 2005. doi:10.1103/PhysRevE.71.016703. [8] S. A. Chin and C. R. Chen. Fourth order gradient symplectic integrator methods for solving the time-dependent Schr¨odinger equation. J. Chem. Phys., 114:7338–7341, 2001. doi:10.1063/1.1362288. [9] S. A. Chin and C. R. Chen. Gradient symplectic algorithms for solving the Schr¨ odinger equation with time-dependent potentials. J. Chem. Phys., 117:1409–1415, 2002. doi:10.1063/1.1485725. 15
[10] B. Fornberg. A Practical Guide to Pseudospectral Methods. Cambridge University Press, Cambridge, 1996. [11] E. Hairer, C. Lubich, and G. Wanner. Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations. Springer, Berlin, 2002. [12] E. Hairer, S. P. Nørsett, and G. Wanner. Solving Ordinary Differential Equations I. Nonstiff Problems. Springer, Berlin, 2nd edition, 1993. [13] E. Hansen and A. Ostermann. Exponential splitting for unbounded operators. Math. Comput., 78:1485–1496, 2009. doi:10.1090/ S0025-5718-09-02213-3. [14] E. Hansen and A. Ostermann. High order splitting methods for analytic semigroups exist. BIT, 49:527–542, 2009. doi:10.1007/ s10543-009-0236-x. [15] T. Jahnke and C. Lubich. Error bounds for exponential operator splittings. BIT, 40:735–744, 2000. doi:10.1023/A:1022396519656. [16] O. Koch, C. Neuhauser, and M. Thalhammer. Error analysis of high-order splitting methods for nonlinear evolutionary Schr¨odinger equations and application to the MCTDHF equations in electron dynamics. M2AN Math. Model. Numer. Anal., 47:1265–1286, 2013. doi:10.1051/m2an/2013067. [17] H.-O. Kreiss and J. Oliger. Comparison of accurate methods for the integration of hyperbolic equations. Tellus, 24:199–215, 1972. doi: 10.1111/j.2153-3490.1972.tb01547.x. [18] I. P. Omelyan, I. M. Mryglod, and R. Folk. Construction of high-order forcegradient algorithms for integration of motion in classical and quantum systems. Phys. Rev. E, 66:026701, 2002. doi:10.1103/PhysRevE.66.026701. [19] G. Russo and P. Smereka. The Gaussian wave packet transform: Efficient computation of the semi-classical limit of the Schr¨odinger equation. Part 2 – multidimensional case. J. Comput. Phys., 257:1022–1038, 2014. doi: 10.1016/j.jcp.2013.09.023. [20] M. Suzuki. Hybrid exponential product formulas for unbounded operators with possible applications to Monte Carlo simulations. Phys. Lett. A, 201:425–428, 1995. doi:10.1016/0375-9601(95)00266-6. [21] M. Suzuki. New scheme of hybrid exponential product formulas with applications to quantum Monte-Carlo simulations. In D. P. Landau, K. K. Mon, and H.-B. Sch¨ uttler, editors, Computer Simulation Studies in CondensedMatter Physics VIII, volume 80 of Springer Proceedings in Physics, pages 169–174. Springer, Berlin, 1995. doi:10.1007/978-3-642-79991-4_21. [22] M. Thalhammer. High-order exponential operator splitting methods for time-dependent Schr¨odinger equations. SIAM J. Numer. Anal., 46:2022– 2038, 2008. doi:10.1137/060674636. 16
[23] M. Thalhammer. Convergence analysis of high-order time-splitting pseudospectral methods for nonlinear Schr¨odinger equations. SIAM J. Numer. Anal., 50:3231–3258, 2012. doi:10.1137/120866373.
17
