3.4 Further investigations of the operator splittings
3.4.2 Local error analysis
In Section 3.2 we gave the conditions under which the local splitting error vanishes for the sequential splitting. As it was proven, for two operators the commutativity was a necessary and sufficient condition. In other words, the operator norm of the commutator should equal zero to get zero local splitting error. In the following we analyze whether the dependence of the magnitude of the local splitting error is continuous on the norm of the commutator. (Is it true that to a small commutator norm corresponds a small local splitting error?)
When the order of the local splitting error is p, then this means the following: for sufficiently small values ofτ the magnitude of Errspl(τ) is defined by the leading term of the error, i.e., by the coefficient ofτp+1. We observed that for the operator splittings under investigation this coefficient can be expressed by the commutators of the operators. This implies the direct dependence of the splitting error on the commutator norm. However, we should emphasize that this is true only for a sufficiently small τ. The following example shows that the magnitude of the norm of the commutator does not define the magnitude of the local splitting error in any case.
Example 3.4.4 Assume that
τ 12τ2k[A2, A1]ek∞ 1
6τ3kRAek∞
1 1.25(−1) 6.25(−2)
0.1 1.25(−3) 6.25(−5)
0.01 1.25(−5) 6.25(−8)
0.001 1.25(−7) 6.25(−11)
Table 3.4.1: The second and third terms in the splitting error for A.
A=A1+A2, p∈IR(p6= 0) arbitrary ,B1 =pA1, B2 = 1
pA2, B =B1+B2, v0 = (1,1,1,1) and τ = 1. We consider the following problems:
u0(t) =Au(t), t∈(0,1]
u(0) =v0
¾
(3.4.25)
and w0(t) = Bw(t), t∈(0,1]
w(0) =v0
¾
. (3.4.26)
Then
k[A1, A2]k∞ =k[B1, B2]k∞= 1 4, that is the norm of the two different commutators are equal.
Let us select p = 1000, and apply the sequential splitting for both problems (3.4.25) and (3.4.26). Then for the local splitting error we get
°°ErrAseq(1)°
°∞ =°
°(e(A1+A2)−eA2eA1)v0°
°∞= 0.125
and °
°ErrseqB (1)°
°∞=°
°(e(B1+B2)−eB2eB1)v0°
°∞= 20.8334.
One can see that in the second example we got a much (166.67 times) bigger local splitting error than in the first example. What is more, choosing bigger values of p, we can obtain arbitrarily big differences. This example suggests that the size of the splitting error is not completely determined by the commutator norm.
The reason is that in this example the valueτ = 1 is still quite big, and the contribution of the higher-order terms in the local splitting error is still important. When we decrease the value of τ, then the difference between the local splitting errors disappears. In order to analyze this phenomenon, we compute the coefficient at τ3 for the error. The direct calculation shows that
Errseq(τ) = τ2
2 [A2, A1]w0+ τ3
6RAw0+O(τ4), (3.4.27) where the operator RA has the form
RA=A1A2A1+A2A1A2+A21A2+A1A22−2A22A1−2A2A21. (3.4.28) (The same expression is valid for B.) In our example, for the two different splittings, the second and the third order terms in the local splitting error for the operators A and B are shown in Table 3.4.1 and Table 3.4.2, respectively.
The total splitting errors for the different values of τ are shown in Table 3.4.3.
τ 12τ2k[B2, B1]ek∞ 1
6τ3kRBek∞
1 1.25(−1) 20.833
0.1 1.25(−3) 2.08(−2)
0.01 1.25(−5) 2.08(−5)
0.001 1.25(−7) 2.08(−8)
Table 3.4.2: The second and third terms in the splitting error for B.
τ splitting error for A splitting error for B
1 1.25(−1) 20.833
0.1 1.25(−3) 2.08(−2)
0.01 1.25(−5) 2.08(−5)
0.001 1.25(−7) 1.25(−7)
Table 3.4.3: The total splitting errors for A and B.
We can see that the local splitting error for the operator A is determined by the second order term for each τ; however, for the splitting of B it does only for the last value (τ = 0.001), for the bigger values the third order term defines the error.
The above example raises the question of a suitable choice of the splitting time dis-cretization parameter τ: when it is too big, we lose the order; otherwise, when it is too small, we have to execute too many computations in order to get an approximation at some fixed time level. One can get a good hint by using the magnitude ofkRAk. However, formula (3.4.28) shows that the expressions for the coefficients at τk become more and more complicated with increasing k. To avoid this problem, i.e., to write operator RA in a more compact form, we can use the Baker-Campbell-Hausdorff (BCH) formula [143], which is very useful for the sequential splitting with two operators. The idea is to write the product of two matrix exponentials as a matrix exponential of one matrix.6 We seek such Cn =Cn(A2, A1), n = 0,1, . . . for which the relation
eτ A2eτ A1 =eP∞n=0τnCn (3.4.31)
6The Baker-Campbell-Hausdorff formula is the solution to Z = log(eXeY) for non-commuting X and Y. It was first noted in print by Campbell [21], elaborated by Poincar´e [109] and Baker [6], and systematized by Hausdorff [66]. The formula below was introduced by Dynkin. For an arbitrary Lie algebra the formula is cumbersome and in practice it is not perspicuous. However, for a matrix Lie group we obtain a simpler formula:
Z=X
n>0
(−1)n−1 n
X
ri+si>0 1≤i≤n
Xr1Ys1· · ·XrnYsn
r1!s1!· · ·rn!sn! . (3.4.29)
The Zassenhaus formula [155] is
et(X+Y)=etX etY e−t22[X,Y] et63(2[Y,[X,Y]]+[X,[X,Y]]) et4.... . . (3.4.30) which shows the necessity and sufficiency of commutatitivity of the operators for the relationetX etY = etY etX.For some more details we refer to [60].
holds. It is not difficult to show that, for the first coefficients, we have C0 =C0(A2, A1) = 0, C1 =C1(A2, A1) =A1+A2,
C2 =C2(A2, A1) = 12[A2, A1], C3 =C3(A2, A1) = 16[A2−A1,12[A2, A1]].
(3.4.32) For the further valuesn(n= 4,5. . .) for the matricesCn=Cn(A2, A1) there is a recursion formula, see (3.4.29) in the footnote. We remark that from this formula follows also the fact that the sequential splitting error for two operators vanishes if and only if the operators commute. According to the definition of the exponential of operators, and using the formulae (3.4.32) in (3.4.31), we obtain
eτ(A1+A2)=I+τ(A1+A2) + 1
2!τ2(A1+A2)2 + 1
3!τ3(A1+A2)3+O(τ4), (3.4.33) and
eτ A2eτ A1 =eτ(A1+A2)+P∞n=2τnCn =I +τ(A1+A2) + X∞
n=2
τnCn+
+1
2![τ(A1+A2) + X∞
n=2
τnCn]2+ 1
3![τ(A1 +A2) + X∞
n=2
τnCn]3+O(τ4). (3.4.34) It is easy to see that expressions (3.4.33) and (3.4.34) are equal up to the first order in τ.
Moreover, the coefficient of the second order term in the difference (3.4.34)-(3.4.33) reads 1
2![A2, A1]. (3.4.35)
The resulting difference operator at τ3 can be written as RA =C3+ 1
2!(A1+A2)C2+ 1
2!C2(A1+A2).
Applying the corresponding expressions forC2 and C3, we get RA= 1
12[A2−A1,[A2, A1]] + 1
4(A2+A1) [A2, A1] + 1
4[A2, A1] (A2+A1). (3.4.36) Using the trivial upper bound of the norm for the commutator, we get
kRAk ≤(1
6kA2−A1k∞+1
2kAk∞)k[A2, A1]k∞. (3.4.37) From the above formula one can draw the following conclusions. In any non-trivial case, the estimating expression vanishes only if the operators commute. However, if the norm of the original operator or the norm of the difference of the split operators are big, then the obtained upper bound can also be big. Therefore the error may be (but not necessarily is) significant, even if the commutator norm is relatively small. This is in accordance with the results of our numerical experiment, since
kAk∞= 1, kBk∞= 500.0005, kA2 −A1k∞= 1
2, kB2−B1k∞ = 749.99925, so kB2−B1k∞ >>kA2−A1k∞ and kBk∞>> kAk∞.
The dominance of the second order term, by using the estimation (3.4.37) means the condition
τ ≤ 3
kA2−A1k∞+ 3kAk∞. (3.4.38) In our example, for the splitting of operator B, this gives the bound τ ≤ 1.33·10−3. However, we can give another estimation, since, RA can be rewritten as
RA = 1 6
©(2A2+A1) [A2, A1] + [A2, A1] (2A1+A2)ª
. (3.4.39)
Hence, the condition (3.4.38) can be modified as
τ ≤ 3
2©
k2A2+A1k∞+k2A1+A2k∞ª. (3.4.40) In our example, for the splitting of operator B, this gives the bound τ ≤1.0·10−3, which is sharper (c.f. numerical results in Table 3.4.3).