# Nonlinear Extension of Multiproduct Expansion Schemes and Applications to Rigid Bodies

Authors: Jürgen Geiser
Journal: International Journal of Differential Equations
 (2013)

---

## Abstract

In this paper we discuss time integrators for nonlinear differential equations. In recent years, splitting approaches have become an important tool for reducing the computational time needed to solve differential equations. Moreover, nonlinearity is a challenge to splitting schemes, while one has to extend the exp-functions in terms of a nonlinear Magnus expansion. Here we discuss a novel extension of the so-called multiproduct expansion methods, which is used to improve the standard Strang splitting schemes as to their nonlinearity. We present an extension of linear splitting schemes and concentrate on nonlinear systems of differential equations and generalise in this respect the recent MPE method; see (Chin and Geiser, 2011). Some first numerical examples, of rigid body problems, are given as benchmarks.


---

## Body

## 1. Introduction

In recent years, applications to nonlinear differential equations of multiscale problems, for example, the rigid body (see [1]), have arisen and become important.Here, splitting schemes are important for decoupling the different nonlinear scales, for example, a Hamiltonian with kinetic and potential operators, and treating them with the best solver schemes; see [2, 3].Theoretically, we deal with subproblems which can be solved independently, but we have taken into account their nonlinearity, so that standard fundamental solutions, for example,exp-functions, of the subproblems have to be extended in respect to their nonlinearity.Here we apply the nonlinear Magnus expansion; see [4, 5].Our contributions are to use the first and second order schemes and extend them with an extrapolation scheme, in our case with a multiproduct expansion, to gain more accurate and efficient higher order schemes; see [6].In this paper we concentrate on approximate solutions of the nonlinear evolution equation; for example,(1)∂tu=F(u(t))=A(u(t))+B(u(t)),u(0)=u0,
with unbounded operators A:D(A)⊂X→X and B:D(B)⊂X→X. We have further F(v)=A(v)+B(v), v∈D(A)∩D(B).We assume there are suitable chosen subspaces of the underlying Banach space(X,∥·∥X) such that D(F)=D(A)∩D(B)≠∅.For such equations, we concentrate on applying nonlinear splitting schemes to extrapolation ideas to obtain higher order schemes. Here, we deal with Suzuki's methods and apply factorised symplectic algorithms with forward derivatives; see [6, 7].The exact solution of the evolution problem (1) is
(2)u(t)=ℰF(t,u(0)),0≤t≤T,
with the evolution operator ℰF depending on the actual time t and the initial value u(0). We use a formal notation:
(3)u(t)=exp(tDF)u(0),0≤t≤T.
Here the evolution operator exp(tDF) and the Lie derivative DF associated with F are
(4)exp(tDF)Gv=G(ℰF(t,v)),0≤t≤T,DFGv=G′(V)F(v)
for any unbounded nonlinear operator G:D(G)⊂X→X with Frechet’s derivative G′.The paper is outlined as follows. In Section2, we discuss the standard splitting schemes. The extended nonlinear multiproduct expansion is presented in Section 3. In Section 4, we present the improved methods due to the initialisation process. The alternative methods with respect to the linearisation schemes are presented in Section 5. The numerical experiments are given in Section 6. The conclusion is in Section 7.


## 2. Nonlinear Operator Splitting Methods

In the literature, there are various types of splitting methods. We mainly consider the following operators splitting schemes in this study.(1)
Sequential operator splitting:A-B splitting
(5)∂c*(t)∂t=A(c*(t))witht∈[tn,tn+1],c*(tn)=cspn,(6)∂c**(t)∂t=B(c**(t))witht∈[tn,tn+1],c**(tn)=c*(tn+1),for n=0,1,…,N-1, whereby cspn=c0 is given from (1). The approximated split solution at the point t=tn+1 is defined by cspn+1=c**(tn+1).(2)
The Strang-Marchuk operator splitting:A-B-A splitting
(7)∂c*(t)∂t=A(c*(t))witht∈[tn,tn+1/2],c*(tn)=cspn,∂c**(t)∂t=B(c**(t))witht∈[tn,tn+1/2],c**(tn)=c*(tn+1/2),∂c***(t)∂t=A(c*(t))witht∈[tn+1/2,tn+1],c***(tn+1/2)=c**(tn+1),where tn+1/2=tn+τ/2, τ is the local time-step. The approximate split solution at the point t=tn+1 is defined by cspn+1=c***(tn+1).(3)
Iterative splitting with respect to one operator(8)∂ci(t)∂t=A(ci(t))+B(ci-1(t)),withci(tn)=cn,i=1,2,…,m,where c0(tn)=cn,c-1=0, and 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=cm+1(tn+1).


### 2.1. Error Analysis of the Classical Splitting Schemes

In the following, we discuss the abstract error analysis with respect to the Lie operator conventions.We associate a new operatorF^ which is a linear Lie operator.Definition 1.
For the given operatorF we associate a new operator, denoted by F^.
This Lie operator acts on the space of the differentiable operators of the type(9)S⟶S
and maps each operator G into the new operator F^(G), such that for any element c∈S(10)(F^(G))(c)=(G′(c)∘F)(c).Here the derivatives are given by the following.Definition 2.
Thekth power of the Lie operator F^ applied to some operator G can be expressed as the kth derivative of G; that is, the relation
(11)(F^k(G))(c)=∂kG(c)∂tk
is valid for all k=1,2,….We deal with the problem (1) and split the operator F into the sum A+B; then we have the following results of the abstract splitting errors. We have an appropriate norm ∥·∥ in a general Banach space X.Theorem 3.
The splitting error for theA-B splitting is given by
(12)∥ErrA-B,1,2∥=τ((B′(c)∘A)(c)-(A′(c)∘B)(c))+O(τ2),
The splitting error of the Strang splitting is given by(13)ErrStrang,1,2=124τ2([B,[B,A]](c)-2[A,[A,B]](c))+O(τ4).
with the commutator [A,B](c)=(B′(c)∘A)(c)-(A′(c)∘B)(c).
The splitting error of the iterative splitting is given by(14)∥Erriter,i(t)∥=τn∥F^2(I)∥∥Erriter,i-1(t)∥,forτn=t-tn,
where one assumes that the previous error ∥Erriter(tn)∥=0.Proof.
The splitting error for theA-B splitting is given by
(15)ErrA-B=exp(τA^+B^)(I)-(exp(τA^)exp(τB^))(I)(c).
We evaluate the exp-function and apply the commutator [A,B](c) and obtain
(16)ErrA-B,1,2=τ((B′(c)∘A)(c)-(A′(c)∘B)(c))+O(τ2).
The splitting error of the Strang Splitting is given by(17)ErrStrang=exp(τA^+B^)(I)-(exp(τ2A^)exp(τB^)exp(τ2A^))(I)(c).
We evaluate the exp-function and apply twice the commutator [A,B](c) based on the symmetry of the scheme and obtain
(18)ErrStrang,1,2=124τ2([B,[B,A]](c)-2[A,[A,B]](c))+O(τ4),
where the commutator is given as [A,B](c)=(B′(c)∘A)(c)-(A′(c)∘B)(c).
The splitting error of the iterative splitting is given by(19)∥Erriter,i(t)∥=∥exp((τA^)(I))Err(tn)+∫tnt‍exp(((t-s)A^)(I))B^(I)Erriter,i-1(s)ds∥,≤τn∥B^(I)∥∥Erriter,i-1(t)∥,forτn=t-tn,
where we have applied the method of variation of constants to the linearised operators. Furthermore, we assume that the previous error ∥Erriter(tn)∥=0 and Erriter,i(t)=u(t)-ui(t) is the error of the iterative scheme.


## 2.1. Error Analysis of the Classical Splitting Schemes

In the following, we discuss the abstract error analysis with respect to the Lie operator conventions.We associate a new operatorF^ which is a linear Lie operator.Definition 1.
For the given operatorF we associate a new operator, denoted by F^.
This Lie operator acts on the space of the differentiable operators of the type(9)S⟶S
and maps each operator G into the new operator F^(G), such that for any element c∈S(10)(F^(G))(c)=(G′(c)∘F)(c).Here the derivatives are given by the following.Definition 2.
Thekth power of the Lie operator F^ applied to some operator G can be expressed as the kth derivative of G; that is, the relation
(11)(F^k(G))(c)=∂kG(c)∂tk
is valid for all k=1,2,….We deal with the problem (1) and split the operator F into the sum A+B; then we have the following results of the abstract splitting errors. We have an appropriate norm ∥·∥ in a general Banach space X.Theorem 3.
The splitting error for theA-B splitting is given by
(12)∥ErrA-B,1,2∥=τ((B′(c)∘A)(c)-(A′(c)∘B)(c))+O(τ2),
The splitting error of the Strang splitting is given by(13)ErrStrang,1,2=124τ2([B,[B,A]](c)-2[A,[A,B]](c))+O(τ4).
with the commutator [A,B](c)=(B′(c)∘A)(c)-(A′(c)∘B)(c).
The splitting error of the iterative splitting is given by(14)∥Erriter,i(t)∥=τn∥F^2(I)∥∥Erriter,i-1(t)∥,forτn=t-tn,
where one assumes that the previous error ∥Erriter(tn)∥=0.Proof.
The splitting error for theA-B splitting is given by
(15)ErrA-B=exp(τA^+B^)(I)-(exp(τA^)exp(τB^))(I)(c).
We evaluate the exp-function and apply the commutator [A,B](c) and obtain
(16)ErrA-B,1,2=τ((B′(c)∘A)(c)-(A′(c)∘B)(c))+O(τ2).
The splitting error of the Strang Splitting is given by(17)ErrStrang=exp(τA^+B^)(I)-(exp(τ2A^)exp(τB^)exp(τ2A^))(I)(c).
We evaluate the exp-function and apply twice the commutator [A,B](c) based on the symmetry of the scheme and obtain
(18)ErrStrang,1,2=124τ2([B,[B,A]](c)-2[A,[A,B]](c))+O(τ4),
where the commutator is given as [A,B](c)=(B′(c)∘A)(c)-(A′(c)∘B)(c).
The splitting error of the iterative splitting is given by(19)∥Erriter,i(t)∥=∥exp((τA^)(I))Err(tn)+∫tnt‍exp(((t-s)A^)(I))B^(I)Erriter,i-1(s)ds∥,≤τn∥B^(I)∥∥Erriter,i-1(t)∥,forτn=t-tn,
where we have applied the method of variation of constants to the linearised operators. Furthermore, we assume that the previous error ∥Erriter(tn)∥=0 and Erriter,i(t)=u(t)-ui(t) is the error of the iterative scheme.


## 3. Nonlinear Multiproduct Decomposition

We apply the abstract standard splitting schemes to the multiproduct decomposition.We have to carry out the following steps:(i)
apply the nonlinear Strang splitting scheme,(ii)
embed the Strang splitting scheme into the multiproduct expansion.To apply the abstract setting of a nonlinear Magnus expansion, we deal with the following modified nonlinear equation:(20)∂tu=A(t,u(t))u(t)+B(t,u(t))u(t),u(0)=u0,
where A(t,u),B(t,u)∈[0,T]×X are noncommuting operators and X is a general Banach space; for example, X⊂ℝm, where m is the rank of the matrices.To apply the nonlinear Magnus expansion, we deal with.(21)∂tu=F(t,u(t))u(t),u(0)=u0,u(t)=exp(ΩF(t,u0)),
where the first order Magnus operator is given by Euler's formula
(22)ΩF,1(t,u0)=∫0t‍F(s,u(s))ds=tF(0,u0),
and the second order Magnus operator is given by the midpoint rule:
(23)ΩF,2(t,u0)=∫0t‍F(s,exp(ΩF,1(s,u0)))ds,=tF(t2,exp(ΩF,1(t2,u0))),=tF(t2,t2,F(0,u0)),
to generalise this to higher order schemes; see [5].In the following, we apply the schemes for the multiproduct expansions to the time sequences,Δt=tn+1-tn and n=1,…,N.For the sequential splitting (Lie–Trotter) product we have(24)𝒯1,AB(Δt)=exp(ΩA,2(Δt,exp(ΩA,1(Δt,u0))))×exp(ΩB,1(Δt,u0)),
where ΩA,i,ΩB,i are the Magnus expansions (see (22) and (23)) of order i.Furthermore, for the Strang splitting scheme [8], we approximate the scheme based on the symmetrical splitting by
(25)𝒯2(Δt)=12(𝒯1,AB(Δt)+𝒯1,BA(Δt)).Based on the nonlinear approximations of Strang's, we generalise to higher order schemes the idea of the multiproduct expansion (MPE):(26)exp(ΩF,1(Δt,u0))=exp(ΩA,1(Δt,u0)+ΩB,1(Δt,u0))=∑k‍ck∏i‍exp(akiΩA,1(Δt,u0))×exp(bkiΩB,1(Δt,u0)),
where we assume that a linearisation with the first order integration scheme is sufficient.For the applications, we assume we have to solve time-irreversible problem, andaki and bki have to be positive.We straightforwardly follow the derivation, based on the linear multiproduct expansion; see [6].We construct𝒯S(Δt) as a product with ∑i‍aki=1 and ∑i‍bki=1. Based on the symmetry of the linearised scheme, we have
(27)𝒯S(-h)𝒯S(h)=1;
we skip the even orders based on the error truncation of the Strang splitting (see [9]) and obtain
(28)𝒯S(Δt)=exp(ΩA,1(Δt,u0)+ΩB,1(Δt,u0))+Δt3E3+Δt5E5+⋯,
where the error terms Ei are nested commutators of ΩA,1 and ΩB,1 depending on the specific form of 𝒯S.We apply𝒯S at step size Δt/k and obtain
(29)𝒯Sk(Δtk)=exp(ΩA,1(Δt,u0)+ΩB,1(Δt,u0))+k-2Δt3E3+k-4Δt5E5+⋯;
this serves as the kernel of a linearised multiproduct expansion; see [6].The first schemes are given by(30)𝒯4(Δt)=-13𝒯2(Δt)+43𝒯22(Δt2),(31)𝒯6(Δt)=124𝒯2(Δt)-1615𝒯22(Δt2)+8140𝒯23(Δt3).


### 3.1. Error Analysis of the Nonlinear Multiproduct Expansion

We deal with the following error of a linearised MPE based on the first order nonlinear Magnus expansion operators. In Theorem4, we discuss the improvement of a lower order Strang splitting scheme, given as 𝒯2 (see (25)), to a higher order scheme based on linearized multiproduct expansion (MPE) methods. The idea is based on an iterative evaluation of the extrapolation methods, see [6], which allows to receive higher accuracy.Theorem 4.
For the numerical solution of (20), one considers the linearised MPE algorithm (26) of linear order 2n. One has the following convergence result:
(32)∥Sm-exp(ΩA,1(mΔt,u0)+ΩB,1(mΔt,u0))∥≤CO(Δt2n+1),mΔt≤tend,
where S=∑i=1n‍ci𝒯2ki(Δt/ki) and C is to be chosen uniformly on bounded time intervals and independently of m and Δt.Proof.
We apply the telescoping identity and the assumption to split the linearised operators and obtain(33)(Sm-exp(ΩA,1(mΔt,u0)+ΩB,1(mΔt,u0)))=∑ν=0m-1‍Sm-ν-1(S-exp(ΩA,1(Δt,u0)+ΩB,1(Δt,u0)))×exp(ΩA,1(νΔt,u0)+ΩB,1(νΔt,u0)),
where S=∑i=1n‍ci𝒯2ki(Δt/ki).
The linearised MPE expansion can be estimated (see [6])
(34)∥∑i=1n‍ci𝒯2ki(Δtki)∥≤exp(cωΔt),
and the consistency is given by estimating the second part of (33):
(35)∥∑i=1n‍ci𝒯2ki(Δtki)-exp(ΩA,1(νΔt,u0)+ΩB,1(νΔt,u0))∥≤CO(Δt2n+1);
such a result is obtained by the linear MPE version. So for the first order linearisation, we obtain a consistent scheme with higher accuracy.Remark 5.
The results in Theorem4 combine splitting of a nonlinear equation in two operators, see (20), and extrapolating a standard splitting scheme, here given as second order Strang splitting scheme, to a higher order approach with the order 2n, while n are the stages of the extrapolation scheme. Here the novel result is based on the extension to nonlinear higher order results, while the operators can be linearized and approached with linearized MPE schemes.


## 3.1. Error Analysis of the Nonlinear Multiproduct Expansion

We deal with the following error of a linearised MPE based on the first order nonlinear Magnus expansion operators. In Theorem4, we discuss the improvement of a lower order Strang splitting scheme, given as 𝒯2 (see (25)), to a higher order scheme based on linearized multiproduct expansion (MPE) methods. The idea is based on an iterative evaluation of the extrapolation methods, see [6], which allows to receive higher accuracy.Theorem 4.
For the numerical solution of (20), one considers the linearised MPE algorithm (26) of linear order 2n. One has the following convergence result:
(32)∥Sm-exp(ΩA,1(mΔt,u0)+ΩB,1(mΔt,u0))∥≤CO(Δt2n+1),mΔt≤tend,
where S=∑i=1n‍ci𝒯2ki(Δt/ki) and C is to be chosen uniformly on bounded time intervals and independently of m and Δt.Proof.
We apply the telescoping identity and the assumption to split the linearised operators and obtain(33)(Sm-exp(ΩA,1(mΔt,u0)+ΩB,1(mΔt,u0)))=∑ν=0m-1‍Sm-ν-1(S-exp(ΩA,1(Δt,u0)+ΩB,1(Δt,u0)))×exp(ΩA,1(νΔt,u0)+ΩB,1(νΔt,u0)),
where S=∑i=1n‍ci𝒯2ki(Δt/ki).
The linearised MPE expansion can be estimated (see [6])
(34)∥∑i=1n‍ci𝒯2ki(Δtki)∥≤exp(cωΔt),
and the consistency is given by estimating the second part of (33):
(35)∥∑i=1n‍ci𝒯2ki(Δtki)-exp(ΩA,1(νΔt,u0)+ΩB,1(νΔt,u0))∥≤CO(Δt2n+1);
such a result is obtained by the linear MPE version. So for the first order linearisation, we obtain a consistent scheme with higher accuracy.Remark 5.
The results in Theorem4 combine splitting of a nonlinear equation in two operators, see (20), and extrapolating a standard splitting scheme, here given as second order Strang splitting scheme, to a higher order approach with the order 2n, while n are the stages of the extrapolation scheme. Here the novel result is based on the extension to nonlinear higher order results, while the operators can be linearized and approached with linearized MPE schemes.


## 4. Improving the Initialisation of Splitting Methods

A delicate problem in splitting methods is to achieve sufficient accuracy in the first splitting steps; see [10]. Based on the accuracy of initial starting solution, the results of all the next methods, for example MPE schemes or iterative scheme, are dependent and influenced by the order of the initialization process; see [11]. Therefore, it is important to take into account the accuracy of the initial solution. In the following, we discuss such on idea to improve the initialization process with the help of the Zassenhaus formula.The improvement of the initial solution is based on the combination of the novel commutators, that are given in Section2. Such a new combination allows to reduce the splitting error between the two novel operators A and B; see (1).In Theorem6, we discuss the improvement of the initial starting solution via the Zassenhaus exponents.Theorem 6.
One solves the initial value problem by using the method given in (5) and (6). One assumes bounded and nonlinear operators A and B.
The consistency error of theA-B splitting is 𝒪(t); then one can improve the error of the A-B splitting scheme to 𝒪(tp),p>1, by improving the starting conditions c0 by
(36)c0=(πj=2pexp(C^jtj))c0,
where the C^j are called the nonlinear Zassenhaus exponents, which are given by the novel commutator, for example, the linear case given in [12].
The local splitting error ofA-B splitting method can be expressed as follows:
(37)ρn=(exp(τn(A^+B^))(I)-exp(τnB^)(I)exp(τnA^)(I))(cspn)=C^Tτnp+1+𝒪(τnp+2),
where C^T is a function of the nonlinear Lie brackets of A^ and B^.Proof.
Let us consider the subinterval[0,t], where τ=t. Then the solution of the subproblem (5) is
(38)c*(t)=exp(tA^)(I)(c0)
and after improving the initialisation, we have
(39)c*(t)=exp(tA^)(I)(πj=2pexp(C^jtj))(I)(c0)
and so the solution of the subproblem (6) becomes
(40)c**(t)=exp(tB^)(I)exp(tA^)(I)(πj=2pexp(C^jtj))(I)c0=exp(τn(B^+A^))(I)(c0)+𝒪(tp+1)
with the help of the Zassenhaus product formula.Remark 7.
The novel result in Theorem6 is an improvement of the first splitting step to a higher order result of 𝒪(tp),p>1, where p is the number of the Zassenhaus exponents. Such an additional prestep can be done for all splitting schemes, for example, A-B splitting, Strang splitting, iterative splitting and so on. Therefore the additional prestep can be done parallel to the existing splitting schemes and can be improved at the end of the computation the splitting results to the order 𝒪(tp1+p2),p1,p2>1, where p1 is the order of the splitting approach and p2 is the order of the Zassenhaus formula.Remark 8.
For example, the second orderA-B splitting after improving the initialisation is
(41)c**(t)=((-12t2[B^,A^])exp(tB^)(I)exp(tA^)(I)×exp(-12t2[B^,A^])(I))(c0)=exp(t(B^+A^))(I)(c0)+𝒪(t3),
where the novel order of the new scheme is 𝒪(tp1+p2),p1,p2>1, with p1=1 and p2=2.A next higher order is a third orderA-B splitting after improving the initialisation is(42)c**(t)=(exp(tB^)(I)exp(tA^)(I)exp(-12t2[B,A])×(I)exp(16t3[B^,[B^,A^]](I)-13[A^,[A^,B^]](I)))(c0)=exp(t(B^+A^))(I)(c0)+𝒪(t4),
where the novel order of the new scheme is 𝒪(tp1+p2),p1,p2>1, with p1=1 and p2=3, and the commutator is given by [A,B](c)=(B′(c)∘A)(c)-(A′(c)∘B)(c).Remark 9.
The same idea can also be carried out using the Strang splitting method; see the linear case in [10]. We achieve a novel order of the new scheme to 𝒪(tp1+p2),p1,p2>1, with p1=2 and p2 being the number of Zassenhaus exponents.


## 5. Alternative Approaches: Iterative Schemes for Linearisation

We next discuss fixed-point iteration and Newton's method as alternative approaches to linearise nonlinear problems.We are to solve the nonlinear problem(43)F(x)=0,
where F:ℝn→ℝn.


### 5.1. Fixed-Point Iteration

The nonlinear equations can be formulated as fixed-point problems:(44)x=K(x),
where K is the fixed-point map and is nonlinear; for example, K(x)=x-F(x).A solution of (45) is called a fixed point of the map K.The fixed-point iteration is given by(45)xi+1=K(xi)
and is called a  nonlinear Richardson iteration, a  Picard iteration, or  the method of successive substitution.Definition 10.
Suppose thatΩ≤ℝn and G:Ω→ℝm. G is Lipschitz continuous on Ω with Lipschitz constant γ if
(46)∥G(x)-G(y)∥≤γ∥x-y∥,
for all x,y∈Ω.For the convergence, we have to assume thatK is a contraction map on Ω with Lipschitz constant γ<1.In the following, we present new results, in Algorithm14, to an application to a nonlinear Hamiltonian problem, see [3]. The model equations are numerically simulated in Section 6.Algorithm 11.
We apply the fixed-point iterative scheme to decouple the nonseparable Hamiltonian problem(47)q˙i=∂H∂p(pi,qi-1),(48)p˙i=-∂H∂q(pi-1,qi),
while we have the initial condition for the fixed-point iteration:
(49)(p0,q0)=(p(tn),q(tn)).
The initial starting solutions for theith iterative steps are given as follows: pi-1(t),qi-1(t) are the solutions of the (i-1)th iterative step and we have the initial condition for the fixed-point iteration: (p0(t),q0(t))t=(p(tn),q(tn))t.
We assume that we have convergent results afteri=1…,m iterative steps or with the stopping criterion:
(50)max(∥pi+1-pi∥,∥qi+1-qi∥)≤err,
while ∥·∥ is the Euclidean norm (or a simple vector-norm, e.g., L2) and the error bound is given, for example, as err=10-4.Example 12.
In the following, we present the application of the iterative Verlet as fixed-point scheme to the nonlinear Hamiltonian problem (47)-(48).
(
1
) We start with (q0,p0)t=(q(tn),p(tn))t.
(
2
) The iterative scheme based on the iterative Verlet, see [3], is given as
(51)qi(t)=q(tn)+h∂H∂p(p(tn)-h2∂H∂q(pi-1(t),qi-1(t)),q(tn)),(52)pi(t)=p(tn)-h2∂H∂q(pi-1(t),qi-1(t))-h2∂H∂q((p(tn)-h2∂H∂q(pi-1(t),qi-1(t))),q(tn)+h∂H∂p×(p(tn)-h2∂H∂q×(pi-1(t),qi-1(t)),q(tn)∂H∂q)),(53)fort∈[tn,tn+1],h=tn+1-tn,n=0,1,…,N,(54)i=1,2,3…,I,
where the stopping criterion is I=3 or 4 or with (62).Remark 13.
For the alternative scheme, we can also improve the initialisation process; see Section4 for the standard schemes.
To improve the initial solution we can start with the following preprocesses.
(
1
) We start with a result of the explicit Euler method:
(55)(q0,p0)t=(q(tn+1)Euler1st,p(tn+1)Euler1st)t.
(
2
) We start with a result of the explicit RK method:
(56)(q0,p0)t=(q(tn+1)RK4th,p(tn+1)RK4th)t.
Such an improvement allows to obtain higher order results and reduce the reduction of the order due to lower order initial solutions; see [11].


### 5.2. Newton's Method

We solve the nonlinear operator equation (43).Here,F:D⊂X→Y with the Banach spaces X,Y is given with the norms ∥·∥X and ∥·∥Y. Let F be at least continuously differentiable; furthermore, we assume that x0 is a starting solution of the unknown solution x*.Then thesuccessive linearisation leads to general Newton's method:
(57)F′(xi)Δxi=-F(xi),
where Δxi=xi+1-xi and i=0,1,2,….The method finds the solution to a nonlinear problem by solving a sequence of linear problems of the same kind.Algorithm 14.
We apply Newton's method to the non-separable Hamiltonian problem:(58)q˙=∂H∂p(p,q),p˙=-∂H∂q(p,q),(59)F1(p,q)=q˙-∂H∂p(p,q)=0,F2(p,q)=p˙+∂H∂q(p,q)=0,
where 0 is the zero vector.
Further, we have the initial condition, which are given as follows:(60)(p0,q0)=(p(tn),q(tn)).
The initial starting solutions for theith iterative steps are given as follows: pi-1(t),qi-1(t) are the solutions of the (i-1)th iterative step and we have the initial condition for the fixed-point iteration: (p0(t),q0(t))t=(p(tn),q(tn))t.
We apply Newton's method in its discrete form, which is given as(61)qi=qi-1-D(F1(pi-1,qi-1))-1(F1(pi-1,qi-1)),withD(F1(pi-1,qi-1))=∂∂qi-1(q˙i-1-∂H∂p(pi-1,qi-1)),pi=pi-1-D(F2(pi-1,qi-1))-1(F2(pi-1,qi-1)),withD(F2(pi-1,qi-1))=∂∂pi-1(p˙i-1-∂H∂q(pi-1,qi-1)),i=1,2,3,…,I,
where D is the Jacobian matrix (see [13]) and pi-1(t),qi-1(t) are the given solutions of the (i-1)th iterative step.
We assume that we have convergent results afteri=1…,I iterative steps or with the stopping criterion:
(62)max(∥pi-pi-1∥,∥qi-qi-1∥)≤err,
while ∥·∥ is the Euclidean norm (or a simple vector-norm, e.g., L2) and the error bound is given, for example, as err=10-4.Example 15.
We apply the iterative Verlet scheme to the underlying fixed-point problem (59) of nonlinear Hamiltonian problem (58).
(
1
) We start with (q0,p0)t=(q(tn),p(tn))t.
(
2
) The iterative scheme based on the iterative Verlet (see [3]) is given in (51)–(53).Remark 16.
For the improvement method, we can apply the weighted Newton method. We try to skip the delicate outer diagonals in the Jacobian matrixD and apply
(63)qi=qi-1-D(F1(pi-1,qi-1)+δ1qi-1)-1×(F1(pi-1,qi-1)+ϵ1qi-1),pi=pi-1-D(F2(pi-1,qi-1)+δ2pi-1)-1×(F2(pi-1,qi-1)+ϵ2pi-1),i=1,2,3,…,I,
where the function δ1,δ2 can be applied as a scalar, for example δ=10-6, and the same with ϵ1,ϵ2; see [13]. It is important to ensure that δ is small enough to preserve the convergence.Remark 17.
Linearisation methods can be applied to MPE methods. Nonseparable Hamiltonian problems can be decoupled into separable Hamiltonian problems; see the ideas in [14].


## 5.1. Fixed-Point Iteration

The nonlinear equations can be formulated as fixed-point problems:(44)x=K(x),
where K is the fixed-point map and is nonlinear; for example, K(x)=x-F(x).A solution of (45) is called a fixed point of the map K.The fixed-point iteration is given by(45)xi+1=K(xi)
and is called a  nonlinear Richardson iteration, a  Picard iteration, or  the method of successive substitution.Definition 10.
Suppose thatΩ≤ℝn and G:Ω→ℝm. G is Lipschitz continuous on Ω with Lipschitz constant γ if
(46)∥G(x)-G(y)∥≤γ∥x-y∥,
for all x,y∈Ω.For the convergence, we have to assume thatK is a contraction map on Ω with Lipschitz constant γ<1.In the following, we present new results, in Algorithm14, to an application to a nonlinear Hamiltonian problem, see [3]. The model equations are numerically simulated in Section 6.Algorithm 11.
We apply the fixed-point iterative scheme to decouple the nonseparable Hamiltonian problem(47)q˙i=∂H∂p(pi,qi-1),(48)p˙i=-∂H∂q(pi-1,qi),
while we have the initial condition for the fixed-point iteration:
(49)(p0,q0)=(p(tn),q(tn)).
The initial starting solutions for theith iterative steps are given as follows: pi-1(t),qi-1(t) are the solutions of the (i-1)th iterative step and we have the initial condition for the fixed-point iteration: (p0(t),q0(t))t=(p(tn),q(tn))t.
We assume that we have convergent results afteri=1…,m iterative steps or with the stopping criterion:
(50)max(∥pi+1-pi∥,∥qi+1-qi∥)≤err,
while ∥·∥ is the Euclidean norm (or a simple vector-norm, e.g., L2) and the error bound is given, for example, as err=10-4.Example 12.
In the following, we present the application of the iterative Verlet as fixed-point scheme to the nonlinear Hamiltonian problem (47)-(48).
(
1
) We start with (q0,p0)t=(q(tn),p(tn))t.
(
2
) The iterative scheme based on the iterative Verlet, see [3], is given as
(51)qi(t)=q(tn)+h∂H∂p(p(tn)-h2∂H∂q(pi-1(t),qi-1(t)),q(tn)),(52)pi(t)=p(tn)-h2∂H∂q(pi-1(t),qi-1(t))-h2∂H∂q((p(tn)-h2∂H∂q(pi-1(t),qi-1(t))),q(tn)+h∂H∂p×(p(tn)-h2∂H∂q×(pi-1(t),qi-1(t)),q(tn)∂H∂q)),(53)fort∈[tn,tn+1],h=tn+1-tn,n=0,1,…,N,(54)i=1,2,3…,I,
where the stopping criterion is I=3 or 4 or with (62).Remark 13.
For the alternative scheme, we can also improve the initialisation process; see Section4 for the standard schemes.
To improve the initial solution we can start with the following preprocesses.
(
1
) We start with a result of the explicit Euler method:
(55)(q0,p0)t=(q(tn+1)Euler1st,p(tn+1)Euler1st)t.
(
2
) We start with a result of the explicit RK method:
(56)(q0,p0)t=(q(tn+1)RK4th,p(tn+1)RK4th)t.
Such an improvement allows to obtain higher order results and reduce the reduction of the order due to lower order initial solutions; see [11].


## 5.2. Newton's Method

We solve the nonlinear operator equation (43).Here,F:D⊂X→Y with the Banach spaces X,Y is given with the norms ∥·∥X and ∥·∥Y. Let F be at least continuously differentiable; furthermore, we assume that x0 is a starting solution of the unknown solution x*.Then thesuccessive linearisation leads to general Newton's method:
(57)F′(xi)Δxi=-F(xi),
where Δxi=xi+1-xi and i=0,1,2,….The method finds the solution to a nonlinear problem by solving a sequence of linear problems of the same kind.Algorithm 14.
We apply Newton's method to the non-separable Hamiltonian problem:(58)q˙=∂H∂p(p,q),p˙=-∂H∂q(p,q),(59)F1(p,q)=q˙-∂H∂p(p,q)=0,F2(p,q)=p˙+∂H∂q(p,q)=0,
where 0 is the zero vector.
Further, we have the initial condition, which are given as follows:(60)(p0,q0)=(p(tn),q(tn)).
The initial starting solutions for theith iterative steps are given as follows: pi-1(t),qi-1(t) are the solutions of the (i-1)th iterative step and we have the initial condition for the fixed-point iteration: (p0(t),q0(t))t=(p(tn),q(tn))t.
We apply Newton's method in its discrete form, which is given as(61)qi=qi-1-D(F1(pi-1,qi-1))-1(F1(pi-1,qi-1)),withD(F1(pi-1,qi-1))=∂∂qi-1(q˙i-1-∂H∂p(pi-1,qi-1)),pi=pi-1-D(F2(pi-1,qi-1))-1(F2(pi-1,qi-1)),withD(F2(pi-1,qi-1))=∂∂pi-1(p˙i-1-∂H∂q(pi-1,qi-1)),i=1,2,3,…,I,
where D is the Jacobian matrix (see [13]) and pi-1(t),qi-1(t) are the given solutions of the (i-1)th iterative step.
We assume that we have convergent results afteri=1…,I iterative steps or with the stopping criterion:
(62)max(∥pi-pi-1∥,∥qi-qi-1∥)≤err,
while ∥·∥ is the Euclidean norm (or a simple vector-norm, e.g., L2) and the error bound is given, for example, as err=10-4.Example 15.
We apply the iterative Verlet scheme to the underlying fixed-point problem (59) of nonlinear Hamiltonian problem (58).
(
1
) We start with (q0,p0)t=(q(tn),p(tn))t.
(
2
) The iterative scheme based on the iterative Verlet (see [3]) is given in (51)–(53).Remark 16.
For the improvement method, we can apply the weighted Newton method. We try to skip the delicate outer diagonals in the Jacobian matrixD and apply
(63)qi=qi-1-D(F1(pi-1,qi-1)+δ1qi-1)-1×(F1(pi-1,qi-1)+ϵ1qi-1),pi=pi-1-D(F2(pi-1,qi-1)+δ2pi-1)-1×(F2(pi-1,qi-1)+ϵ2pi-1),i=1,2,3,…,I,
where the function δ1,δ2 can be applied as a scalar, for example δ=10-6, and the same with ϵ1,ϵ2; see [13]. It is important to ensure that δ is small enough to preserve the convergence.Remark 17.
Linearisation methods can be applied to MPE methods. Nonseparable Hamiltonian problems can be decoupled into separable Hamiltonian problems; see the ideas in [14].


## 6. Numerical Examples

In this section, we treat experiments to verify the benefit of our methods: we treat nonseparable Hamiltonian with a nonlinear kinetic and nonlinear potential.The motivation arose from simulating a Levitron.A Levitron is described on the basis of rigid body theory; see the convention of [15] for the Euler angles. The angular velocity ωϕ is along the z-axis of the system, ωθ is along the line of the nodes, and ωψ is along the z′-axis. We transform into body coordinates and obtain
(64)ω=(ϕ˙sinθsinψ+θ˙cosψϕ˙sinθcosψ+θ˙sinψϕ˙cosθ+ψ˙).The kinetic energy can be written as(65)T=12[m(x˙2+y˙2+z˙2)+A(θ˙2+ϕ˙2sin2θ)+C(ψ˙+ϕ˙cosθ)2].
The potential energy U is given by the sum of the gravitational energy and the interaction potential of the Levitron in the magnetic field of the base plate and is
(66)U=mgz-μ(sinψsinθΦx+cosψsinθΦy+cosθΦz)
with mu the magnetic moment of the top and Φ the magnetostatic potential. Furthermore, we use the notation of [16] and the potential of a ring dipole as an approximation for a magnetised plane with a centred unmagnetised hole. Furthermore, we introduced a nondimensionalisation for the variables and the magnetostatic potential:
(67)Ψ=Z(1+Z2)3/2-(X2+Y2)34(2Z2-3)Z(1+Z2)7/2.
The lengths were scaled by the radius R of the base plane, mass is measured in units of m, and energy in units of mgh. Therefore the one time unit is R/g. So the dimensionless Hamiltonian with q=(X,Y,Z,θ,ψ,ϕ) is given by
(68)H=12(p12+p22+p32+p42a+(p5-P-6cosq4)2asin2q4+p6c)-M[sinq4(cosq5∂Ψ∂q1sinq5∂Ψ∂q2)+cosq4∂Ψ∂q3]+q3
with a and c being the nondimensionalised inertial parameters and M being the ratio of the gravitational and magnetic energy.For the nonseparable Hamiltonian of (68), we have(69)q˙=∂H∂p(p,q)=(p1,p2,p3,p4a,(p5-p6cosq4)2asin2q4,p6(cos2q4+(a/c)sin42)-p5cosq4asin2q4).The same is given for(70)p˙=-∂H∂q(p,q)=(M(sinq4cosq5∂2Ψ∂q12+cosq4∂2Ψ∂q1∂q3),M(sinq4cosq5∂2Ψ∂q22+cosq4∂2Ψ∂q2∂q3),M(sinq4(sinq5∂2Ψ∂q2∂q3+cosq5∂2Ψ∂q1∂q3)+cosq4∂2Ψ∂q32)-1,M(cosq4(sinq5Ψq2+cosq5∂Ψ∂q1)-sinq4∂Ψ∂q3)-p6(p5-p6cosq4)asinq4-cosq4(p5-p6cosq4)2asin3q4,M(sinq4(cosq5∂Ψ∂q2-sinq5∂Ψ∂q1)),0).A and B are linearised Lie operators and their vector fields are given by
(71)A=∂H∂p·∂∂q,B=-∂H∂q·∂∂p.The transfer to the operators is given in the following description.The exponential operatorsehA and ehB are then just shift operators, with 𝒯2(h) being a symmetric second order splitting method:
(72)𝒯2,VV(h)(Δt)=e(Δt/2)BeΔtAe(Δt/2)B.
This corresponds to the velocity form of the Verlet algorithm (VV).Furthermore, the splitting scheme(73)𝒯2,PV(h)(Δt)=e(Δt/2)AeΔtBe(Δt/2)A
corresponds to the position form of the Verlet algorithm (PV).Solving the equations of motion (69) and (70) numerically with nonlinear MPE methods based on the Verlet algorithm, we plot the movement of the centre of mass as shown in Figure 1. The axes in the plot show the nondimensional variables X, Y, and Z. The trajectory starts at the equilibrium point (q1,q2,q3)=(0,0,1.72).Trajectory calculated with MPE method supplied with the Verlet algorithm.
(a)
(b)The linearisation is done with two and four iterations per time-step, to see whether how many iterations are reasonable. The results are shown in Figure2.Errors of the numerical scheme for MPE method supplied with the iterative Verlet method.
(a)
(b)In a first comparison, we deal with the second order Verlet algorithm and improve the scheme with iterative steps. The tables should give an impression of the timescales of the problem and the errors; see Table1.Table 1
Numerical results for MPE method supplied with the iterative Verlet algorithm.
Verlet
Timestep
1
0
-
6
1
0
-
6
1
0
-
6
Iterations per step
1
2
4
Stability
Ok
Ok
Ok
Computing time
67 min
120 min
219 min
Mean error
0.068
0.068
0.068
Maximal error
0.0187
0.0188
0.0187Remark 18.
The iterations improve the algorithm, when we apply sufficiently accurate initial solutions. We also obtain a benefit in reducing the computational time instead of applying only expensive Runge-Kutta methods.We improve the solution with an extrapolation scheme with fourth order. We can see the errors this algorithm produces in comparison with the Runge-Kutta method with small time-steps (10-5 time units per step). Figure 3 presents the results of the 4th order MPE method with different time-steps and compares it with the Runge-Kutta solution.Figure 3
Errors of the numerical scheme: 4th order extrapolation scheme with the Verlet method and kernel (h=10-5).Also we tested the 6th order MPE method with different time-steps and compared it with the Runge-Kutta method; see Figure4.Figure 4
Errors of the numerical scheme: 6th order extrapolation scheme with Verlet method as kernel (h=10-5).Table2 gives an impression of the standard higher order Runge-Kutta solver (4th order scheme) with the timescales of the problem and the errors.Table 2
Stability and computational time with 4th order explicit Runge-Kutta.
Runge-Kutta
Timestep
1
0
-
5
Number of steps
1
0
9
Computing time
119 min
Stability
OkWe could improve the accuracy and computational time needed by using higher order extrapolation schemes; see Table3.Table 3
Errors and computational time with 4th and 6th order MPE schemes using Verlet scheme as kernel.
Extrapolation 4th order
Extrapolation 6th order
Timestep
1
0
-
5
1
0
-
6
1
0
-
5
1
0
-
6
Number of steps
1
0
8
1
0
9
1
0
8
1
0
9
Computing time
14 min
142 min
29 min
272 min
Mean error
0.007
0.007
0.0068
0.0068
Maximal error
0.0226
0.0234
0.0188
0.0188Remark 19.
We improved the basic splitting schemes with the Verlet algorithm and extrapolation schemes. At least more accurate solutions were achieved and computational time was saved. The best result was achieved with order 10 andh=10-3, and for such a case we could improve the results and obtain them faster than with standard RK schemes.Remark 20.
We were able to accelerate the computation with extrapolation methods and increase the time-steps. Because of the accuracy of the underlying kernel, which is an iterative Verlet method withi=1, we were able to remain at order 8 and only 10-3 as the time-step.


## 7. Conclusions and Discussion

We presented novel MPE approaches to nonlinear differential equations. Based on the ideas of the nonlinear Magnus expansion and linearisation schemes, we derived a linearised MPE approach possessing higher accuracy. We discussed the numerical errors and different improvements of the underlying splitting methods. Numerical examples confirmed the application to nonlinear equations. In the future, we will focus on the development of improved MPE methods for more highly nonlinear approaches.



---

*Source: 681575-2013-09-18.xml*