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

Xuechuan Wang,Wei He,⋆,Haoyang Feng and Satya N.Atluri

1School of Astronautics,Northwestern Polytechnical University,Xi’an,710072,China

2Department of Mechanical Engineering,Texas Tech University,Lubbock,79409,USA

ABSTRACT Although predictor-corrector methods have been extensively applied, they might not meet the requirements of practical applications and engineering tasks, particularly when high accuracy and efficiency are necessary.A novel class of correctors based on feedback-accelerated Picard iteration (FAPI) is proposed to further enhance computational performance.With optimal feedback terms that do not require inversion of matrices,significantly faster convergence speed and higher numerical accuracy are achieved by these correctors compared with their counterparts; however, the computational complexities are comparably low.These advantages enable nonlinear engineering problems to be solved quickly and accurately, even with rough initial guesses from elementary predictors.The proposed method offers flexibility, enabling the use of the generated correctors for either bulk processing of collocation nodes in a domain or successive corrections of a single node in a finite difference approach.In our method, the functional formulas of FAPI are discretized into numerical forms using the collocation approach.These collocated iteration formulas can directly solve nonlinear problems, but they may require significant computational resources because of the manipulation of high-dimensional matrices.To address this, the collocated iteration formulas are further converted into finite difference forms, enabling the design of lightweight predictor-corrector algorithms for real-time computation.The generality of the proposed method is illustrated by deriving new correctors for three commonly employed finite-difference approaches: the modified Euler approach,the Adams-Bashforth-Moulton approach,and the implicit Runge-Kutta approach.Subsequently,the updated approaches are tested in solving strongly nonlinear problems,including the Matthieu equation,the Duffing equation,and the low-earth-orbit tracking problem.The numerical findings confirm the computational accuracy and efficiency of the derived predictor-corrector algorithms.

KEYWORDS Predictor-corrector method; feedback-accelerated Picard iteration; nonlinear dynamical system; real-time computation

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 pathtracking 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], highdimensional 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

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

It can be rewritten as

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

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

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

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

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

whereJis the Jacobian matrix ofg(x,τ)w.r.tx,andIis 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

which means the error caused will not exceedO2(δ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

whereG(τ)=n(τ)-g(xn,τ).Note that bothJandGdepend onxn.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

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

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

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

Table 1: The two versions of FAPI formula

2.2 Collocated FAPI

2.2.1 The 1st Version of Collocated FAPI

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

We suppose thatx(t)=[x1(t),x2(t),...,xd(t),...,xD(t)]Tis aD-dimensional vector of time-varying variables,and that the variablesxd(t)are approximated by the linear combinations ofNorthogonal basis functionsφd,nb(t),

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

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

where the integration matrixis defined as

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

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

In Eq.(30),g(xn,τ)is denoted asgfor simplicity.As aforementioned,the block diagonal matricesare defined as

Similar approximations can be made to the componentsgd(x,t)ofg(x,t),i.e.,

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

where the integration matrixis defined as

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

Table 2: The two versions of collocated FAPI formula

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

where the superscriptMof matricesdenotes the(D(M-1)+1)-th to(DM)-th row matrix ofthat are related with the collocation pointtM.

Given the values of collocation nodesx(t1),x(t2),...,x(tM-1)being settled,Eq.(38)can be used to update the value of collocation nodex(tM).Similarly, an unknown collocation nodex(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

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

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

The collocation form of Eq.(41)is

It is an implicit finite difference formula.Usually,we letm=Min Eq.(43)to solve forx(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

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

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

wherexn+1andgn,n=0,1,2,...,stand for the(n+1)-th correctedxandn-th correctedg,respectively.Usuallynis 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

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

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

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:

According to Eq.(29),we have

It is rewritten as

wherehis the interval betweent1andt2.

3.2.2 The 2nd Version of ME-FAPI Method

According to Eq.(37),we have

wherehis the interval betweent1andt2.By comparing with Eq.(52),we can see that an acceleration term related with the Jacobian matrix is added at the right hand side.

3.3 Adams-Bashforth-Moulton Method Using FAPI

Taking the 4-th order Adams-Bashforth-Moulton(ABM) method for instance, it solves Eq.(1)with the following formula:

of which Eq.(60) is a predictor using an Adams-Bashforth formula and Eq.(61) is a corrector that uses Picard iteration to solve an Adams-Moulton formula.

3.3.1 The 1st Version of ABM-FAPI Method

To replace the Picard iteration with Feedback-Accelerated Picard iteration in the Adams-Moulton corrector (Eq.(61)), we only need to determine the differentiation matriQ, the integration matriP, the Jacobian matriJ, and matriT.Those matrices simply depend on the collocation pointst2,...,t5and the basis functions of the collocation counterpart of Adams-Moulton formula.Let the basis functions used for interpolation in the Adams-Moulton formula be the Lagrange polynomials.

With Eq.(62),the differentiation matriQis derived as

To keep consistent with the Adams-Moulton formula underlying Eq.(61),we let the integration in Eq.(12)start fromt4.Thus,the matricesPis written as

By substituting Eqs.(63)and(69)into Eq.(29),we have

It is explicitly expressed as

where the terms related withJ(t5)aid to accelerate the convergence of Picard iteration of the Adams-Moulton formula in Eq.(61).

3.3.2 The 2nd Version of ABM-FAPI Method

According to Eq.(37),we have

It is rewritten as

3.4 Implicit Runge-Kutta Method Using FAPI

The equivalence between implicit Runge-Kutta(IRK)method and collocation method has been well discussed and proved[22].For that,the collocated FAPI provides iterative formulae for implicit Runge-Kutta methods.A 4-stage implicit Runge-Kutta method is taken as an example, where its corresponding collocation method takes the first kind of Chebyshev polynomials and Chebyshev-Gauss-Lobatto(CGL)nodes as basis functions and collocation nodes,respectively.

By re-scaling the time segment as [-1,1], the first kind of Chebyshev polynomials are obtained from the recurrence relation

whereξis the rescaled time variable.The CGL nodes areξ1= -1,ξ2= -1/2,ξ3= 1/2,ξ4= 1.The differentiation matriQis obtained as

From Eqs.(23)and(26),the integration matricesandare

The 4-stage implicit Runge-Kutta method using Picard iteration is expressed as

3.4.1 The 1st Version of IRK-FAPI Method

By substituting Eqs.(75)–(77)into the 1st collocated FAPI formula

3.4.2 The 2nd Version of IRK-FAPI Method

The 2nd version of the 4-stage IRK-FAPI formula can be obtained by substituting Eq.(76)into the 2nd collocated FAPI formula as following:

4 Numerical Results and Discussion

In this section, conventional predictor-corrector methods, the corresponding FAPI enhanced methods,and MATLAB-built-in ode45/ode113 were used to solve the dynamical responses of some typical nonlinear systems, such as the Mathieu equation, the forced Duffing equation, and the perturbed two-body problem.The numerical simulation was conducted in MATLAB R2017a using an ASUS laptop with an Intel Core i5-7300HQ CPU.GPU acceleration and parallel processing were not used in the following examples.

In this study, to conduct a comprehensive and rational evaluation of each approach, two approaches were employed to compare the accuracy,convergence speed,and computational efficiency of conventional approaches and the corresponding FAPI-enhanced methods.The first approach involves performing the iterative correction only once, which is the general usage of the predictioncorrection algorithm.The second approach involves repeating the iterative correction until the corresponding terminal conditions are met.

The maximum computational error in the simulation time interval was selected to represent the computational accuracy when evaluating the computational accuracy of each approach.The solution of ode45 or ode113 was employed as the benchmark for calculating this computational error.The relative and absolute errors of ode45 and ode113 were set as 10-15,while the convergence criteria for the predictor-corrector methods are related to specific problems.In each problem,the computational efficiency of different approaches was compared under the same requirement of computational accuracy,with computational time used as the metric for measuring the computational efficiency of the various methods.

Mathieu Equation

The Mathieu equation is

where the parameters are set asδ= 0.5 andε= 0.1.The initial conditions aret= 0,x= 1,and dx/dt= 0.The Mathieu equation is employed to investigate nonlinear vibration problems with periodic forcing and the periodic motion of nonlinear autonomous systems.The stability of the Mathieu equation’s solution depends on the parametersδandε.The set of parameters employed in this study yields stable quasi-periodic motion,which has been confirmed through simulation results obtained from ode45 and other methods.

Forced Duffing Equation

The forced Duffing equation is

where the parameters are set asc= 0.01,k1= 1,k2= 1,f= 7.5,andω= 1.The initial conditions aret= 0,x= 1.5,and dx/dt= 0.The system could exhibit chaotic motion using ode45 and other reliable approaches.Any small state perturbation may result in substantial deviations in the simulation results because of the sensitivity of chaotic systems to the current state.Therefore,when solving chaotic problems,computational errors of general methods with low accuracy accumulate rapidly.

Low-Earth-Orbit Tracking Problem

The orbital problem is considered in the following form:

wherer=[x,y,z]Trepresents the position vector andUrepresents the potential function of the Earth gravity field.The gravity potential is modeled using the 10 deg Earth Gravity Model 2008.The initial condition is set as

The computational accuracy of orbital motion plays a crucial role in determining the success of space missions.In the case of LEO,the long-term trajectory of a spacecraft can be significantly affected by a very small perturbation term of gravity force.A high-order gravitational model should be used to achieve relatively high accuracy.In this case,most computational time is dedicated to assessing the gravitational force terms.The relative error of position is defined asε=‖Δr‖2/‖r‖2,whererrepresents the reference result obtained using ode113 andΔrrepresents the position error at the sampling point.The relative error of the velocity can be modeled similarly.The simulation time for the LEO tracking problem is set to 10 times the orbital period to reliably assess each method.

4.1 Modified Euler Method Using FAPI

Using ode45/ode113 as the benchmark,Figs.2 and 3 show a straightforward comparison of the phase plane portraits, time responses, and computational errors of MEM, 1st ME-FAPI, and 2nd ME-FAPI when solving the Matthieu and forced Duffing equations.Fig.4 shows a comparison of the trajectory of LEO and the computational errors.Table 3 presents the time step size of each method.

Figure 2:Results of Mathieu equation using MEM,ME-FAPI,and ode45:(a)phase plane portrait;(b)time responses;(c)computational error

Table 3: The configurations of MEM and ME-FAPI method

The MEM and ME-FAPI methods exhibited a good match with ode45/ode113.In Fig.2c, the error curve of the 1st ME-FAPI method is relatively steady,whereas the error curves of the 2nd MEFAPI method and MEM gradually increase with time.This indicates that the 1st ME-FAPI method has a significantly lower error accumulation rate than the 2nd ME-FAPI method and MEM when solving the Mathieu equation.In the numerical results of the Duffing equation, an improvement in computational accuracy with the 1st ME-FAPI method is also observed.The error curve of the LEO tracking problem in Fig.4b shows a significant correlation with the current orbital position(or altitude)of the satellite,which is due to the strong association between the macroscopic gravity field and the orbital position (altitude).This observation reveals that the 2nd ME-FAPI method is more accurate in solving the orbit tracking problem.

Figure 3:Results of the forced duffing equation using MEM,ME-FAPI,and ode45:(a)phase plane portrait;(b)time responses;(c)computational error

Figure 4:Results of MEM,ME-FAPI,and ode113 for LEO:(a)trajectory of LEO;(b)computational error

A more comprehensive presentation of the computational accuracy and efficiency of each method is shown in Figs.5–7 when the proposed ME-FAPI method and MEM are only corrected once.The maximum computational error of each method was recorded by sweeping the time step sizeh.This enables us to depict the curve of the computational error to the time step size.In addition,the error tolerance for each method is specified, and the correspondinghvalues are selected to enable a comparison of the computational time for each method under the same computational accuracy.Table 3 shows the duration of the simulated dynamical response.

Figure 5:Performance of MEM and ME-FAPI in solving Mathieu equation when corrected only once:(a)Maximal error;(b)computational time

Figs.5–7 demonstrate that the ME-FAPI method proposed in this study outperforms MEM when corrected once in terms of accuracy and efficiency.The two versions of the ME-FAPI method have varying applicability for different nonlinear systems.For instance,the 1st ME-FAPI method is more suitable for solving the Mathieu equation,where its computational error is 1–2 magnitudes lower than that of the MEM,and its computational time is only 10%of the latter.However,the 2nd ME-FAPI method is more appropriate for the LEO tracking problem, offering 1 order of magnitude higher computational accuracy and a computational speed 1 time faster than the MEM.

Figure 6:Performance of MEM,and ME-FAPI in solving the forced Duffing equation when corrected only once:(a)Maximal error;(b)computational time

Figure 7:Performance of MEM,and ME-FAPI for LEO when corrected only once:(a)Maximal error;(b)computational time

The following section examines the performance indices of each approach after the convergence of the iterative process,including the maximum computational error,average iteration steps,and computational efficiency.Table 3 presents the configurations of these methods, and Figs.8–10 illustrate the simulation results.Since the convergence speed of each method shows similar rules in different cases,only the corresponding simulation results of the Mathieu equation are retained.

Figs.8–10 demonstrate that at least one of the two versions of the ME-FAPI method proposed in this study will have higher computational efficiency than the classical MEM.Only the 1st MEFAPI method is more efficient than MEM when solving the Mathieu equation,as illustrated in Fig.8.However,the 1st and 2nd ME-FAPI methods are more efficient than MEM when solving the Duffing equation.Thus,in practical applications,it is crucial to select the appropriate approach based on the specific problem to be solved.

Figure 8: Performance of MEM, and ME-FAPI in solving Mathieu equation when the algorithms converge:(a)Maximal error;(b)average iteration steps;(c)computational time

Furthermore, Fig.8b illustrates that the 1st ME-FAPI and 2nd ME-FAPI methods converge almost as fast,while the MEM method exhibits the slowest convergence.The advantage of convergence speed substantially enhances the efficiency of the method.For instance,the evaluations of the force model consume most of the computational time when solving the perturbed two-body problem.Therefore,a reduction in the number of iterations plays a significant role in reducing computational cost.

Figure 10:Performance of MEM,and ME-FAPI for LEO when the algorithms converge:(a)Maximal error;(b)computational time

4.2 Adams-Bashforth-Moulton Method Using FAPI

When corrected only once, the simulation results in Figs.11–13 provide an overview of the computational accuracy and efficiency of each method.Table 4 shows the duration of the simulated dynamical response.

Table 4: The configurations of ABM,and ABM-FAPI

Figs.11–13 demonstrate that when corrected only once, the 1st ABM-FAPI method is more accurate and efficient than the ABM method in all cases.For instance, the computational error of the 1st ABM-FAPI method is 1–2 orders of magnitude lower than that of the ABM method at the same step size when solving the Mathieu equation(Fig.11).Furthermore,at the same computational accuracy,the former can save approximately 10%-60%of the computational time of the latter.

Figure 11:Performance of ABM,and ABM-FAPI in solving Mathieu equation when corrected only once:(a)Maximal error;(b)computational time

Figure 12: Performance of ABM, and ABM-FAPI in solving the forced Duffing equation when corrected only once:(a)Maximal error;(b)computational time

New results can be obtained when the algorithms converge,as shown in Figs.14–16.Table 4 shows the configurations of these methods.For the same reason as in the previous section,only the simulation results of the Mathieu equation concerning the convergence speed are retained.

The simulation findings indicate that,in all cases,the 1st ABM-FAPI method significantly reduces computational time compared with the classical ABM method.For instance,as shown in Fig.14c,the 1st ABM-FAPI method can save approximately 28%-65%of the computational time of the ABM method.

Figure 13:Performance of ABM,and ABM-FAPI for LEO when corrected only once:(a)Maximal error;(b)computational time

Figure 14: (Continued)

Figure 15: Performance of ABM, and ABM-FAPI in solving the forced Duffing equation when the algorithms converge:(a)Maximal error;(b)computational time

Fig.14b demonstrates that the convergence speeds of the two versions of the ABM-FAPI method are almost the same.The ABM and ABM-FAPI methods converge at the same speed when the step size is small.However,the convergence speed of the ABM method lags behind that of the ABM-FAPI method as the step size increases.

Furthermore,the algorithm complexity of the 2nd ABM-FAPI method,which incorporates more information about the Jacobi matrix at time nodes,is significantly higher than that of the 1st ABMFAPI method.However,as shown in Figs.11–16,the performance of the 2nd ABM-FAPI method does not significantly outperform that of other methods,whereas the relatively“simpler”1st ABM-FAPI method exhibits a significant advantage in computational accuracy and efficiency.

Figure 16: Performance of ABM, and ABM-FAPI for LEO when the algorithms converge: (a)Maximal error;(b)computational time

4.3 Implicit Runge-Kutta Method Using FAPI

The initial approximation in each step is selected as a straight line without additional computation of the derivatives when using the IRK method and the corresponding IRK-FAPI method to solve the nonlinear system.This initial approximation serves as a“cold start”for the methods.Generally,the“cold start”is a very rough prediction method.

Figs.17–19 show the simulation results when correcting once.Table 5 presents the duration of the simulated dynamical response.Due to the notably poor accuracy of the IRK method, to the extent that it can barely satisfy the corresponding accuracy requirements in some cases, its computational time is not recorded.

Table 5: The configurations of IRK,and IRK-FAPI

Figure 18:Performance of IRK,and IRK-FAPI in solving the forced Duffing equation when corrected only once:(a)Maximal error;(b)computational time

Figure 19:Performance of IRK,and IRK-FAPI for LEO when corrected only once:(a)Maximal error;(b)computational time

Figs.17–19 demonstrate that the IRK-FAPI method outperforms the IRK method with extremely rough initial guesses in terms of accuracy and efficiency.As shown in Fig.18,even with a very small step size, the IRK method cannot predict the long-term dynamical response of the chaotic system,whereas the IRK-FAPI method is very accurate by correcting once.

Similarly,Figs.20–22 show the simulation results when the algorithm converges.Table 5 presents the configurations of these methods.

Figure 20: Performance of IRK, and IRK-FAPI in solving Mathieu equation when the algorithms converge:(a)computational time;(b)average iteration steps

Figure 21: Performance of IRK, and IRK-FAPI in solving the forced Duffing equation when the algorithms converge:(a)computational time;(b)average iteration steps

Figs.20–22 demonstrate that at least one of the two versions of the IRK-FAPI method will have higher computational efficiency than the IRK method, although the difference is not obvious in the case of the Duffing-Holmes’s oscillator.For instance, when solving the Mathieu equation,the 1st ABM-FAPI method can save more than half the computational time of the ABM method.Furthermore,the two versions of the IRK-FAPI method converge at almost the same speed,and both converge faster than the IRK method.

Figure 22:Performance of IRK,and IRK-FAPI for LEO when the algorithms converge:(a)computational time;(b)average iteration steps

5 Conclusion

This study proposes a novel approach for deriving efficient predictor-corrector approaches based on FAPI.Three classical approaches,including the modified Euler approach,the Adams-Bashforth-Moulton approach,and the implicit Runge-Kutta approach,are used to demonstrate how FAPI can be employed to modify and improve commonly used predictor-corrector approaches.For each classical method,two versions of the FAPI correctors are derived.The resulting six FAPI-featured approaches are further used to numerically solve three different types of nonlinear dynamical systems originating from various research fields, ranging from mechanical vibration to astrodynamics.The numerical findings indicate that these enhanced correctors can effectively solve these problems with high accuracy and efficiency,with the 1st version of FAPI featured methods displaying superior performance in most cases.For instance, the 1st FAPI correctors achieve 1–2 orders of magnitude higher accuracy and efficiency than conventional correctors in solving the Mathieu equation.Furthermore, at least one version of the FAPI correctors shows superiority over the original correctors in all these problems,particularly when the correctors are used only once per step.In future studies,the application of the proposed approach to bifurcation problems and high-dimensional nonlinear structural problems will be presented.

The proposed approach offers a general framework for enhancing existing predictor-corrector methods by deriving new correctors using FAPI.Although three existing methods were retrofitted and their improvements demonstrated with three problems,it is anticipated that similar enhancements can be achieved by applying the proposed approach to various other predictor-corrector methods and nonlinear problems.For clarity,methods with low approximation orders(<5)were selected to illustrate our approach.However, methods of higher order can also be derived using the same procedure.In addition,our method can be directly combined with existing adaptive computational approaches for predictor-corrector methods to adaptively solve real-world problems.

Acknowledgement:The authors are grateful for the support by Northwestern Polytechnical University and Texas Tech University.Additionally, we would like to express our appreciation to anonymous reviewers and journal editors for assistance.

Funding Statement:This work is supported by the Fundamental Research Funds for the Central Universities(No.3102019HTQD014)of Northwestern Polytechnical University,Funding of National Key Laboratory of Astronautical Flight Dynamics, and Young Talent Support Project of Shaanxi State.

Author Contributions:The authors confirm contribution to the paper as follows: study conception and design:Xuechuan Wang,Wei He;data collection:Wei He;analysis and interpretation of results:Xuechuan Wang, Wei He, Haoyang Feng; draft manuscript preparation: Xuechuan Wang, Wei He,Haoyang Feng,Satya N.Atluri.All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials:The data supporting the conclusions of this article are included within the article.Any queries regarding these data may be directed to the corresponding author.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.