Nonlinear Extension of Multiproduct Expansion Schemes and Applications to Rigid Bodies

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.

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, 𝜕_𝑡𝑢=𝐹(𝑢(𝑡))=𝐴(𝑢(𝑡))+𝐵(𝑢(𝑡)),𝑢(0)=𝑢₀, (1) with unbounded operators 𝐴∶𝐷(𝐴)⊂𝐗→𝐗 and 𝐵∶𝐷(𝐵)⊂𝐗→𝐗. We have further 𝐹(𝑣)=𝐴(𝑣)+𝐵(𝑣), 𝑣∈𝐷(𝐴)∩𝐷(𝐵).

We assume there are suitable chosen subspaces of the underlying Banach space (𝑋,‖⋅‖_𝑋) such that 𝐷(𝐹)=𝐷(𝐴)∩𝐷(𝐵)≠∅.

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 𝑢(𝑡)=ℰ_𝐹(𝑡,𝑢(0)),0≤𝑡≤𝑇, (2) with the evolution operator ℰ_𝐹 depending on the actual time 𝑡 and the initial value 𝑢(0). We use a formal notation: 𝑢(𝑡)=exp𝑡𝐷_𝐹𝑢(0),0≤𝑡≤𝑇. (3) Here the evolution operator exp(𝑡𝐷_𝐹) and the Lie derivative 𝐷_𝐹 associated with 𝐹 are exp𝑡𝐷_𝐹𝐺𝑣=𝐺ℰ_𝐹()𝐷𝑡,𝑣,0≤𝑡≤𝑇,_𝐹𝐺𝑣=𝐺^′(𝑉)𝐹(𝑣) (4) for any unbounded nonlinear operator 𝐺∶𝐷(𝐺)⊂𝑋→𝑋 with Frechet’s derivative 𝐺^′.

The paper is outlined as follows. In Section 2, 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 𝜕𝑐^∗(𝑡)𝑐𝜕𝑡=𝐴^∗(𝑡)with𝑡∈𝑡^𝑛,𝑡^𝑛+1,𝑐^∗𝑡^𝑛=𝑐^𝑛_𝑠𝑝,𝜕𝑐^∗∗(𝑡)𝑐𝜕𝑡=𝐵^∗∗(𝑡)with𝑡∈𝑡^𝑛,𝑡^𝑛+1,𝑐^∗∗𝑡^𝑛=𝑐^∗𝑡^𝑛+1, (5)(6)for 𝑛=0,1,…,𝑁−1, whereby 𝑐^𝑛_𝑠𝑝=𝑐₀ is given from (1). The approximated split solution at the point 𝑡=𝑡^𝑛+1 is defined by 𝑐^𝑛+1_𝑠𝑝=𝑐^∗∗(𝑡^𝑛+1).(2) The Strang-Marchuk operator splitting: A-B-A splitting 𝜕𝑐^∗(𝑡)𝑐𝜕𝑡=𝐴^∗(𝑡)with𝑡∈𝑡^𝑛,𝑡^𝑛+1/2,𝑐^∗𝑡^𝑛=𝑐^𝑛_𝑠𝑝,𝜕𝑐^∗∗(𝑡)𝑐𝜕𝑡=𝐵^∗∗(𝑡)with𝑡∈𝑡^𝑛,𝑡^𝑛+1/2,𝑐^∗∗𝑡^𝑛=𝑐^∗𝑡^𝑛+1/2,𝜕𝑐^∗∗∗(𝑡)𝑐𝜕𝑡=𝐴^∗(𝑡)with𝑡∈𝑡^𝑛+1/2,𝑡^𝑛+1,𝑐^∗∗∗𝑡^𝑛+1/2=𝑐^∗∗𝑡^𝑛+1, (7)where 𝑡^𝑛+1/2=𝑡^𝑛+𝜏/2, 𝜏 is the local time-step. The approximate split solution at the point 𝑡=𝑡^𝑛+1 is defined by 𝑐^𝑛+1_𝑠𝑝=𝑐^∗∗∗(𝑡^𝑛+1).(3) Iterative splitting with respect to one operator 𝜕𝑐_𝑖(𝑡)𝑐𝜕𝑡=𝐴_𝑖(𝑡)+𝐵𝑐_𝑖−1(𝑡),with𝑐_𝑖𝑡^𝑛=𝑐^𝑛,𝑖=1,2,…,𝑚, (8)where 𝑐₀(𝑡^𝑛)=𝑐^𝑛,𝑐₋₁=0, and 𝑐^𝑛 is the known split approximation at the time level 𝑡=𝑡^𝑛. The split approximation at the time level 𝑡=𝑡^𝑛+1 is defined by 𝑐^𝑛+1=𝑐_𝑚+1(𝑡^𝑛+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 operator 𝐹 which is a linear Lie operator.

Definition 1. For the given operator 𝐹 we associate a new operator, denoted by 𝐹.
This Lie operator acts on the space of the differentiable operators of the type 𝑆⟶𝑆 (9) and maps each operator 𝐺 into the new operator 𝐹(𝐺), such that for any element 𝑐∈𝑆𝐹(𝐺)(𝑐)=𝐺^′(𝑐)∘𝐹(𝑐). (10)

Here the derivatives are given by the following.

Definition 2. The 𝑘th power of the Lie operator 𝐹 applied to some operator 𝐺 can be expressed as the 𝑘th derivative of 𝐺; that is, the relation 𝐹^𝑘(𝐺𝑐)=𝜕)(^𝑘𝐺(𝑐)𝜕𝑡^𝑘 (11) is valid for all 𝑘=1,2,….

We deal with the problem (1) and split the operator 𝐹 into the sum 𝐴+𝐵; then we have the following results of the abstract splitting errors. We have an appropriate norm ‖⋅‖ in a general Banach space 𝐗.

Theorem 3. The splitting error for the A-B splitting is given by Err_{𝐴-𝐵,1,2}=𝜏𝐵^′(𝑐)∘𝐴(𝑐)−𝐴^′(𝑐)(𝑐)𝜏∘𝐵+𝑂², (12)
The splitting error of the Strang splitting is given by Err_Strang,1,2=1𝜏24²[𝑐)[[𝑐𝜏([𝐵,𝐵,𝐴]](−2𝐴,𝐴,𝐵]]())+𝑂⁴. (13) with the commutator [𝐴,𝐵](𝑐)=(𝐵^′(𝑐)∘𝐴)(𝑐)−(𝐴^′(𝑐)∘𝐵)(𝑐).
The splitting error of the iterative splitting is given by Err_iter,𝑖(𝑡)=𝜏_𝑛𝐹₂(𝐼)Err_{iter,𝑖−1}(𝑡),for𝜏_𝑛=𝑡−𝑡^𝑛, (14) where one assumes that the previous error ‖Err_iter(𝑡^𝑛)‖=0.

Proof. The splitting error for the A-B splitting is given by Err_𝐴-𝐵=exp𝜏𝐴+𝐵(𝐼)−exp𝜏𝐴exp𝜏𝐵(𝐼)(𝑐). (15) We evaluate the exp-function and apply the commutator [𝐴,𝐵](𝑐) and obtain Err_{𝐴-𝐵,1,2}=𝜏𝐵^′(𝑐)∘𝐴(𝑐)−𝐴^′(𝑐)∘𝐵(𝑐)+𝑂𝜏². (16)
The splitting error of the Strang Splitting is given by Err_Strang𝜏𝐵(𝐼)−𝜏=exp𝐴+exp2𝐴𝜏𝐵𝜏expexp2𝐴(𝐼)(𝑐). (17) We evaluate the exp-function and apply twice the commutator [𝐴,𝐵](𝑐) based on the symmetry of the scheme and obtain Err_Strang,1,2=1𝜏24²[𝑐)[[𝑐𝜏([𝐵,𝐵,𝐴]](−2𝐴,𝐴,𝐵]]())+𝑂⁴, (18) where the commutator is given as [𝐴,𝐵](𝑐)=(𝐵^′(𝑐)∘𝐴)(𝑐)−(𝐴^′(𝑐)∘𝐵)(𝑐).
The splitting error of the iterative splitting is given by Err_iter,𝑖(𝑡=𝜏𝐴(𝐼)𝑡)expErr^𝑛+^𝑡_𝑡𝑛exp(𝑡−𝑠)𝐴(𝐼)𝐵(𝐼)Err_{iter,𝑖−1}(𝑠)𝑑𝑠,≤𝜏_𝑛𝐵(𝐼)Err_{iter,𝑖−1}(𝑡),for𝜏_𝑛=𝑡−𝑡^𝑛, (19) where we have applied the method of variation of constants to the linearised operators. Furthermore, we assume that the previous error ‖Err_iter(𝑡^𝑛)‖=0 and Err_iter,𝑖(𝑡)=𝑢(𝑡)−𝑢_𝑖(𝑡) 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: 𝜕_𝑡𝑢=𝐴(𝑡,𝑢(𝑡))𝑢(𝑡)+𝐵(𝑡,𝑢(𝑡))𝑢(𝑡),𝑢(0)=𝑢₀, (20) where 𝐴(𝑡,𝑢),𝐵(𝑡,𝑢)∈[0,𝑇]×𝐗 are noncommuting operators and 𝐗 is a general Banach space; for example, 𝐗⊂ℝ^𝑚, where 𝑚 is the rank of the matrices.

To apply the nonlinear Magnus expansion, we deal with. 𝜕_𝑡((𝑡𝑢(𝑡),𝑢(0)𝑢=𝐹𝑡,𝑢))=𝑢₀,𝑢(𝑡)=expΩ_𝐹𝑡,𝑢₀, (21) where the first order Magnus operator is given by Euler's formula Ω_𝐹,1𝑡,𝑢₀=^𝑡₀𝐹(𝑠,𝑢(𝑠))𝑑𝑠=𝑡𝐹0,𝑢₀, (22) and the second order Magnus operator is given by the midpoint rule: Ω_𝐹,2𝑡,𝑢₀=^𝑡₀𝐹𝑠,expΩ_𝐹,1𝑠,𝑢₀𝑡𝑑𝑠,=𝑡𝐹2,expΩ_𝐹,1𝑡2,𝑢₀,𝑡=𝑡𝐹2,𝑡2,𝐹0,𝑢₀, (23) to generalise this to higher order schemes; see [5].

In the following, we apply the schemes for the multiproduct expansions to the time sequences, Δ𝑡=𝑡^𝑛+1−𝑡^𝑛 and 𝑛=1,…,𝑁.

For the sequential splitting (Lie–Trotter) product we have 𝒯_1,𝐴𝐵(Δ𝑡)=expΩ_𝐴,2Δ𝑡,expΩ_𝐴,1Δ𝑡,𝑢₀×expΩ_𝐵,1Δ𝑡,𝑢₀, (24) where Ω_𝐴,𝑖,Ω_𝐵,𝑖 are the Magnus expansions (see (22) and (23)) of order 𝑖.

Furthermore, for the Strang splitting scheme [8], we approximate the scheme based on the symmetrical splitting by 𝒯₂()=1Δ𝑡2𝒯_1,𝐴𝐵(Δ𝑡)+𝒯_1,𝐵𝐴(Δ𝑡). (25)

Based on the nonlinear approximations of Strang's, we generalise to higher order schemes the idea of the multiproduct expansion (MPE): expΩ_𝐹,1Δ𝑡,𝑢₀=expΩ_𝐴,1Δ𝑡,𝑢₀+Ω_𝐵,1Δ𝑡,𝑢₀=_𝑘𝑐_𝑘𝑖exp𝑎_𝑘𝑖Ω_𝐴,1Δ𝑡,𝑢₀×exp𝑏_𝑘𝑖Ω_𝐵,1Δ𝑡,𝑢₀, (26) 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, and 𝑎_𝑘𝑖 and 𝑏_𝑘𝑖 have to be positive.

We straightforwardly follow the derivation, based on the linear multiproduct expansion; see [6].

We construct 𝒯_𝑆(Δ𝑡) as a product with ∑_𝑖𝑎_𝑘𝑖=1 and ∑_𝑖𝑏_𝑘𝑖=1. Based on the symmetry of the linearised scheme, we have 𝒯_𝑆(−ℎ)𝒯_𝑆(ℎ)=1; (27) we skip the even orders based on the error truncation of the Strang splitting (see [9]) and obtain 𝒯_𝑆(Δ𝑡)=expΩ_𝐴,1Δ𝑡,𝑢₀+Ω_𝐵,1Δ𝑡,𝑢₀+Δ𝑡³𝐸₃+Δ𝑡⁵𝐸₅+⋯, (28) where the error terms 𝐸_𝑖 are nested commutators of Ω_𝐴,1 and Ω_𝐵,1 depending on the specific form of 𝒯_𝑆.

We apply 𝒯_𝑆 at step size Δ𝑡/𝑘 and obtain 𝒯^𝑘_𝑆Δ𝑡𝑘=expΩ_𝐴,1Δ𝑡,𝑢₀+Ω_𝐵,1Δ𝑡,𝑢₀+𝑘⁻²Δ𝑡³𝐸₃+𝑘⁻⁴Δ𝑡⁵𝐸₅+⋯; (29) this serves as the kernel of a linearised multiproduct expansion; see [6].

The first schemes are given by 𝒯₄()1Δ𝑡=−3𝒯₂()+4Δ𝑡3𝒯²₂Δ𝑡2,𝒯₆()=1Δ𝑡𝒯24²(Δ𝑡)−16𝒯15²²Δ𝑡2+81𝒯40³²Δ𝑡3. (30)(31)

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 Theorem 4, we discuss the improvement of a lower order Strang splitting scheme, given as 𝒯₂ (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 2𝑛. One has the following convergence result: 𝑆^𝑚−expΩ_𝐴,1𝑚Δ𝑡,𝑢₀+Ω_𝐵,1𝑚Δ𝑡,𝑢₀≤𝐶𝑂Δ𝑡^2𝑛+1,𝑚Δ𝑡≤𝑡_end, (32) where 𝑆=∑^𝑛_𝑖=1𝑐_𝑖𝒯^𝑘𝑖₂(Δ𝑡/𝑘_𝑖) and 𝐶 is to be chosen uniformly on bounded time intervals and independently of 𝑚 and Δ𝑡.

Proof. We apply the telescoping identity and the assumption to split the linearised operators and obtain 𝑆^𝑚−expΩ_𝐴,1𝑚Δ𝑡,𝑢₀+Ω_𝐵,1𝑚Δ𝑡,𝑢₀=^𝑚−1_𝜈=0𝑆^{𝑚−𝜈−1}𝑆−expΩ_𝐴,1Δ𝑡,𝑢₀+Ω_𝐵,1Δ𝑡,𝑢₀×expΩ_𝐴,1𝜈Δ𝑡,𝑢₀+Ω_𝐵,1𝜈Δ𝑡,𝑢₀, (33) where 𝑆=∑^𝑛_𝑖=1𝑐_𝑖𝒯^𝑘𝑖₂(Δ𝑡/𝑘_𝑖).
The linearised MPE expansion can be estimated (see [6]) ^𝑛_𝑖=1𝑐_𝑖𝒯^𝑘𝑖₂Δ𝑡𝑘_𝑖≤exp(𝑐𝜔Δ𝑡), (34) and the consistency is given by estimating the second part of (33): ^𝑛_𝑖=1𝑐_𝑖𝒯^𝑘𝑖₂Δ𝑡𝑘_𝑖−expΩ_𝐴,1𝜈Δ𝑡,𝑢₀+Ω_𝐵,1𝜈Δ𝑡,𝑢₀≤𝐶𝑂Δ𝑡^2𝑛+1; (35) 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 Theorem 4 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 2𝑛, while 𝑛 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 Section 2. Such a new combination allows to reduce the splitting error between the two novel operators 𝐴 and 𝐵; see (1).

In Theorem 6, 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 𝐴 and 𝐵.
The consistency error of the A-B splitting is 𝒪(𝑡); then one can improve the error of the A-B splitting scheme to 𝒪(𝑡^𝑝),𝑝>1, by improving the starting conditions 𝑐₀ by 𝑐₀=𝜋^𝑝_𝑗=2exp𝐶_𝑗𝑡^𝑗𝑐₀, (36) where the 𝐶_𝑗 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 of A-B splitting method can be expressed as follows: 𝜌_𝑛=exp𝜏_𝑛𝐴+𝐵(𝐼)−exp𝜏_𝑛𝐵(𝐼)exp𝜏_𝑛𝐴(𝐼)𝑐^𝑛_𝑠𝑝=𝐶_𝑇𝜏_𝑛^𝑝+1+𝒪𝜏_𝑛^𝑝+2, (37) where 𝐶_𝑇 is a function of the nonlinear Lie brackets of 𝐴 and 𝐵.

Proof. Let us consider the subinterval [0,𝑡], where 𝜏=𝑡. Then the solution of the subproblem (5) is 𝑐^∗(𝑡)=exp𝑡𝐴(𝐼)𝑐₀ (38) and after improving the initialisation, we have 𝑐^∗(𝑡)=exp𝑡𝐴(𝐼)𝜋^𝑝_𝑗=2exp𝐶_𝑗𝑡^𝑗(𝐼)𝑐₀ (39) and so the solution of the subproblem (6) becomes 𝑐^∗∗(𝑡)=exp𝑡𝐵(𝐼)exp𝑡𝐴(𝐼)𝜋^𝑝_𝑗=2exp𝐶_𝑗𝑡^𝑗(𝐼)𝑐₀=exp𝜏_𝑛𝐵+𝐴(𝐼)𝑐₀+𝒪𝑡^𝑝+1 (40) with the help of the Zassenhaus product formula.

Remark 7. The novel result in Theorem 6 is an improvement of the first splitting step to a higher order result of 𝒪(𝑡^𝑝),𝑝>1, where 𝑝 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 𝒪(𝑡^𝑝1+𝑝2),𝑝₁,𝑝₂>1, where 𝑝₁ is the order of the splitting approach and 𝑝₂ is the order of the Zassenhaus formula.

Remark 8. For example, the second order A-B splitting after improving the initialisation is 𝑐^∗∗(𝑡)=𝑡𝐵(𝐼)𝑡𝐴(𝐼)−1expexp×exp2𝑡²𝐵,𝐴(𝐼)𝑐₀=exp𝑡𝐵+𝐴(𝐼)𝑐₀+𝒪𝑡³, (41) where the novel order of the new scheme is 𝒪(𝑡^𝑝1+𝑝2),𝑝₁,𝑝₂>1, with 𝑝₁=1 and 𝑝₂=2.

A next higher order is a third order A-B splitting after improving the initialisation is𝑐^∗∗(𝑡)=𝑡𝐵(𝐼)𝑡𝐴(𝐼)−1expexpexp2𝑡²[×(𝐼)1𝐵,𝐴]exp6𝑡³𝐴𝐼)−1𝐵,𝐵,(3𝐴,𝐴,𝐵(𝐼)𝑐₀=exp𝑡𝐵+𝐴(𝐼)𝑐₀+𝒪𝑡⁴, (42) where the novel order of the new scheme is 𝒪(𝑡^𝑝1+𝑝2),𝑝₁,𝑝₂>1, with 𝑝₁=1 and 𝑝₂=3, and the commutator is given by [𝐴,𝐵](𝑐)=(𝐵^′(𝑐)∘𝐴)(𝑐)−(𝐴^′(𝑐)∘𝐵)(𝑐).

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 𝒪(𝑡^𝑝1+𝑝2),𝑝₁,𝑝₂>1, with 𝑝₁=2 and 𝑝₂ 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 𝐹(𝑥)=0, (43) where 𝐹∶ℝ^𝑛→ℝ^𝑛.

5.1. Fixed-Point Iteration

The nonlinear equations can be formulated as fixed-point problems: 𝑥=𝐾(𝑥), (44) where 𝐾 is the fixed-point map and is nonlinear; for example, 𝐾(𝑥)=𝑥−𝐹(𝑥).

A solution of (45) is called a fixed point of the map 𝐾.

The fixed-point iteration is given by 𝑥_𝑖+1=𝐾𝑥_𝑖 (45) and is called a nonlinear Richardson iteration, a Picard iteration, or the method of successive substitution.

Definition 10. Suppose that Ω≤ℝ^𝑛 and 𝐺∶Ω→ℝ^𝑚. 𝐺 is Lipschitz continuous on Ω with Lipschitz constant 𝛾 if 𝐺(𝑥)−𝐺𝑦≤𝛾𝑥−𝑦, (46) for all 𝑥,𝑦∈Ω.

For the convergence, we have to assume that 𝐾 is a contraction map on Ω with Lipschitz constant 𝛾<1.

In the following, we present new results, in Algorithm 14, 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 𝐪̇_𝑖=𝜕𝐻𝐩𝜕𝐩^𝑖,𝐪_𝑖−1,𝐩̇_𝑖=−𝜕𝐻𝐩𝜕𝐪^𝑖−1,𝐪_𝑖, (47)(48) while we have the initial condition for the fixed-point iteration: 𝐩₀,𝐪₀=𝐩𝑡^𝑛,𝐪𝑡^𝑛. (49)
The initial starting solutions for the ith iterative steps are given as follows: 𝐩_𝑖−1(𝑡),𝐪_𝑖−1(𝑡) are the solutions of the (𝑖−1)th iterative step and we have the initial condition for the fixed-point iteration: (𝐩₀(𝑡),𝐪₀(𝑡))^𝑡=(𝐩(𝑡^𝑛),𝐪(𝑡^𝑛))^𝑡.
We assume that we have convergent results after 𝑖=1…,𝑚 iterative steps or with the stopping criterion: max𝐩_𝑖+1−𝐩_𝑖,𝐪_𝑖+1−𝐪_𝑖≤err, (50) while ‖⋅‖ is the Euclidean norm (or a simple vector-norm, e.g., 𝐿₂) and the error bound is given, for example, as err=10⁻⁴.

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 (𝐪₀,𝐩₀)^𝑡=(𝐪(𝑡^𝑛),𝐩(𝑡^𝑛))^𝑡.
(2) The iterative scheme based on the iterative Verlet, see [3], is given as 𝐪_𝑖(𝑡)=𝐪𝑡^𝑛+ℎ𝜕𝐻𝐩𝑡𝜕𝐩^𝑛−ℎ2𝜕𝐻𝐩𝜕𝐪^𝑖−1(𝑡),𝐪_𝑖−1(𝑡),𝐪𝑡^𝑛,𝐩_𝑖(𝑡)=𝐩𝑡^𝑛−ℎ2𝜕𝐻𝐩𝜕𝐪^𝑖−1(𝑡),𝐪_𝑖−1(𝑡)−ℎ2𝜕𝐻𝐩𝑡𝜕𝐪^𝑛−ℎ2𝜕𝐻𝐩𝜕𝐪^𝑖−1(𝑡),𝐪_𝑖−1(𝑡),𝐪𝑡^𝑛+ℎ𝜕𝐻×𝐩𝑡𝜕𝐩_𝑛−ℎ2𝜕𝐻×𝐩𝜕𝐪_𝑖−1(𝑡),𝐪_𝑖−1(𝑡),𝐪𝑡^𝑛,for𝑡∈𝑡^𝑛,𝑡^𝑛+1,ℎ=𝑡^𝑛+1−𝑡^𝑛,𝑛=0,1,…,𝑁,𝑖=1,2,3…,𝐼, (51)(52)(53)(54) where the stopping criterion is 𝐼=3 or 4 or with (62).

Remark 13. For the alternative scheme, we can also improve the initialisation process; see Section 4 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: 𝐪₀,𝐩₀^𝑡=𝐪𝑡^𝑛+1_Euler1st,𝐩𝑡^𝑛+1_Euler1st^𝑡. (55)
(2) We start with a result of the explicit RK method: 𝐪₀,𝐩₀^𝑡=𝐪𝑡^𝑛+1_RK4th,𝐩𝑡^𝑛+1_RK4th^𝑡. (56)
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, 𝐹∶𝐷⊂𝑋→𝑌 with the Banach spaces 𝑋,𝑌 is given with the norms ‖⋅‖_𝑋 and ‖⋅‖_𝑌. Let 𝐹 be at least continuously differentiable; furthermore, we assume that 𝑥₀ is a starting solution of the unknown solution 𝑥^∗.

Then the successive linearisation leads to general Newton's method: 𝐹^′𝑥_𝑖Δ𝑥_𝑖=−𝐹𝑥_𝑖, (57) where Δ𝑥_𝑖=𝑥_𝑖+1−𝑥_𝑖 and 𝑖=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)(59) where 𝟎 is the zero vector.
Further, we have the initial condition, which are given as follows: 𝐩₀,𝐪₀=𝐩𝑡^𝑛,𝐪𝑡^𝑛. (60)
The initial starting solutions for the ith iterative steps are given as follows: 𝐩_𝑖−1(𝑡),𝐪_𝑖−1(𝑡) are the solutions of the (𝑖−1)th iterative step and we have the initial condition for the fixed-point iteration: (𝐩₀(𝑡),𝐪₀(𝑡))^𝑡=(𝐩(𝑡^𝑛),𝐪(𝑡^𝑛))^𝑡.
We apply Newton's method in its discrete form, which is given as 𝐪_𝑖=𝐪_𝑖−1−𝐷𝐅₁𝐩_𝑖−1,𝐪_𝑖−1⁻¹𝐅₁𝐩_𝑖−1,𝐪_𝑖−1,with𝐷𝐅₁𝐩_𝑖−1,𝐪_𝑖−1=𝜕𝜕𝐪_𝑖−1𝐪̇_𝑖−1−𝜕𝐻𝐩𝜕𝐩^𝑖−1,𝐪_𝑖−1,𝐩_𝑖=𝐩_𝑖−1−𝐷𝐅₂𝐩_𝑖−1,𝐪_𝑖−1⁻¹𝐅₂𝐩_𝑖−1,𝐪_𝑖−1,with𝐷𝐅₂𝐩_𝑖−1,𝐪_𝑖−1=𝜕𝜕𝐩_𝑖−1𝐩̇_𝑖−1−𝜕𝐻𝐩𝜕𝐪^𝑖−1,𝐪_𝑖−1,𝑖=1,2,3,…,𝐼, (61) where 𝐷 is the Jacobian matrix (see [13]) and 𝐩_𝑖−1(𝑡),𝐪_𝑖−1(𝑡) are the given solutions of the (𝑖−1)th iterative step.
We assume that we have convergent results after 𝑖=1…,𝐼 iterative steps or with the stopping criterion: max𝐩_𝑖−𝐩_𝑖−1,𝐪_𝑖−𝐪_𝑖−1≤err, (62) while ‖⋅‖ is the Euclidean norm (or a simple vector-norm, e.g., 𝐿₂) and the error bound is given, for example, as err=10⁻⁴.

Example 15. We apply the iterative Verlet scheme to the underlying fixed-point problem (59) of nonlinear Hamiltonian problem (58).
(1) We start with (𝐪₀,𝐩₀)^𝑡=(𝐪(𝑡^𝑛),𝐩(𝑡^𝑛))^𝑡.
(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 matrix 𝐷 and apply 𝐪_𝑖=𝐪_𝑖−1−𝐷𝐅₁𝐩_𝑖−1,𝐪_𝑖−1+𝛿₁𝐪_𝑖−1⁻¹×𝐅₁𝐩_𝑖−1,𝐪_𝑖−1+𝜖₁𝐪_𝑖−1,𝐩_𝑖=𝐩_𝑖−1−𝐷𝐅₂𝐩_𝑖−1,𝐪_𝑖−1+𝛿₂𝐩_𝑖−1⁻¹×𝐅₂𝐩_𝑖−1,𝐪_𝑖−1+𝜖₂𝐩_𝑖−1,𝑖=1,2,3,…,𝐼, (63) where the function 𝛿₁,𝛿₂ can be applied as a scalar, for example 𝛿=10⁻⁶, and the same with 𝜖₁,𝜖₂; 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 𝑧-axis of the system, 𝜔_𝜃 is along the line of the nodes, and 𝜔_𝜓 is along the 𝑧^′-axis. We transform into body coordinates and obtain ̇̇̇̇̇𝜓̇𝜔=𝜙sin𝜃sin𝜓+𝜃cos𝜓𝜙sin𝜃cos𝜓+𝜃sin𝜓𝜙cos𝜃+. (64)

The kinetic energy can be written as 1𝑇=2𝑚𝑥̇²+𝑦̇²+𝑧̇²+𝐴𝜃̇²+𝜙̇²sin²𝜃+𝐶𝜓̇+𝜙̇cos𝜃². (65) The potential energy 𝑈 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 Φ𝑈=𝑚𝑔𝑧−𝜇sin𝜓sin𝜃𝑥Φ+cos𝜓sin𝜃𝑦Φ+cos𝜃𝑧 (66) with 𝑚𝑢 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: 𝑍Ψ=1+𝑍^23/2−𝑋²+𝑌²342𝑍²−3𝑍1+𝑍^27/2. (67) The lengths were scaled by the radius 𝑅 of the base plane, mass is measured in units of 𝑚, and energy in units of 𝑚𝑔ℎ. Therefore the one time unit is 𝑅/𝑔. So the dimensionless Hamiltonian with 𝐪=(𝐗,𝐘,𝐙,𝜃,𝜓,𝜙) is given by 1𝐻=2𝑝²₁+𝑝²₂+𝑝²₃+𝑝²₄𝑎+𝑝₅−𝑃−6cos𝑞₄²𝑎sin²𝑞₄+𝑝₆𝑐−𝑀sin𝑞₄cos𝑞₅𝜕Ψ𝜕𝑞₁sin𝑞₅𝜕Ψ𝜕𝑞₂+cos𝑞₄𝜕Ψ𝜕𝑞₃+𝑞₃ (68) with 𝑎 and 𝑐 being the nondimensionalised inertial parameters and 𝑀 being the ratio of the gravitational and magnetic energy.

For the nonseparable Hamiltonian of (68), we have𝐪̇=𝜕𝐻()=𝑝𝜕𝐩𝐩,𝐪₁,𝑝₂,𝑝₃,𝑝₄𝑎,𝑝₅−𝑝₆cos𝑞₄²𝑎sin²𝑞₄,𝑝₆cos²𝑞₄+(𝑎/𝑐)sin²₄−𝑝₅cos𝑞4𝑎sin²𝑞₄. (69)

The same is given for 𝐩̇=−𝜕𝐻()=𝑀𝜕𝐪𝐩,𝐪sin𝑞₄cos𝑞₅𝜕²Ψ𝜕𝑞²₁+cos𝑞₄𝜕²Ψ𝜕𝑞₁𝜕𝑞₃,𝑀sin𝑞₄cos𝑞₅𝜕²Ψ𝜕𝑞²₂+cos𝑞₄𝜕²Ψ𝜕𝑞₂𝜕𝑞₃,𝑀sin𝑞₄sin𝑞₅𝜕²Ψ𝜕𝑞₂𝜕𝑞₃+cos𝑞₅𝜕²Ψ𝜕𝑞₁𝜕𝑞₃+cos𝑞₄𝜕²Ψ𝜕𝑞²₃𝑀−1,cos𝑞₄sin𝑞₅Ψ𝑞₂+cos𝑞₅𝜕Ψ𝜕𝑞₁−sin𝑞₄𝜕Ψ𝜕𝑞₃−𝑝₆𝑝₅−𝑝₆cos𝑞₄𝑎sin𝑞₄−cos𝑞₄𝑝₅−𝑝₆cos𝑞₄²𝑎sin³𝑞₄,𝑀sin𝑞₄cos𝑞₅𝜕Ψ𝜕𝑞₂−sin𝑞₅𝜕Ψ𝜕𝑞₁,0. (70)

𝐴 and 𝐵 are linearised Lie operators and their vector fields are given by 𝐴=𝜕𝐻⋅𝜕𝜕𝐩𝜕𝐪,𝐵=−𝜕𝐻⋅𝜕𝜕𝐪.𝜕𝐩 (71)

The transfer to the operators is given in the following description.

The exponential operators e^ℎ𝐴 and e^ℎ𝐵 are then just shift operators, with 𝒯₂(ℎ) being a symmetric second order splitting method: 𝒯_2,VV(ℎ)(Δ𝑡)=𝑒^{(Δ𝑡/2)𝐵}𝑒^Δ𝑡𝐴𝑒^{(Δ𝑡/2)𝐵}. (72) This corresponds to the velocity form of the Verlet algorithm (VV).

Furthermore, the splitting scheme 𝒯_2,PV(ℎ)(Δ𝑡)=𝑒^{(Δ𝑡/2)𝐴}𝑒^Δ𝑡𝐵𝑒^{(Δ𝑡/2)𝐴} (73) 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 𝑋, 𝑌, and 𝑍. The trajectory starts at the equilibrium point (𝑞₁,𝑞₂,𝑞₃)=(0,0,1.72).

(a)

(b)

Figure 1 

Trajectory calculated with MPE method supplied with the Verlet algorithm.

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 Figure 2.

(a)

(b)

Figure 2 

Errors of the numerical scheme for MPE method supplied with the iterative Verlet method.

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 Table 1.

Table 1 

Numerical results for MPE method supplied with the iterative Verlet algorithm.

Remark 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⁻⁵ 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 (ℎ=10−5).

Also we tested the 6th order MPE method with different time-steps and compared it with the Runge-Kutta method; see Figure 4.

Figure 4 

Errors of the numerical scheme: 6th order extrapolation scheme with Verlet method as kernel (ℎ=10−5).

Table 2 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.

We could improve the accuracy and computational time needed by using higher order extrapolation schemes; see Table 3.

Table 3 

Errors and computational time with 4th and 6th order MPE schemes using Verlet scheme as kernel.

Remark 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 and ℎ=10⁻³, 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 with 𝑖=1, we were able to remain at order 8 and only 10⁻³ 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.