# Iterative Splitting Methods for Integrodifferential Equations: Theory and Applications

Authors: Jürgen Geiser
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We present novel iterative splitting methods to solve integrodifferential equations. Such integrodifferential equations are applied, for example, in scattering problems of plasma simulations. We concentrate on a linearised integral part and a reformulation to a system of first order differential equations. Such modifications allow for applying standard iterative splitting schemes and for extending the schemes, respecting the integral operator. A numerical analysis is presented of the system of semidiscretised differential equations as abstract Cauchy problems.
In the applications, we present benchmark and initial realistic applications
to transport problems with scattering terms. We also discuss the benefits of such iterative schemes as fast solver methods.


---

## Body

## 1. Introduction

The motivation is to apply fast splitting schemes that have been developed for systems of differential equations, to integrodifferential equations. The applications are related to transport problems with scattering terms, for example, plasma simulations with collision terms; see [1, 2]. Due to the physics of the situation, the numerical implementation allows decoupling the system into transport and scattering terms, so that the underlying methods can employ optimal splitting schemes.We concentrate on a semidiscretised transport-scattering equation, given as a Cauchy problem. In the following, we consider the abstract homogeneous Cauchy problem as a system of ordinary integrodifferential equations:(1)du(t)dt=A(u(t))+∫0t‍B(u(s))ds,t∈[0,T],(2)u(0)=u0,
where A,B∈X→X are bounded and nonlinear operators, which are derived from the semidiscretised transport and scattering operators. X∈Rn is a Banach space and ∥·∥ is the corresponding norm in X. The operator A represents the transport processes, and the operator B is a given scattering operator.This paper is organised as follows.In Section2 we present the mathematical model and a possible reduced model for further approximations. The numerical schemes are presented in Section 3. The results of some numerical experiments are given in Section 4. The contents of Section 5 contain a summary of the results.


## 2. Analytical Approach to the Integrodifferential Equation

In the following, we discuss the transformation to a system of differential equations in order to apply iterative splitting schemes; see also the decomposition of matrices in [3].For the linearisation, we make the following assumptions.Assumption 1.
We assume a weak nonlinearity and employ the following simplification.(1)
The Jacobian ofA, given as dA/du, is linearised at t=0, which means (dA/du)|t=0=A, where A∈X×X.(2)
The operatorB is linearised at t=0, which means B(u)|t=0=B, where B∈X×X.(3)
The initial condition of the first derivative is given bydu0/dt=Au0.(4)
For the linearised case, the following matrix operations can be applied:(3)(AAt4+B)1/2=V(D14+D2)V-1,
where AAt and B are simultaneous and diagonalisable, with the same matrix V, to the diagonal matrices D1 and D2. Furthermore, AAt and B have n eigenvectors as a basis for X∈Rn. Further AAt and B commute.Corollary 2.
The integrodifferential equation (1) can be transformed, under Assumption 1, to the following differential equation:
(4)d2u(t)d2t=Adu(t)dt+Bu(t),(5)u(0)=u0,(6)dudt|t=0=Au0,
where the analytical solution is
(7)u(t)=exp⁡(A2t)(cosh⁡(Qt)-A2Q-1sinh⁡(Qt))u0,(8)+exp⁡(A2t)Q-1sinh⁡(Qt)Au0,
where Q=(1/2)AAt/4+B.
WhileAAtand B commute, we also have the commutation of A and Q and, for small t, we have the solution
(9)u(t)=exp⁡(A2t)cosh⁡(Qt)u0,=12(exp⁡((A2+Q)t)+exp⁡((A2-Q)t))u0.Proof.
The analytical solution of the second order differential equation (4) is given as follows. We rewrite (4) into
(10)du(t)dt=v(t),dvdt=Av(t)+Bu(t),u(0)=u0,v(0)=dudt|t=0=Au0.
The exponential matrixA∈R2n×2n can be decoupled via the Cayley Hamilton theorem as
(11)exp⁡(A)=s0(t)I+s1(t)A,
where I∈R2n×2n is the identity matrix.
The characteristic polynomial ofA is given as
(12)z2-Az-B=(z-λ1)(z-λ2)=0,
and the roots are given as
(13)λ1,2=A2±AAt4+B.
Thens0(t)=exp⁡((A/2)t)(cosh⁡(Qt)-(A/2)Q-1sinh⁡(Qt)) and s1(t)=exp⁡((A/2)t)Q-1sinh⁡(Qt); see also [4]. Therefore the analytical solution is given by (7).
We apply this solution to our integrodifferential equation (1) and, with Assumption 1, the equation is satisfied.A further simplification can be made to rewrite the integral-differential equation as two first order differential equations. Later such a reduction will allow us to apply fast iterative splitting methods.Corollary 3.
The integrodifferential equation (1) can be transformed, with the help of Assumption 1, to two first order differential equations:
(14)du1(t)dt=A~u1(t)+B~u1(t),(15)u1(0)=u0,(16)du2(t)dt=A~u2(t)-B~u2(t),(17)u2(0)=u0,
where A~=A/2 and B~=(1/2)AAt/4+B.
The analytical solution is(18)u(t)=12(u1(t)+u2(t))=12(exp⁡((A2+12AAt4+B)t)+exp⁡((A2-12AAt4+B)t))u0.Proof.
The analytical solution of the first order differential equations (14) and (16) is obtained using their characteristic polynomials:
(19)λ1-(A~-B~)=0,λ2-(A~+B~)=0;
hence the solution is, using the notations of A~ and B~,
(20)λ1,2=A~±B~,
and therefore the analytical solution is given by (18) with A~=A/2,B~=(1/2)AAt/4+B.
Therefore this is the solution of our integrodifferential equation (1) under Assumption 1.


## 3. Iterative Splitting Approach to the Integrodifferential Equation

Operator splitting methods are used to solve complex models in geophysical and environmental physics. They have been developed and applied in [5–7]. Such ideas are the basis for this paper: solving the simpler equations but using higher order discretisation methods for the remaining, more complicated equations. With this aim, we use the operator splitting method and decouple the system of equations as will be described in the following.We will concentrate on the iterative splitting method. This algorithm is based on an iteration with fixed splitting discretisation step-sizeτ; namely, on the time interval [tn,tn+1] we solve the following subproblems consecutively for i=0,2,…,2m (cf. [8, 9]):
(21)∂ci(t)∂t=Aci(t)+B∫0t‍ci-1(s)ds,withci(tn)=cn,c0(tn)=cn,c-1=0.0,∂ci+1(t)∂t=Aci(t)+B∫0t‍ci+1(s)ds,withci+1(tn)=cn,
where cn is the known split approximation at the time level t=tn. The split approximation at the time level t=tn+1 is defined by cn+1=c2m+1(tn+1). (Clearly, the function ci+1(t) depends on the interval [tn,tn+1], too, but, for the sake of simplicity, in our notation, we omit the dependence on n.)For the simplification, we apply two standard iterative splitting schemes to the first order differential equations.The first iterative process is(22)∂c1,i(t)∂t=A~c1,i(t)+B~c1,i-1(t),withc1,i(tn)=cn,c1,0(tn)=cn,c1,-1=0.0,∂c1,i+1(t)∂t=A~c1,i(t)+B~c1,i+1(t),withc1,i+1(tn)=cn,
and the second iterative process is
(23)∂c2,i(t)∂t=A~c2,i(t)-B~c2,i-1(t),withc2,i(tn)=cn,c2,0(tn)=cn,c1,-1=0.0,∂c2,i+1(t)∂t=A~c2,i(t)-B~c2,i+1(t),withc2,i+1(tn)=cn,
where cn is the known split approximation at the time level t=tn. The split approximation at the time level t=tn+1 is defined by cn+1=(1/2)(c1,2m+1(tn+1)+c2,2m+1(tn+1)).In the following, we will analyse the convergence and the rate of convergence of the methods (22) and (23) for m tending to infinity for the linear operators A~,B~:X→X using Assumption 1.It suffices to discuss one case, since the other case is the same with only a change in the sign of the operatorB~.Theorem 4.
Let us consider the abstract Cauchy problem in a Banach spaceX:(24)∂tc(t)=A~c(t)+B~c(t),0<t≤T,c(0)=c0,
where A~,B~,A~+B~:X→X are given linear operators which are generators of the C0-semigroup and c0∈X is a given element. Then the iteration process (22) is convergent and the rate of convergence is of second order.Proof.
The proof is given in [10].In the following Algorithm5, we have combined the first and second iterative processes.Algorithm 5.
We apply(-B~i)=(-1)iB~ and exp⁡(B~t)exp⁡(-B~t)=I. The time-steps are tj,j=0,…,N, with tj+1-tj=Δt, and where I∈N+ is a given number of iterative steps.
We have the following algorithm:∂c1/2i(t)/∂t=A~c1/2i(t)±B~ci-1(t),withc1/2i(tn)=c1/2i-1(tn+1),the starting valuesc1/20(tn)=c(tn),results of the last iteration,c1/2-1(tn)=0.0,∂c1/2i+1(t)/∂t=A~c1/2i(t)±B~c1/2i+1(t),withc1/2i+1(tn)=c1/2i(tn+1),ifϵ>max⁡j=1,2⁡|cji+1(tn+1)-cji-1(tn+1)|ori+1=Iwestopc(tn+1)=(1/2)(c1i+1(tn+1)+c2i+1(tn+1))we go to the next time-stepn=n+1tilln+1=N,elsei=i+1.


## 4. Numerical Experiments

For the applications part, we will now treat test problems and real-life applications to transport-scattering problems.The first examples are to validate our novel splitting and solver methods; then, in the subsequent examples, we use more delicate coupled differential equations.


### 4.1. First Example: A Matrix Problem with Integral Term

We deal with simpler integrodifferential equations:(25)c′(t)=c+∫0t‍c(t)dt,t∈[0,1],
where we take ∫0t‍c(t)dt=tc(t) as a first order approximation of the integral and treat
(26)c′(t)=c+tc(t),t∈[0,T],c(0)=c(t0)=1,
where T=10.0 and we have the analytical solution for the approximation
(27)c(t)=exp⁡(t+t22)u(0).We split this into(28)A=1,B(t)=t.We have the following solutions for the iterative scheme:(29)c1(t)=exp⁡(A(tn+1-t))c(tn),t∈(tn,tn+1],c2(t)=exp⁡(∫tntn+1‍B(t)dt)c(tn)+∫tntn+1‍exp⁡(∫stn+1‍B(tn+1-t))dtAc1(s)ds,0000000000000000000000000t∈(tn,tn+1],
where n=0,1,…,N and tN=T while the time-steps are given by Δt=tn+1-tn.We employ the following recurrence relations, distinguishing between even and odd iterations.For the odd iterations,i=2m+1, for m=0,1,2,…,
(30)ci(t)=exp⁡(A(t-tn))c(tn)+∫tnt‍exp⁡(sA)B(s)ci-1(tn+1-s)ds,0000000000000i00000t∈(tn,tn+1].For the even iterations,i=2m, for m=1,2,…,
(31)ci(t)=exp⁡(∫tntn+1‍B(s)ds)c(tn)+∫tnt‍exp⁡(∫0s‍B(t)dt)Aci-1(tn+1-s)ds,00000000000000000i000000t∈(tn,tn+1].
We have a blow-up solution of exponential character in long time behavior; see Figure 1. Therefore, it is advisable that the integrator should be highly accurate for long time behavior.Short time and long time behavior of the solution ((a)T=1, (b) T=10).
(a)
(b)Table1 gives the numerical results of the iterative splitting scheme.Table 1
Numerical experiment with 10 iterative steps for the first example.
Δ
t = 1
Δ
t = 0.5
Δ
t = 0.25
Δ
t = 2−3
Δ
t = 2−4
c
1
1.7634
0.4958
0.1793
0.0753
0.0343
c
2
0.8628
0.1444
0.0282
0.0061
0.0014
c
3
0.2220
0.0104
5.2127
e
-
04
2.8455
e
-
05
1.6511
e
-
06
c
4
0.1116
0.0041
1.8660
e
-
04
9.7846
e
-
06
5.5769
e
-
07
c
5
0.0971
0.0039
1.8367
e
-
04
9.7418
e
-
06
5.5644
e
-
07
c
6
0.0956
0.0039
1.8365
e
-
04
9.7418
e
-
06
5.5644
e
-
07
c
7
0.0954
0.0039
1.8365
e
-
04
9.7416
e
-
06
5.5644
e
-
07
c
8
0.0954
0.0039
1.8366
e
-
04
9.7343
e
-
06
5.5556
e
-
07Figure2 presents the one-sided and two-sided iterative results.Numerical errors of the one-sided and two-sided splitting schemes: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,8 iterative steps.
(a)
(b)Remark 6.
In the experiments, we obtain improved results with each additional step. By the way, the solution blows up and we have to use also very fine time-steps to control the errors. Optimal results are obtained by using the integral part (stiff part) as the implicit part in the iteration (one-sided overB).


### 4.2. Second Example: Integrodifferential Equation

We treat the system of integrodifferential equations(32)c′(t)=Ac+B∫0t‍c(t)dt,t∈[0,1],
where we have the analytical solution
(33)c(t)=exp⁡(At+12Bt2)c(0),t∈[0,T],c(0)=c(t0)=(1,…,1)t,
where A,B∈IRm×m and, for example, we take m=10.We have employed the integral approximation(34)c′(t)=Ac+B∫0t‍c(t)dtc′(t)=Ac+Btc(t).The matrices are given by(35)A=Λ1+Λ2=(-λ10…00-λ10⋮⋮⋱⋱⋱0…0-λ1)+(00…0λ200⋮⋮⋱⋱⋱0…λ20),B=(-λ30…00-λ30⋮⋮⋱⋱⋱0…0-λ3),
where the parameters are λ1=0.02,λ2=0.02, and λ3=0.01.The splitting into the two operatorsA and B is made in order to accelerate the computation of the exponential matrices exp⁡(A) and exp⁡(B).(1) The operatorA is computed as
(36)exp⁡(At)=exp⁡(Λ1t)exp⁡(Λ2t),exp⁡(Λ1t)=(exp⁡(-λ1t)0⋯00exp⁡(-λ1t)0⋮⋮⋱⋱⋮0⋯0exp⁡(-λ1t)),exp⁡⁡(Λ2t)=I+Λ2t+12Λ22t2+16Λ23t3+⋯+1(q-1)!Λ2q-1tq-1,
where Λ2q=0.We have the following solutions for the iterative scheme:(37)c1(t)=exp⁡(A(tn+1-t))c(tn),t∈(tn,tn+1],c2(t)=exp⁡(B(tn+1-t)22)c(tn)+∫tntn+1‍exp⁡⁡(B(t-s)22)Ac1(s)ds,00000000000000i000t∈(tn,tn+1],
where n=0,1,…,N and tN=T while the time-steps are given by Δt=tn+1-tn.We employ the following recurrence relations, distinguishing between even and odd iterations.For the odd iterations,i=2m+1, for m=0,1,2,…,
(38)ci(t)=exp⁡(A(t-tn))c(tn)+∫tnt‍exp⁡⁡((t-s)A)Bsci-1(s)ds,00000000000000i000t∈(tn,tn+1].For the even iterations,i=2m, for m=1,2,…,
(39)ci(t)=exp⁡(B(tn+1-t)22)c(tn)+∫tnt‍exp⁡⁡(B(t-s)22)Aci-1(s)ds,00000000000000i000t∈(tn,tn+1].For the integral computations, we apply higher order schemes, the Newton-Cotes formula (e.g., Simpson’s rule and Boole’s rule); see [11].We acchieve for the second example (vectorial benchmark problem) the same accuracy for its first component as in the first example (scalar benchmark problem).the same accuracy for its first component as in the first example (scalar benchmark problem).Figure3 presents the one-sided and two-sided iterative results.Numerical errors of the one-sided and two-sided splitting schemes: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)The computational results are given in Figure4, which presents the one-sided and two-sided iterative results.The computational time of the one-sided and two-sided splitting schemes: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)Remark 7.
In the experiments, we obtain improved results with each additional step. By the way, the solution blows up and so we have to use very fine time-steps to control the errors. The optimal results are obtained by using the integral part (stiff part) as the implicit part in the iteration (one-sided overB).


#### 4.2.1. Exact Integral Part Based on the Splitting Approach

We treat a system of integrodifferential equations(40)c′(t)=Ac+B∫0t‍c(t)dt,t∈[0,1],
where we assume that the primitive of c is C(0)=0.We have the analytical solution(41)u(t)=12(exp⁡((A2+AAt4+B)t)+exp⁡((A2-AAt4+B)t))u0,
where A~=A/2,B~=AAt/4+B∈Rm×m and, for example, we take m=10.The treatment of the square root for matrices is given in [12, 13].We have the following matrices:(42)A=Λ1+Λ2=(-λ10…00-λ10⋮⋮⋱⋱⋱0…0-λ1)+(00…0λ200⋮⋮⋱⋱⋱0…λ20),B=(-λ30…00-λ30⋮⋮⋱⋱⋱0…0-λ3),
where the parameters are λ1=0.02,λ2=0.02, and λ3=0.01.The exponentials of these operators are computed using Padé approximants; see the discussion in [11].We have the following solutions for the two iterative schemes.Algorithm 8.
Iterative process forc1:
(43)c1,1(t)=exp⁡(A~(tn+1-t))c(tn),t∈(tn,tn+1],c1,2(t)=exp⁡⁡(B~(tn+1-t))c(tn)+∫tntn+1‍exp⁡⁡((tn+1-s)B~)A~c1,1(s)ds,00000000000i00000000t∈(tn,tn+1].
For the odd iterations, i=2m+1, we have, for m=0,1,2,…,
(44)c1,i(t)=exp⁡(A~(t-tn))c(tn)+∫tnt‍exp⁡⁡((t-s)A~)B~c1,i-1(s)ds,00000000000000000t∈(tn,tn+1].
For the even iterations, i=2m, we have, for m=1,2,…,
(45)c1,i(t)=exp⁡(B~(t-tn))c(tn)+∫tnt‍exp⁡⁡(B~(t-s))A~c1,i-1(s)ds,00000000000000000t∈(tn,tn+1],
where n=0,1,…,N, tN=T, and the time-steps are given by Δt=tn+1-tn.Algorithm 9.
Iterative process forc2:
(46)c2,1(t)=exp⁡(A~(tn+1-t))c(tn),t∈(tn,tn+1],c2,2(t)=exp⁡(-B~(tn+1-t))c(tn)+∫tntn+1‍exp⁡⁡(-(tn+1-s)B~)A~c1,1(s)ds,000000000000000000000t∈(tn,tn+1].
For the odd iterations, i=2m+1, we have, for m=0,1,2,…,
(47)c2,i(t)=exp⁡(A~(t-tn))c(tn)+∫tnt‍exp⁡((t-s)A~)(-B~)c2,i-1(s)ds,00000000000000000000t∈(tn,tn+1].
For the even iterations, i=2m, we have, for m=1,2,…,
(48)c2,i(t)=exp⁡(-B~(t-tn))c(tn)+∫tnt‍exp⁡⁡(-B~(t-s))A~c2,i-1(s)ds,0000000000000000I0t∈(tn,tn+1],
where n=0,1,…,N, tN=T, and the time-steps are given by Δt=tn+1-tn.The solution after each iterative step with the two single solutions is(49)ci(t)=12(c1,i(t)+c2,i(t)).For the integral computations, we apply higher order schemes, Newton-Cotes formulas (e.g., Simpson’s rule and Boole’s rule); see [11].Table2 shows the numerical results of the iterative split-ting scheme. We also acchieve similar results for the first component of the vectorial benchmark problem as in its scalar form in the first experiment.Table 2
Numerical experiment with three iterative steps for the second example and the first component.
Δ
t = 1
Δ
t = 0.5
Δ
t = 0.25
Δ
t = 2−3
Δ
t = 2−4
c
1
1.7634
0.4958
0.1793
0.0753
0.0343
c
2
0.8628
0.1444
0.0282
0.0061
0.0014
c
3
0.2220
0.0104
5.2127
e
-
04
2.8455
e
-
05
1.6511
e
-
06Figure5 presents the one-sided and two-sided iterative results for an optimised approach.Numerical errors of the one-sided and two-sided splitting schemes for an optimised approach: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)The computational results are given in Figure6; we present the one-sided and two-sided iterative results for an optimised approach.The computational time of the one-sided and two-sided splitting schemes for an optimised approach: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)Remark 10.
In the experiments, we obtain improved results with each additional step. By the way, the solution blows up and we have to use very fine time-steps to control the errors. The optimal results are obtained by using the integral part (stiff part) as the implicit part in the iteration (one-sided overB).Remark 11.
To optimise the algorithms, we combine both of the iterative cycles with the idea of-Bi=(-1)iB and in the Padé approximants we apply exp⁡(Bt)exp⁡(-Bt)=I; see Algorithm 5.Remark 12.
Table3 presents the computational times of the three experiments. While the exact approach based on the splitting idea of higher order differential equations obtains more accurate solutions, its drawback is its higher computational costs. We reduce such costs by using the symmetry of the underlying iterative schemes and reach nearly the same computational time as in the numerical approach scheme.Table 3
The computational time of all three experiments.
Algorithm
Computational time in sec
num. approach—twoSide
24.3976
num. approach—oneSideA
24.7610
num. approach—oneSideB
24.7252
exact approach—twoSide
67.0159
exact approach—oneSideA
68.4638
exact approach—oneSideB
69.1291
exact optimised approach—twoSide
36.3361
exact optimised approach—oneSideA
36.9736
exact optimised approach—oneSideB
35.8903


### 4.3. Real-Life Example: One-Phase Problem

The next example is a simplified real-life problem with a neutron transport equation, which includes the gain and loss of a neutron.We concentrate on the computational benefits of a fast computation of the iterative scheme, given with matrix exponential functions.The equation is(50)∂tc+∇·Fc=-λ1c+∫0t‍λ2c(x,t)dt,inΩ×[0,t],F=v-D∇,c(x,t)=c0(x),on⁡Ω,c(x,t)=c1(x,t),on⁡∂Ω×[0,t].In the following, we deal with the semidiscretised equation given by the matrices(51)∂tC=(A-Λ1)C+∫0t‍Λ2C(s)ds,
where C=(c1,…,cm)T is the solution of the species in the mobile phase in each spatial discretisation point (i=1,…,m).We have the following two operators for the splitting method:(52)A=DΔx2·(-211-21⋱⋱⋱1-211-2)+vΔx·(1-11⋱⋱-11-11)∈Rm×m.
The reaction and scattering matrices are
(53)Λ1=(λ100λ10⋱⋱⋱0λ100λ1)∈Rm×m,Λ2=(λ200λ20⋱⋱⋱0λ200λ2)∈Rm×m,
where m is the number of spatial points.In the first algorithm, we use the following approximation of the integral term:(54)∫0tλ2c(x,t)dt≈λ2tc(x,t).We use the following recurrence relations, distinguishing between even and odd iterations.For the odd iterations,i=2m~+1, with m~=0,1,2,…, we have
(55)Ci(t)=exp⁡(A(t-tn))C(tn)+∫tnt‍exp⁡⁡((t-s)A)B(s)Ci-1(s)ds,00000000000i0000000t∈(tn,tn+1].For the even iterations,i=2m~, with m~=0,1,2,…, we have
(56)Ci(t)=exp⁡(∫0t-tn‍B(s~)ds~)C(tn)+∫tnt‍exp⁡⁡(∫0t-s‍B(s~)ds~)ACi-1(s)ds,00000000000i000000000t∈(tn,tn+1],
where ∫0t-tn‍B(s~)ds~=Λ1(t-tn)+Λ2((t-tn)2/2).For the reference solution, we apply a fine time-spatial scale without decoupling the equations.For the algorithms, we employ novel Algorithm5 and we obtain A~=A-Λ1 and B~=Λ2.The matrices for the novel algorithm are(57)A~~=A~2,B~~=12A~A~t4+B~.Figure7 presents the one-sided and two-sided iterative results.Numerical errors of the one-sided and two-sided splitting schemes: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)The computational results are given in Figure8, which presents the one-sided and two-sided iterative results.The computational time of the one-sided and two-sided splitting schemes: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)Figure9 presents the one-sided and two-sided iterative results for an optimised approach.Numerical errors of the one-sided and two-sided splitting schemes for an optimised approach: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)The computational results are given in Figure10, which presents the one-sided and two-sided iterative results for an optimised approach.The computational time of the one-sided and two-sided splitting schemes for an optimised approach: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)Table4 presents the computational times of the three experiments.Table 4
The computational time of all three experiments.
Algorithm
Computational time in sec
numerical approach—twoSide
23.8826
numerical approach—oneSideA
25.4946
numerical approach—oneSideB
24.1386
exact approach—twoSide
67.1759
exact approach—oneSideA
62.7140
exact approach—oneSideB
70.8573
exact optimised approach—twoSide
36.4432
exact optimised approach—oneSideA
29.7267
exact optimised approach—oneSideB
38.1619Remark 13.
For all iterative schemes, we can reach faster results using the optimised approach. Such iterative schemes with fast computations of the exponential matrices accelerate standard splitting schemes. With four or five iterative steps we obtain more accurate results than we did using the expensive standard schemes. One-sided iterative schemes yield the best convergence results.


## 4.1. First Example: A Matrix Problem with Integral Term

We deal with simpler integrodifferential equations:(25)c′(t)=c+∫0t‍c(t)dt,t∈[0,1],
where we take ∫0t‍c(t)dt=tc(t) as a first order approximation of the integral and treat
(26)c′(t)=c+tc(t),t∈[0,T],c(0)=c(t0)=1,
where T=10.0 and we have the analytical solution for the approximation
(27)c(t)=exp⁡(t+t22)u(0).We split this into(28)A=1,B(t)=t.We have the following solutions for the iterative scheme:(29)c1(t)=exp⁡(A(tn+1-t))c(tn),t∈(tn,tn+1],c2(t)=exp⁡(∫tntn+1‍B(t)dt)c(tn)+∫tntn+1‍exp⁡(∫stn+1‍B(tn+1-t))dtAc1(s)ds,0000000000000000000000000t∈(tn,tn+1],
where n=0,1,…,N and tN=T while the time-steps are given by Δt=tn+1-tn.We employ the following recurrence relations, distinguishing between even and odd iterations.For the odd iterations,i=2m+1, for m=0,1,2,…,
(30)ci(t)=exp⁡(A(t-tn))c(tn)+∫tnt‍exp⁡(sA)B(s)ci-1(tn+1-s)ds,0000000000000i00000t∈(tn,tn+1].For the even iterations,i=2m, for m=1,2,…,
(31)ci(t)=exp⁡(∫tntn+1‍B(s)ds)c(tn)+∫tnt‍exp⁡(∫0s‍B(t)dt)Aci-1(tn+1-s)ds,00000000000000000i000000t∈(tn,tn+1].
We have a blow-up solution of exponential character in long time behavior; see Figure 1. Therefore, it is advisable that the integrator should be highly accurate for long time behavior.Short time and long time behavior of the solution ((a)T=1, (b) T=10).
(a)
(b)Table1 gives the numerical results of the iterative splitting scheme.Table 1
Numerical experiment with 10 iterative steps for the first example.
Δ
t = 1
Δ
t = 0.5
Δ
t = 0.25
Δ
t = 2−3
Δ
t = 2−4
c
1
1.7634
0.4958
0.1793
0.0753
0.0343
c
2
0.8628
0.1444
0.0282
0.0061
0.0014
c
3
0.2220
0.0104
5.2127
e
-
04
2.8455
e
-
05
1.6511
e
-
06
c
4
0.1116
0.0041
1.8660
e
-
04
9.7846
e
-
06
5.5769
e
-
07
c
5
0.0971
0.0039
1.8367
e
-
04
9.7418
e
-
06
5.5644
e
-
07
c
6
0.0956
0.0039
1.8365
e
-
04
9.7418
e
-
06
5.5644
e
-
07
c
7
0.0954
0.0039
1.8365
e
-
04
9.7416
e
-
06
5.5644
e
-
07
c
8
0.0954
0.0039
1.8366
e
-
04
9.7343
e
-
06
5.5556
e
-
07Figure2 presents the one-sided and two-sided iterative results.Numerical errors of the one-sided and two-sided splitting schemes: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,8 iterative steps.
(a)
(b)Remark 6.
In the experiments, we obtain improved results with each additional step. By the way, the solution blows up and we have to use also very fine time-steps to control the errors. Optimal results are obtained by using the integral part (stiff part) as the implicit part in the iteration (one-sided overB).


## 4.2. Second Example: Integrodifferential Equation

We treat the system of integrodifferential equations(32)c′(t)=Ac+B∫0t‍c(t)dt,t∈[0,1],
where we have the analytical solution
(33)c(t)=exp⁡(At+12Bt2)c(0),t∈[0,T],c(0)=c(t0)=(1,…,1)t,
where A,B∈IRm×m and, for example, we take m=10.We have employed the integral approximation(34)c′(t)=Ac+B∫0t‍c(t)dtc′(t)=Ac+Btc(t).The matrices are given by(35)A=Λ1+Λ2=(-λ10…00-λ10⋮⋮⋱⋱⋱0…0-λ1)+(00…0λ200⋮⋮⋱⋱⋱0…λ20),B=(-λ30…00-λ30⋮⋮⋱⋱⋱0…0-λ3),
where the parameters are λ1=0.02,λ2=0.02, and λ3=0.01.The splitting into the two operatorsA and B is made in order to accelerate the computation of the exponential matrices exp⁡(A) and exp⁡(B).(1) The operatorA is computed as
(36)exp⁡(At)=exp⁡(Λ1t)exp⁡(Λ2t),exp⁡(Λ1t)=(exp⁡(-λ1t)0⋯00exp⁡(-λ1t)0⋮⋮⋱⋱⋮0⋯0exp⁡(-λ1t)),exp⁡⁡(Λ2t)=I+Λ2t+12Λ22t2+16Λ23t3+⋯+1(q-1)!Λ2q-1tq-1,
where Λ2q=0.We have the following solutions for the iterative scheme:(37)c1(t)=exp⁡(A(tn+1-t))c(tn),t∈(tn,tn+1],c2(t)=exp⁡(B(tn+1-t)22)c(tn)+∫tntn+1‍exp⁡⁡(B(t-s)22)Ac1(s)ds,00000000000000i000t∈(tn,tn+1],
where n=0,1,…,N and tN=T while the time-steps are given by Δt=tn+1-tn.We employ the following recurrence relations, distinguishing between even and odd iterations.For the odd iterations,i=2m+1, for m=0,1,2,…,
(38)ci(t)=exp⁡(A(t-tn))c(tn)+∫tnt‍exp⁡⁡((t-s)A)Bsci-1(s)ds,00000000000000i000t∈(tn,tn+1].For the even iterations,i=2m, for m=1,2,…,
(39)ci(t)=exp⁡(B(tn+1-t)22)c(tn)+∫tnt‍exp⁡⁡(B(t-s)22)Aci-1(s)ds,00000000000000i000t∈(tn,tn+1].For the integral computations, we apply higher order schemes, the Newton-Cotes formula (e.g., Simpson’s rule and Boole’s rule); see [11].We acchieve for the second example (vectorial benchmark problem) the same accuracy for its first component as in the first example (scalar benchmark problem).the same accuracy for its first component as in the first example (scalar benchmark problem).Figure3 presents the one-sided and two-sided iterative results.Numerical errors of the one-sided and two-sided splitting schemes: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)The computational results are given in Figure4, which presents the one-sided and two-sided iterative results.The computational time of the one-sided and two-sided splitting schemes: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)Remark 7.
In the experiments, we obtain improved results with each additional step. By the way, the solution blows up and so we have to use very fine time-steps to control the errors. The optimal results are obtained by using the integral part (stiff part) as the implicit part in the iteration (one-sided overB).


### 4.2.1. Exact Integral Part Based on the Splitting Approach

We treat a system of integrodifferential equations(40)c′(t)=Ac+B∫0t‍c(t)dt,t∈[0,1],
where we assume that the primitive of c is C(0)=0.We have the analytical solution(41)u(t)=12(exp⁡((A2+AAt4+B)t)+exp⁡((A2-AAt4+B)t))u0,
where A~=A/2,B~=AAt/4+B∈Rm×m and, for example, we take m=10.The treatment of the square root for matrices is given in [12, 13].We have the following matrices:(42)A=Λ1+Λ2=(-λ10…00-λ10⋮⋮⋱⋱⋱0…0-λ1)+(00…0λ200⋮⋮⋱⋱⋱0…λ20),B=(-λ30…00-λ30⋮⋮⋱⋱⋱0…0-λ3),
where the parameters are λ1=0.02,λ2=0.02, and λ3=0.01.The exponentials of these operators are computed using Padé approximants; see the discussion in [11].We have the following solutions for the two iterative schemes.Algorithm 8.
Iterative process forc1:
(43)c1,1(t)=exp⁡(A~(tn+1-t))c(tn),t∈(tn,tn+1],c1,2(t)=exp⁡⁡(B~(tn+1-t))c(tn)+∫tntn+1‍exp⁡⁡((tn+1-s)B~)A~c1,1(s)ds,00000000000i00000000t∈(tn,tn+1].
For the odd iterations, i=2m+1, we have, for m=0,1,2,…,
(44)c1,i(t)=exp⁡(A~(t-tn))c(tn)+∫tnt‍exp⁡⁡((t-s)A~)B~c1,i-1(s)ds,00000000000000000t∈(tn,tn+1].
For the even iterations, i=2m, we have, for m=1,2,…,
(45)c1,i(t)=exp⁡(B~(t-tn))c(tn)+∫tnt‍exp⁡⁡(B~(t-s))A~c1,i-1(s)ds,00000000000000000t∈(tn,tn+1],
where n=0,1,…,N, tN=T, and the time-steps are given by Δt=tn+1-tn.Algorithm 9.
Iterative process forc2:
(46)c2,1(t)=exp⁡(A~(tn+1-t))c(tn),t∈(tn,tn+1],c2,2(t)=exp⁡(-B~(tn+1-t))c(tn)+∫tntn+1‍exp⁡⁡(-(tn+1-s)B~)A~c1,1(s)ds,000000000000000000000t∈(tn,tn+1].
For the odd iterations, i=2m+1, we have, for m=0,1,2,…,
(47)c2,i(t)=exp⁡(A~(t-tn))c(tn)+∫tnt‍exp⁡((t-s)A~)(-B~)c2,i-1(s)ds,00000000000000000000t∈(tn,tn+1].
For the even iterations, i=2m, we have, for m=1,2,…,
(48)c2,i(t)=exp⁡(-B~(t-tn))c(tn)+∫tnt‍exp⁡⁡(-B~(t-s))A~c2,i-1(s)ds,0000000000000000I0t∈(tn,tn+1],
where n=0,1,…,N, tN=T, and the time-steps are given by Δt=tn+1-tn.The solution after each iterative step with the two single solutions is(49)ci(t)=12(c1,i(t)+c2,i(t)).For the integral computations, we apply higher order schemes, Newton-Cotes formulas (e.g., Simpson’s rule and Boole’s rule); see [11].Table2 shows the numerical results of the iterative split-ting scheme. We also acchieve similar results for the first component of the vectorial benchmark problem as in its scalar form in the first experiment.Table 2
Numerical experiment with three iterative steps for the second example and the first component.
Δ
t = 1
Δ
t = 0.5
Δ
t = 0.25
Δ
t = 2−3
Δ
t = 2−4
c
1
1.7634
0.4958
0.1793
0.0753
0.0343
c
2
0.8628
0.1444
0.0282
0.0061
0.0014
c
3
0.2220
0.0104
5.2127
e
-
04
2.8455
e
-
05
1.6511
e
-
06Figure5 presents the one-sided and two-sided iterative results for an optimised approach.Numerical errors of the one-sided and two-sided splitting schemes for an optimised approach: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)The computational results are given in Figure6; we present the one-sided and two-sided iterative results for an optimised approach.The computational time of the one-sided and two-sided splitting schemes for an optimised approach: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)Remark 10.
In the experiments, we obtain improved results with each additional step. By the way, the solution blows up and we have to use very fine time-steps to control the errors. The optimal results are obtained by using the integral part (stiff part) as the implicit part in the iteration (one-sided overB).Remark 11.
To optimise the algorithms, we combine both of the iterative cycles with the idea of-Bi=(-1)iB and in the Padé approximants we apply exp⁡(Bt)exp⁡(-Bt)=I; see Algorithm 5.Remark 12.
Table3 presents the computational times of the three experiments. While the exact approach based on the splitting idea of higher order differential equations obtains more accurate solutions, its drawback is its higher computational costs. We reduce such costs by using the symmetry of the underlying iterative schemes and reach nearly the same computational time as in the numerical approach scheme.Table 3
The computational time of all three experiments.
Algorithm
Computational time in sec
num. approach—twoSide
24.3976
num. approach—oneSideA
24.7610
num. approach—oneSideB
24.7252
exact approach—twoSide
67.0159
exact approach—oneSideA
68.4638
exact approach—oneSideB
69.1291
exact optimised approach—twoSide
36.3361
exact optimised approach—oneSideA
36.9736
exact optimised approach—oneSideB
35.8903


## 4.2.1. Exact Integral Part Based on the Splitting Approach

We treat a system of integrodifferential equations(40)c′(t)=Ac+B∫0t‍c(t)dt,t∈[0,1],
where we assume that the primitive of c is C(0)=0.We have the analytical solution(41)u(t)=12(exp⁡((A2+AAt4+B)t)+exp⁡((A2-AAt4+B)t))u0,
where A~=A/2,B~=AAt/4+B∈Rm×m and, for example, we take m=10.The treatment of the square root for matrices is given in [12, 13].We have the following matrices:(42)A=Λ1+Λ2=(-λ10…00-λ10⋮⋮⋱⋱⋱0…0-λ1)+(00…0λ200⋮⋮⋱⋱⋱0…λ20),B=(-λ30…00-λ30⋮⋮⋱⋱⋱0…0-λ3),
where the parameters are λ1=0.02,λ2=0.02, and λ3=0.01.The exponentials of these operators are computed using Padé approximants; see the discussion in [11].We have the following solutions for the two iterative schemes.Algorithm 8.
Iterative process forc1:
(43)c1,1(t)=exp⁡(A~(tn+1-t))c(tn),t∈(tn,tn+1],c1,2(t)=exp⁡⁡(B~(tn+1-t))c(tn)+∫tntn+1‍exp⁡⁡((tn+1-s)B~)A~c1,1(s)ds,00000000000i00000000t∈(tn,tn+1].
For the odd iterations, i=2m+1, we have, for m=0,1,2,…,
(44)c1,i(t)=exp⁡(A~(t-tn))c(tn)+∫tnt‍exp⁡⁡((t-s)A~)B~c1,i-1(s)ds,00000000000000000t∈(tn,tn+1].
For the even iterations, i=2m, we have, for m=1,2,…,
(45)c1,i(t)=exp⁡(B~(t-tn))c(tn)+∫tnt‍exp⁡⁡(B~(t-s))A~c1,i-1(s)ds,00000000000000000t∈(tn,tn+1],
where n=0,1,…,N, tN=T, and the time-steps are given by Δt=tn+1-tn.Algorithm 9.
Iterative process forc2:
(46)c2,1(t)=exp⁡(A~(tn+1-t))c(tn),t∈(tn,tn+1],c2,2(t)=exp⁡(-B~(tn+1-t))c(tn)+∫tntn+1‍exp⁡⁡(-(tn+1-s)B~)A~c1,1(s)ds,000000000000000000000t∈(tn,tn+1].
For the odd iterations, i=2m+1, we have, for m=0,1,2,…,
(47)c2,i(t)=exp⁡(A~(t-tn))c(tn)+∫tnt‍exp⁡((t-s)A~)(-B~)c2,i-1(s)ds,00000000000000000000t∈(tn,tn+1].
For the even iterations, i=2m, we have, for m=1,2,…,
(48)c2,i(t)=exp⁡(-B~(t-tn))c(tn)+∫tnt‍exp⁡⁡(-B~(t-s))A~c2,i-1(s)ds,0000000000000000I0t∈(tn,tn+1],
where n=0,1,…,N, tN=T, and the time-steps are given by Δt=tn+1-tn.The solution after each iterative step with the two single solutions is(49)ci(t)=12(c1,i(t)+c2,i(t)).For the integral computations, we apply higher order schemes, Newton-Cotes formulas (e.g., Simpson’s rule and Boole’s rule); see [11].Table2 shows the numerical results of the iterative split-ting scheme. We also acchieve similar results for the first component of the vectorial benchmark problem as in its scalar form in the first experiment.Table 2
Numerical experiment with three iterative steps for the second example and the first component.
Δ
t = 1
Δ
t = 0.5
Δ
t = 0.25
Δ
t = 2−3
Δ
t = 2−4
c
1
1.7634
0.4958
0.1793
0.0753
0.0343
c
2
0.8628
0.1444
0.0282
0.0061
0.0014
c
3
0.2220
0.0104
5.2127
e
-
04
2.8455
e
-
05
1.6511
e
-
06Figure5 presents the one-sided and two-sided iterative results for an optimised approach.Numerical errors of the one-sided and two-sided splitting schemes for an optimised approach: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)The computational results are given in Figure6; we present the one-sided and two-sided iterative results for an optimised approach.The computational time of the one-sided and two-sided splitting schemes for an optimised approach: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)Remark 10.
In the experiments, we obtain improved results with each additional step. By the way, the solution blows up and we have to use very fine time-steps to control the errors. The optimal results are obtained by using the integral part (stiff part) as the implicit part in the iteration (one-sided overB).Remark 11.
To optimise the algorithms, we combine both of the iterative cycles with the idea of-Bi=(-1)iB and in the Padé approximants we apply exp⁡(Bt)exp⁡(-Bt)=I; see Algorithm 5.Remark 12.
Table3 presents the computational times of the three experiments. While the exact approach based on the splitting idea of higher order differential equations obtains more accurate solutions, its drawback is its higher computational costs. We reduce such costs by using the symmetry of the underlying iterative schemes and reach nearly the same computational time as in the numerical approach scheme.Table 3
The computational time of all three experiments.
Algorithm
Computational time in sec
num. approach—twoSide
24.3976
num. approach—oneSideA
24.7610
num. approach—oneSideB
24.7252
exact approach—twoSide
67.0159
exact approach—oneSideA
68.4638
exact approach—oneSideB
69.1291
exact optimised approach—twoSide
36.3361
exact optimised approach—oneSideA
36.9736
exact optimised approach—oneSideB
35.8903


## 4.3. Real-Life Example: One-Phase Problem

The next example is a simplified real-life problem with a neutron transport equation, which includes the gain and loss of a neutron.We concentrate on the computational benefits of a fast computation of the iterative scheme, given with matrix exponential functions.The equation is(50)∂tc+∇·Fc=-λ1c+∫0t‍λ2c(x,t)dt,inΩ×[0,t],F=v-D∇,c(x,t)=c0(x),on⁡Ω,c(x,t)=c1(x,t),on⁡∂Ω×[0,t].In the following, we deal with the semidiscretised equation given by the matrices(51)∂tC=(A-Λ1)C+∫0t‍Λ2C(s)ds,
where C=(c1,…,cm)T is the solution of the species in the mobile phase in each spatial discretisation point (i=1,…,m).We have the following two operators for the splitting method:(52)A=DΔx2·(-211-21⋱⋱⋱1-211-2)+vΔx·(1-11⋱⋱-11-11)∈Rm×m.
The reaction and scattering matrices are
(53)Λ1=(λ100λ10⋱⋱⋱0λ100λ1)∈Rm×m,Λ2=(λ200λ20⋱⋱⋱0λ200λ2)∈Rm×m,
where m is the number of spatial points.In the first algorithm, we use the following approximation of the integral term:(54)∫0tλ2c(x,t)dt≈λ2tc(x,t).We use the following recurrence relations, distinguishing between even and odd iterations.For the odd iterations,i=2m~+1, with m~=0,1,2,…, we have
(55)Ci(t)=exp⁡(A(t-tn))C(tn)+∫tnt‍exp⁡⁡((t-s)A)B(s)Ci-1(s)ds,00000000000i0000000t∈(tn,tn+1].For the even iterations,i=2m~, with m~=0,1,2,…, we have
(56)Ci(t)=exp⁡(∫0t-tn‍B(s~)ds~)C(tn)+∫tnt‍exp⁡⁡(∫0t-s‍B(s~)ds~)ACi-1(s)ds,00000000000i000000000t∈(tn,tn+1],
where ∫0t-tn‍B(s~)ds~=Λ1(t-tn)+Λ2((t-tn)2/2).For the reference solution, we apply a fine time-spatial scale without decoupling the equations.For the algorithms, we employ novel Algorithm5 and we obtain A~=A-Λ1 and B~=Λ2.The matrices for the novel algorithm are(57)A~~=A~2,B~~=12A~A~t4+B~.Figure7 presents the one-sided and two-sided iterative results.Numerical errors of the one-sided and two-sided splitting schemes: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)The computational results are given in Figure8, which presents the one-sided and two-sided iterative results.The computational time of the one-sided and two-sided splitting schemes: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)Figure9 presents the one-sided and two-sided iterative results for an optimised approach.Numerical errors of the one-sided and two-sided splitting schemes for an optimised approach: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)The computational results are given in Figure10, which presents the one-sided and two-sided iterative results for an optimised approach.The computational time of the one-sided and two-sided splitting schemes for an optimised approach: one-sided splitting overA (a) and two-sided splitting scheme alternating between A and B (b) with 1,…,6 iterative steps.
(a)
(b)Table4 presents the computational times of the three experiments.Table 4
The computational time of all three experiments.
Algorithm
Computational time in sec
numerical approach—twoSide
23.8826
numerical approach—oneSideA
25.4946
numerical approach—oneSideB
24.1386
exact approach—twoSide
67.1759
exact approach—oneSideA
62.7140
exact approach—oneSideB
70.8573
exact optimised approach—twoSide
36.4432
exact optimised approach—oneSideA
29.7267
exact optimised approach—oneSideB
38.1619Remark 13.
For all iterative schemes, we can reach faster results using the optimised approach. Such iterative schemes with fast computations of the exponential matrices accelerate standard splitting schemes. With four or five iterative steps we obtain more accurate results than we did using the expensive standard schemes. One-sided iterative schemes yield the best convergence results.


## 5. Conclusions and Discussion

We presented a coupled model for a transport and kinetic model for deposition species in a plasma environment. We proposed numerical schemes to decompose the flow field and the reaction-collision term. A generalised iterative scheme was discussed and their analysis and benefits were presented.Such combined splitting schemes help to reduce the computational time required for delicate plasma problems. In the future, we will deal with nonlinear schemes of collision problems in plasma reactors.



---

*Source: 812137-2014-08-24.xml*