Fast and Accurate Predictor-Corrector Methods Using Feedback-Accelerated Picard Iteration for Strongly Nonlinear Problems

1  Introduction

In real-world engineering tasks and scientific research, predictor-corrector methods are fundamental to numerically solving nonlinear differential equations [1–3]. Currently, there are hundreds of variants of these methods, which often involve the use of explicit formulas to obtain rough predictions and implicit formulas for further corrections. In general, the predictor-corrector method is widely used in asymptotic and weighted residual methods. This study focuses on finite difference approaches, in which correctors are typically employed once per step. Although multiple corrections can enhance computational accuracy, this approach significantly slows down the algorithm. Therefore, the preference is to reduce the step size rather than repeat corrections. Conventional correctors are often derived from Picard iteration (PI) to solve nonlinear algebraic equations (NAEs) associated with implicit formulas [4,5]. Despite their ease of implementation, these correctors may be insufficient in generating accurate results and ensuring the reliability of long-term simulations for strongly nonlinear dynamical systems. In practical applications, it can result in destructive time delays in tasks that aim for highly accurate results, or it can fail to achieve highly accurate results in real-time computation. Such instances are found in large-scale dynamics simulations [6], autonomous navigation [7], and path-tracking control [8].

This study presents three main contributions. First, an approach to derive fast and accurate correctors to enhance the performance of predictor-corrector methods for more advanced tasks in real life is proposed. This is achieved by replacing PI with feedback-accelerated Picard iteration (FAPI) [9]. These two iteration approaches are categorized as asymptotic approaches that iteratively solve analytical solutions of nonlinear systems starting from an initial guess [10]. The authors in [11] have demonstrated that PI can be considered a special case of FAPI. The latter approach is based on the concept of correcting an approximated solution through the integration of an optimally weighted residual function. FAPI degenerates into PI if the optimal weighting function is specified as a unit matrix instead of being derived from the optimal condition [12]. Therefore, FAPI converges faster and is more stable. Second, since FAPI cannot be employed directly to solve NAEs, the iterative formula is further discretized with appropriate basis functions and collocation nodes [13] to obtain numerical correctors corresponding to PI and FAPI. The PI correctors obtained in this method are equivalent to conventional correctors that directly employ PI to solve implicit finite difference formulas. Third, the proposed approach can be directly applied to other existing predictor-corrector approaches to enhance performance.

The applications of the enhanced predictor-corrector methods and FAPI span various domains, including simulations of n-body problems in celestial mechanics [6], molecular dynamics [14], high-dimensional nonlinear structural dynamics [15], and other problems in fluid mechanics, solid mechanics, heat transfer, and related areas. Particularly, these methods can offer straightforward applications in aerospace engineering for the accurate on-board prediction of spacecraft flight trajectories to achieve autonomous precision navigation [7].

This study is structured as follows. In Section 2, the relationship between PI and FAPI is revealed, and the methodology is introduced by constructing numerical correctors from the discretized collocation formulas of FAPI, which are derived for a general nonlinear dynamical system. Section 3 briefly reviews the underlying relationship between finite difference and collocation, followed by the derivation of FAPI correctors corresponding to three classical finite difference approaches: the modified Euler approach [16], the Adams-Bashforth-Moulton approach [17], and the implicit Runge-Kutta approach [18]. In Section 4, three nonlinear problems (the Matthieu equation [19], the Duffing equation [20], and the low-earth-orbit (LEO) tracking problem [21]) are used as numerical examples to verify the effectiveness of the proposed FAPI correctors.

2  Methodology

2.1 From PI to FAPI

Consider a nonlinear dynamical system in the form of

dxdt=g(x,t),t∈[t0,tf](1)

Without loss of generality, we suppose x(t0)=0. Picard iteration method solves this system via

xn+1(t)=∫t0tg(xn,τ)dτ(2)

It can be rewritten as

xn+1(t)=xn(t)−∫t0t{x˙n(τ)−g(xn,τ)}dτ(3)

where x˙n(τ)−g(xn,τ) is the residual error of Eq. (1). To make Picard iteration converge faster, we can add a weighting function in the integral of residual error. Then we obtain

xn+1(t)=xn(t)+∫t0tλ(τ){x˙n(τ)−g(xn,τ)}dτ(4)

Suppose Π[x(τ),λ(τ)] is a vector function of x(τ) and λ(τ),

∏[x(t),λ(t)]=x(τ)|τ=t+∫t0tλ(τ){x.(τ)−g(xn,τ)}dτ(5)

Let x^(τ) be the exact solution of Eq. (1). Naturally, it satisfies the express

Π[x^(τ),λ(t)]=x^(τ)|τ=t+∫t0tλ(τ){x^.(τ)−g(x^,τ)}dτ=x^(τ)|τ=t(6)

Now we want to make all the components of the function Π[x(t),λ(t)] stationary about x at x(t)=x^(t). First, the variation of Π[x(t),λ(t)] is derived as

δΠ[x(t),λ(t)]=δx(τ)|τ=t+δ∫t0tλ(τ){x.(τ)−g(x,τ)}dτ                                  =δx(τ)|τ=t+∫t0tδλ(τ){x.(τ)−g(x,τ)}dτ+∫t0tλ(τ)δ{x.(τ)−g(x,τ)}dτ                                  =∫t0tδλ(τ){x.(τ)−g(x,τ)}dτ+δx(τ)|τ=t+λ(τ)δx(τ)|τ=t0τ=t                                     −∫t0t[∂λ(τ)∂τ+λ(τ)∂g(x,τ)∂x]δx(τ)dτ−∫t0tλ(τ)∂g(x,τ)∂τδτdτ(7)

Then we collect the terms including δx(τ)|τ=t and δx(τ),

δx(τ)|τ=t+λ(τ)δx(τ)|τ=t,    ∫t0t[∂λ(τ)∂τ+λ(τ)∂g(x,τ)∂x]δx(τ)dτ(8)

Note that the boundary value of x(τ) at τ=t0 is prescribed, that is to say δx(τ)|τ=t0=0. Thus, the stationary condition for Π[x(τ),λ(τ)] is obtained as

{I+λ(τ)|τ=t=0∂λ∂τ+λ(τ)J(τ)=0(9)

from which we can derive the first order appoximation of λ.

λ≈−I+J(t)(τ−t)(10)

where J is the Jacobian matrix of g(x,τ) w.r.t x, and I is an identity matrix. In Picard iteration, the stationarity conditions cannot be satisfied and thus δΠ[x(t),λ(t)]=O(δx). However, in FAPI, since the stationarity conditions can be satisfied, then

δx(τ)|τ=t+λ(τ)δx(τ)|τ=t0τ=t−∫t0t[∂λ(τ)∂τ+λ(τ)∂g(x,τ)∂x]δx(τ)dτ=0(11)

which means the error caused will not exceed O2(δx). Therefore, we can conclude that the newly proposed method based on FAPI is more accurate than the existing methods based on PI, and to reach the given tolerance, the former needs fewer iterations.

By substituting Eq. (10) into Eq. (4), the correctional formula becomes

xn+1(t)=xn(t)−{I+J(t)t}∫t0tG(τ)dτ+J(t)∫t0tτG(τ)dτ(12)

where G(τ)=x˙n(τ)−g(xn,τ). Note that both J and G depend on xn. Since Eq. (12) involves feedback terms of integration of residual errors that accelerate the convergence of Picard iteration, it is named as Feedback-Accelerated Picard Iteration (FAPI).

Previously, the formula of FAPI takes the form of Eq. (12). In this work, we propose another form of it. By differentiating Eq. (4), we have

dxn+1(t)dt=dxn(t)dt+λ(τ)|τ=t{x.n(τ)|τ=t−g(xn,τ)|τ=t}                            +∫t0t∂λ(τ)∂t{x.n(τ)−g(xn,τ)}dτ(13)

with Eqs. (9), (13) can be rewritten as

dxn+1(t)dt=g(xn,t)+J(t)∫t0tλ(τ){x˙n(τ)−g(xn,τ)}dτ(14)

Simply approximating λ(τ) as −I, which is its zeroth order approximation, Eq. (14) becomes

dxn+1(t)dt=g(xn,t)−J(t)[xn(t)−x(t0)−∫t0tg(xn,τ)dτ](15)

By integrating both sides of Eq. (15), we have

xn+1(t)=x(t0)+∫t0t{g(xn,τ)−J(τ)[xn(τ)−x(t0)−∫t0τg(xn,ξ)dξ]}dτ(16)

For clarity, the two versions of FAPI fomulae are listed in Table 1.

2.2 Collocated FAPI

2.2.1 The 1st Version of Collocated FAPI

By collocating M nodes in time domain, Eq. (12) can be discretized. It is written as

[xn+1(t1)xn+1(t2)⋮xn+1(tM)]=[xn(t1)xn(t2)⋮xn(tM)]− (I~+J~T~)[∫t0t1G(τ)dτ∫t0t2G(τ)dτ⋮∫t0tMG(τ)dτ]+J~[∫t0t1τG(τ)dτ∫t0t2τG(τ)dτ⋮∫t0tMτG(τ)dτ](17)

The block diagonal matrices I~, J~, T~ are defined as

I˜=[II⋱I],J˜=[J(t1)J(t2)⋱J(tM)],T˜=[t1It2I⋱tMI](18)

We suppose that x(t)=[x1(t),x2(t),...,xd(t),…,xD(t)]T is a D-dimensional vector of time-varying variables, and that the variables xd(t) are approximated by the linear combinations of N orthogonal basis functions ϕd,nb(t),

xd(t)=∑nb=1Nαd,nbϕnb(t)=Φ(t)Ad(19)

Using the approximation in Eq. (19), it is obtained that

[G(t1)G(t2)⋮G(tM)]=[x˙n(t1)x˙n(t2)⋮x˙n(tM)]−[g(t1)g(t2)⋮g(tM)]=R~Q~R~−1[xn(t1)xn(t2)⋮xn(tM)]−[g(t1)g(t2)⋮g(tM)](20)

where R~ is a row-rearranging-matrix. The differentiation matrix Q~ is defined as

Q˜=[QQ⋱Q],Q=[Φ.(t1)Φ.(t2)⋮Φ.(tM)][Φ(t1)Φ(t2)⋮Φ(tM)]−1(21)

Further, by making the basis functions in Eq. (19) universal, it is found that

[∫t0t1G(τ)dτ∫t0t2G(τ)dτ⋮∫t0tMG(τ)dτ]=R~P~R~−1[G(t1)G(t2)⋮G(tM)](22)

where the integration matrix P~ is defined as

P˜=[PP⋱P],P=[∫t0t1Φ(τ)dτ∫t0t2Φ(τ)dτ⋮∫t0tMΦ(τ)dτ][Φ(t1)Φ(t2)⋮Φ(tM)]−1(23)

For simplicity, we denote Q~=R~Q~R~−1 and P~=R~P~R~−1. In our previous work, we simply let

[∫t0t1τG(τ)dτ∫t0t2τG(τ)dτ⋮∫t0tMτG(τ)dτ]=R~P~R~−1[t1G(t1)t2G(t2)⋮tMG(tM)](24)

However, it may cause deterioration of the numerical solution, because τG(τ) is approximated with the same basis functions of G(τ). It means that the approximation of τG(τ) is one order lower than it should be. A more rigorous approximation is made as following:

[∫t0t1τG(τ)dτ∫t0t2τG(τ)dτ⋮∫t0tMτG(τ)dτ]=R~P~τR~−1[G(t1)G(t2)⋮G(tM)](25)

where the modified integration matrix P~τ is derived as

P˜τ=[PτPτ⋱Pτ],Pτ=[∫t0t1τΦ(τ)dτ∫t0t2τΦ(τ)dτ⋮∫t0tMτΦ(τ)dτ][Φ(t1)Φ(t2)⋮Φ(tM)]−1(26)

Denoting P~τ=R~P~τR~, and substituting Eqs. (20), (22), and (25) into Eq. (17), it can be found that

[xn+1(t1)xn+1(t2)⋮xn+1(tM)]=[xn(t1)xn(t2)⋮xn(tM)]+{J~(P~τ−T~P~)−P~}{Q~[xn(t1)xn(t2)⋮xn(tM)]−[g(t1)g(t2)⋮g(tM)]}(27)

We further denote the constant matrix P~τ−T~P~ as H~, thus Eq. (27) becomes

x~n+1=x~n+(J~H~−P~)(Q~x~n−g~n)(28)

where x~=[x(t1),x(t2),...,x(tM)]T, g~=[g(t1),g(t2),...,g(tM)]T. Denote J~H~−P~ as J~HP, the algorithm is further simplified as

[xn+1(t1)xn+1(t2)⋮xn+1(tM)]=[xn(t1)xn(t2)⋮xn(tM)]+J~HP{Q~[xn(t1)xn(t2)⋮xn(tM)]−[g(t1)g(t2)⋮g(tM)]}(29)

2.2.2 The 2nd Version of Collocated FAPI

Further numerical discretization is made to the collocation form of Eq. (15). It is expressed as

[x˙n+1(t1)x˙n+1(t2)⋮x˙n+1(tM)]=[g(xn,t1)g(xn,t2)⋮g(xn,tM)]−J~{[xn(t1)xn(t2)⋮xn(tM)]−[x(t0)x(t0)⋮x(t0)]−[∫t0t1gdτ∫t0t2gdτ⋮∫t0tMgdτ]}(30)

In Eq. (30), g(xn,τ) is denoted as g for simplicity. As aforementioned, the block diagonal matrices I~, J~, T~ are defined as

I˜=[II⋱I],J˜=[J(t1)J(t2)⋱J(tM)],T˜=[t1It2I⋱tMI](31)

We suppose that x˙(t)=[x˙1(t),x˙2(t),...,x˙d(t),…,x˙D(t)]T is a D-dimensional vector of time-varying variables. The components x˙d(t) are approximated by the linear combinations of N orthogonal basis functions ϕd,nb(t),

x˙d(t)=∑nb=1Nαd,nbϕnb(t)=Φ(t)Ad(32)

Similar approximations can be made to the components gd(x,t) of g(x,t), i.e.,

gd(x,t)=∑nb=1Nβd,nbϕnb(t)=Φ(t)Bd(33)

Using the approximations in Eqs. (32), (33), it is obtained that

[x(t1)x(t2)⋮x(tM)]=R~P~R~−1[x˙(t1)x˙(t2)⋮x˙(tM)]+[x(t0)x(t0)⋮x(t0)](34)

and

[∫t0t1g(τ)dτ∫t0t2g(τ)dτ⋮∫t0tMg(τ)dτ]=R~P~R~−1[g(t1)g(t2)⋮g(tM)](35)

where the integration matrix P~ is defined as

P˜=[PP⋱P],P=[∫t0t1Φ(τ)dτ∫t0t2Φ(τ)dτ⋮∫t0tMΦ(τ)dτ][Φ(t1)Φ(t2)⋮Φ(tM)]−1(36)

For simplicity, we denote P~=R~P~R~−1. By substituting Eqs. (30), (35) into Eq. (34), it can be found that

x~n+1=x~(t0)+P~g~n⏞Picard iteration−P~J~{x~n−x~(t0)−P~g~n}⏟Acceleration terms(37)

where x~=[x(t1),x(t2),...,x(tM)]T, x~(t0)=[x(t0),x(t0),...,x(t0]T, and g~=[g(t1),g(t2),...,g(tM)]T. It will be shown later that the Jacobian matrix J~ can be regarded as constant in many cases, without harming the convergence of algorithm.

For clarity, the two versions of collocated FAPI formulae are displayed in Table 2.

2.3 FAPI as Numerical Corrector

Each row of Eqs. (29) or (37) can be used independently as a numerical corrector. For instance, the last row of Eq. (37) can be written as

xn+1(tM)=x(t0)+P~M[g(t1)g(t2)⋮g(tM)]−P~MJ~{[xn(t1)−x(t0)xn(t2)−x(t0)⋮xn(tM)−x(t0)]−P~[g(t1)g(t2)⋮g(tM)]}(38)

where the superscript M of matrices P~M denotes the (D(M−1)+1)-th to (DM)-th row matrix of P~ that are related with the collocation point tM.

Given the values of collocation nodes x(t1), x(t2),…, x(tM−1) being settled, Eq. (38) can be used to update the value of collocation node x(tM). Similarly, an unknown collocation node x(tm) can be updated using the (D(m−1)+1)-th to (Dm)-th row of Eq. (37) once the values of the other collocation nodes are settled. The methodology of FAPI corrector is illustrated in Fig. 1.

Figure 1: The process to derive FAPI corrector for predictor-corrector method when dealing general nonlinear ODEs

3  Finite Difference Integrators Featured by FAPI

To apply FAPI as numerical corrector in finite difference method, we need to use the same collocation nodes and basis functions underlying the finite difference discretization. The relationship between finite difference and collocation is briefly reviewed herein.

3.1 From Collocation to Finite Difference

The collocation form of Eq. (1) can be written as

Q~[x(t1)x(t2)⋮x(tM)]=[g(t1)g(t2)⋮g(tM)](39)

where Q~ is the rearranged differentiation matrix. The (D(m−1)+1)-th to (Dm)-th row of Eq. (39) is

Q~m[x(t1)x(t2)⋮x(tM)]=g(tm)(40)

Eq. (40) is a finite difference form of Eq. (1). If m≠M, it is an explicit finite difference formula, since x(tM) can be explicitly obtained using Eq. (40) once the other values of collocation points are known. If m=M, it becomes an implicit formula.

Besides, Eq. (1) can be equivalently expressed as

x(t)=∫t0tg(x,τ)dτ+x(t0)(41)

The collocation form of Eq. (41) is

[x(t1)x(t2)⋮x(tM)]=P~[g(t1)g(t2)⋮g(tM)]+[x(t0)x(t0)⋮x(t0)](42)

where P~ is the rearranged integration matrix. The (D(m−1)+1)-th to (Dm)-th row of Eq. (42) is

x(tm)=P~m[g(t1)g(t2)⋮g(tM)]+x(t0)(43)

It is an implicit finite difference formula. Usually, we let m=M in Eq. (43) to solve for x(tM).

Since Eqs. (43) and (40) are obtained by numerical differentiation and integration, respectively, using collocation method, we can find the counterpart of a finite difference method as a collocation method by choosing proper basis functions and collocation nodes, or derive a finite difference formula from collocation method.

3.2 Modified Euler Method Using FAPI

Herein, we consider the Modified Euler method (MEM) that only uses 2 nodes, of which one is known and the other is unknown. The collocation counterpart of this method has 2 collocation nodes, and takes the Lagrange polynomials ϕ0(t)=(t−t2)/(t1−t2) and ϕ1(t)=(t−t1)/(t2−t1) (or simply ϕ0(t)=1 and ϕ1(t)=t)) as basis functions. From Eq. (21), we have

Q=[Φ˙(t1)Φ˙(t2)][Φ(t1)Φ(t2)]−1=1t2−t1[−11−11](44)

and

Q~=R~Q~R~−1=1t2−t1[−II−II](45)

By substituting Eq. (45) into Eq. (39), we have

1t2−t1[−II−II][x(t1)x(t2)]=[g(t1)g(t2)](46)

The first row of Eq. (46) is an explicit formula of Euler method, while the second row is implicit. They can be used independently or conjointly in a predictor-corrector manner. In the later case, the explicit formula provides prediction for the implicit formula to make further correction. Usually, the correction is made by the Picard iteration method due to its low computational complexity and inversion-free property. Based on Eq. (46), a predictor-corrector formula using Picard iteration can thus be expressed as

1t2−t1[−II][x(t1)x0(t2)]T=g(t1)(47)

1t2−t1[−II][x(t1)xn+1(t2)]T=gn(t2)(48)

where xn+1 and gn, n=0,1,2,..., stand for the (n+1)-th corrected x and n-th corrected g, respectively. Usually n is simply set as 0, implying the correctional formula is only executed for once.

However, the modified Euler method utilizes the trapezoidal formula instead of the implicit Euler formula to make correction. The trapezoidal formula can be derived from Eq. (42). By using the same collocation nodes and basis functions that lead to Eq. (46), we have

P=[∫t1t1Φ(τ)dτ∫t1t2Φ(τ)dτ][Φ(t1)Φ(t2)]−1=t2−t12[0011](49)

and

P~=R~P~R~−1=t2−t12[00II](50)

Note that we start the integration at t1 for simplicity. By substituting Eq. (50) into Eq. (42), we have

[x(t1)x(t2)]=t2−t12[00II][g(t1)g(t2)]+[x(t1)x(t1)](51)

The second row of Eq. (51) is the formula of trapezoidal method. Given x(t1), x(t2) can be solved with Picard iteration, which is

xn+1(t2)=t2−t12[II][g(t1)gn(t2)]T+x(t1)(52)

Eqs. (47) and (52) constitute the modified Euler method.

3.2.1 The 1st Version of ME-FAPI Method

Now we replace Eq. (52) with the feedback-accelerated Picard iteration. The necessary matrices are derived first as following:

J˜=[J(t1)J(t2)],T˜=[t1It2I],P˜τ=t2−t16[00(2t1+t2)I(2t2+t1)I](53)

With Eqs. (50) and (53), H~ is derived as

H~=P~τ−T~P~=−(t2−t1)26[002II](54)

Thus

J~HP=J~H~−P~=−(t2−t1)26[002J(t2)J(t2)]−t2−t12[00II](55)