Error analysis of the combined discretization methods

Now we sketch the second problem, namely, the error analysis of the combined method.

The naive approach suggests two “principles”:

a. If the operator splitting is exact (the splitting error vanishes), then a given numerical method with the same step-size gives the same result without and with splitting;

b. The “principle of the weakest chain link” is valid, i.e., when we apply an operator splitting method and a numerical method of given orders, then the order of the combined method is determined by the lower order.

We show that both conjectures are false.

First we examine the question a. by illustrating the interaction between the splitting procedure and the numerical method on a simple problem applying thesequential splitting and theexplicit Euler method.

Let (x(t), y(t))^T denote a function of type IR → IR², and 0 ≤ t ≤ t^?. The time-evolution of theharmonic oscillator is described by the following equation (see, e.g., [84]):

˙ wherex(t) andy(t) represent the amplitude and the velocity of the oscillator, respectively.

The simplest example of an oscillating system is a mass connected to a rigid foundation with a spring.

Equation (3.5.18) can be written as:

d

We use the decomposition (splitting) of the matrix appearing in (3.5.20)⁷ : µ 0 1

Since the first matrix is the identity matrix, clearly the commutator is zero, i.e., there is no splitting error. If we apply the sequential splitting, the split solution in the points kτ (k = 0,1, . . . , m) with τ := ^t_m^? (m∈IN) are

7The stability of the sequential splitting for this decomposition is proven in [90]. For the more general equation with function coefficients the stability is considered in the recent paper [14].

Let us apply the explicit Euler method with step sizeτfor (3.5.18). Then thenumerical

If the explicit Euler method is applied to the split problems, thenumerical split solution has the following form:

Since the split and the exact solutions do not differ, we expect that if we use the same numerical method for both the split and the unsplit problems, then the numerical solution and the numerical split solution do not differ, either. However, by comparing (3.5.23) and (3.5.24) it can be seen that a certain error appears in the numerical split solution. Hence, we can see that there is a certain interaction between the splitting procedure and the numerical method.

In our numerical experiments we used τ = ₂₀₀^2π while integrating system (3.5.18) on the time interval [0,4π] with and without applying splitting (the period of the harmonic oscillator is 2π). In the left panel of Figure 3.5.1 the effect of the numerical error can be seen: the numerical solution is spiraling outwards compared to the exact solution.

In the right panel of Figure 3.5.1 the effect of the interaction of the errors can be seen:

the numerical split solution is spiraling inwards compared to the exact solution (which coincides with the split solution in this case). Hence, in this case the interaction between the numerical method and the splitting procedure causes the “turn” of the spiral.

-1.5 Effect of the numerical error

exact solution Effect of the interaction error

exact and split solution numerical split solution

Figure 3.5.1: Effects of the numerical error (left panel) and interaction error (right panel).

To show that the conjecture b. is not true, it is enough to refer to the analysis of the Crank-Nicolson method, written by splitting+ numerical method in Section 3.5.1 on p.124: first order numerical methods were applied, however, the combined method had second order accuracy.

The above facts show how complex the error mechanism is in the combination “op-erator splitting + numerical methods to the sub-problems”. Next, we give the way how to choose an appropriate numerical method and its discretization parameter for a given operator splitting.

Theorem 3.5.2 Let us apply an operator splitting procedure of order p < ∞ together with a numerical method of order r, and set ∆t =τ^s (s ≥1, τ < 1). Then the order of the local total error is: %= min{p, rs}.⁸

This means that in order to preserve the accuracy, we can control it through the parame-ters p,r, ands for a givenp. In practice, we apply a given splitting procedure, i.e., pand τ are fixed, and we want to preserve the order p to the local total error (% =p). Hence, in this case we can estimate it, as well.

From Theorem 3.5.2 it follows that the interaction error causes a reduction in the order of Etot(τ) unless rs≥p.

Corollary 3.5.3 From Theorem 3.5.2 it follows that % =p if rs≥p.

Therefore, we shall answer the following two questions for given values of pand τ.

1. How to choose ∆t for a given numerical method of orderr?

2. Fixing ∆t =τ, how to choose the order r of the numerical method?

From Theorem 3.5.2, the following answers can be stated.

Theorem 3.5.4 When a given splitting procedure of order p is applied together with a given numerical method of order r ≤ p, and the numerical step size is set as ∆t = τ^s, then the exponent s has to be chosen as s = ^p_r in order to keep % = p. For r > p the identity %=p holds independently of the choice of ∆t.

Remark 3.5.5 From an algorithmical point of view, it will be much easier to select the case r > p. For this case, clearly, the choice ∆t = τ is optimal, because in this case the integration of the model needs the least computational work.

Theorem 3.5.6 When a given splitting procedure of order p is applied together with a certain numerical method of orderr, and the numerical step size is set as∆t =τ^s (s≥1), then r has to be chosen as r =£_p

s

¤+ 1 ∈IN in order to keep %=p.

Remark 3.5.7 Higher-order numerical methods could be chosen as well, but it would not lead to a higher-order total time discretization. It would only need more computational work.

3.5.3 Richardson-extrapolated sequential splitting with