Improvement of discrete-mechanics-type time-integration schemes by utilizing balance relations in integral form together with Picard-type iterations

of balance, where the time-evolution of the integrands is approximated by established computational techniques of the discrete-mechanics type. In a Picard-type iteration, the outcomes are used for repeating the procedure several times, leading to an increased accuracy. The advantages of the present explicit approach are discussed in the context of linear and nonlinear motions of the mathematical pendulum. We utilize modern sym20


Introduction
The goal of the present contribution is to present a computational strategy for the accurate and efficient numerical time integration of the equations of motion of dynamic mechanical systems, based on the fundamental relations of balance in their most general form, namely the form that involves integrals over time. 30 In the integral form, the relation of balance of some entity ψ that is associated with the mechanical system under consideration can be written as (1.1) Balance relations of momentum, but also of energy and entropy for any mechanical system can be put into a form that is consistent with (1.1). For details on this integral form in the context of the classical field theories of mechanics, see Sect. 194 35 of Truesdell and Toupin [1], who referred to it as the general balance relation. For various examples in the mechanics of solids and fluids, containing the equations of motion of the single particle and the plane pendulum in integral form, see Ziegler [2]. The relation (1.1) can be extended to account for the presence of several nonstandard effects, such as balance of systems with a time-varying mass, see e.g. 40 Irschik and Humer [3]. For instance, in case of balance of linear momentum for a single particle with a constant mass m, there is ψ(t) = mv(t) in (1.1). The absolute velocity vector of the particle at time t is v(t), and ψ(0) = mv(0) denotes the initial linear momentum. In case of balance of moment of momentum for the plane motion of a mathematical 45 pendulum with mass m, a system that we will study subsequently for comparison sake, we have ψ(t) = mL 2 ω(t), where the absolute angular velocity of the pendulum at time t is denoted as ω(t), and L is the distance between the fixed support and the center of gravity of the pendulum mass.
In (1.1), the resultant source that is responsible for the time-wise change of 50 ψ, when written for the intermediate time instant τ , is denoted as R(τ ), where 0 ≤ τ ≤ t. This resultant source represents the instantaneous total of the influx, the supply, and the production quantities of ψ, see [1]. In the special case of balance of linear momentum of a single particle, R stands for the resultant F of the external forces that act upon the particle, while for balance of moment of momentum of 55 a plane pendulum, R represents the resulting moment M of the external forces about the fixed support of the pendulum, see [2]. We consider the integral formulation stated in (1.1) as superior to the wellknown differential form of balance that can be associated with (1.1) by assuming the respective entities to be sufficiently smooth, resulting in an ordinary differential 60 equation of first order. In contrast to the differential form, the integral form (1.1) can directly deal with discontinuous sources R. It also holds when Dirac-type singular functions are present in R. As pointed out in Sec. 194 of [1], relation (1.1) even holds, when there is a jump in ψ, where the integral in (1.1) then has to be understood in the sense of Stieltjes, cf. [1]. Correspondingly, the Stieltjes integral form of (1.1) has been called the fully general form of the relation of balance in [1]. The latter leads to a so-called impulsive balance, a situation, which however also will not be considered in the present contribution.
Important kinematic quantities can be computed from the time-evolution of ψ by further time-wise integration. For instance, the instantaneous position vector 70 r of a single particle follows as (1. 2) The initial position vector is r(0). Analogously, for the plane pendulum, r(t) is to be replaced by the instantaneous inclination angle ϕ(t), measured from the vertical direction, say, where the initial angle is ϕ(0). For the pendulum, one has to replace

75
A main point in the present context of numerical time-initegration is that the integral form of the relation of balance (1.1), as well as the kinematic relation (1.2), already represent time-discrete formulations, which are at our disposal without having introduced any approximation. As a drawback for the usage in numerical time-integration, it however must be noted that the resultant source R(τ ) in (1.1) 80 in general is a function of the yet unknown position r(τ ) and velocity v(τ ) of the system at time τ . For instance, recall that the moment M (τ ) of the weight of the pendulum about the fixed support depends on the angle ϕ(τ ), see [2] and the examples discussed below. Now, since position and velocity at time τ in general are not known in advance, the integrals in (1.1) and (1.2) in general cannot be 85 evaluated directly, and, at a first glance, the latter relations appear to be of little use for the numerical time-integration of mechanical systems.
A strategy from the mathematical theory of ordinary differential equations, which dates back to Emile Picard, and which appears to be suitable for overcoming the above mentioned situation, is the Picard iteration, [4]. In short, substituting the initial condition into the integrand of the integral version of the differential equation of first order at hand as a trial solution, performing the integration and repeating the procedure with the result of this integration as a new trial function, leads to a suitable successive approximation technique. Differential equations of higher orders are solved as a number of simultaneous equations, each of the first 95 order. For some more recent presentations containing concise considerations on existence, convergence and error-bounds, see Chapter 7 of the book by Hirsh, Smale and Devaney [5], Chapter 1.8 of the book by Hairer, Nørsett and Wanner [6], or Chapter 8 of the book by Agarwal and O'Regan [7]. An alternative to the original Picard iteration was presented by Djang [8], who suggested to start the iteration 100 with any convenient function that satisfies the initial conditions, in order to improve convergence of the successive approximation.
It is the goal of the present contribution to specify the suggestion by Djang [8] by utilizing as a first guess well-known numerical schemes that were formulated for the sake of approximately integrating (1.1) or its differential counterparts. Our of the present authors, algorithmic satisfaction of the power law is an important feature; however it does not mean that, in case of conservative mechanical systems, for which a potential of the sources R does exist, the requirement of a conservation of the total of the kinetic and the potential energy is satisfied in an exact sense by the approximate solutions.

130
In the following, we exemplarily demonstrate the validity of our hypothesis, i.e. that improvements can be gained through the repeated usage of the integral balance relation (1.1) in the sense of the Picard iteration, but starting with the established numerical time-integration schemes from the literature as a first guess. We particularly study the motion of the linear and the nonlinear pendulum, for which 135 analytical solutions do exist for the sake of comparison, and which is frequently used in the literature as a benchmark-test. Since we wish to remain as simple as possible, we subsequently use explicit time-integration schemes from the literature for starting the first iteration step.
Our paper is organized as follows: In Sect. 2 below, we introduce some ex-140 plicit time-integration schemes from the literature, which then will be used in order to start the Picard-type iteration procedure. For other explicit time-integration schemes that could be used, the reader is referred to [9] and [10]. Here, we first shortly address two highly valuable explicit time-integration schemes, which were presented by Greenspan in [11], at which occasion he introduced the notion of dis-145 crete mechanics. Subsequently, these two discrete mechanics formulations will be referred to as G1orig and G2orig, respectively. A modification of G2orig using the mid-point rule subsequently will also be studied and will be referred to as G2Morig.
In the course of time, Greenspan and coworkers presented also various implicit time integration schemes, which not only became influential in the field of physics, but 150 which also can be considered as role models for so-called conservative or structure preserving time integrators in nonlinear structural mechanics, see, e. g., [12][13][14][15]. A list of the important contributions of Greenspan can be found in his book on particle modeling [16]. For further implicit time integration techniques, see, e. g., Kim and Reddy [17]. The latter reference will be used subsequently for comparison 155 sake also, likewise to the implicit formulation by Fung [14]. In Sect. 2 below, we also shortly introduce a recent explicit method by Kim and Lee [18], later referred to as KLorig, as well as the explicit fourth-order Runge-Kutta scheme, abbreviated as RK4orig. See, e. g., Chapter 2 of the book by Hairer, Nørsett and Wanner [6] with balance relations" Improvement of discrete-mechanics-type time-integration schemes with balance relations. . . 5 for the explicit foruth-order Runge-Kutta method. With respect to the latter, it 160 is interesting to note that an equivalence of the explicit Runge-Kutta scheme and the original iteration was established by Gur'ianov [19]. Following [19], the solution of the problem by means of the nth order Runge-Kutta method is the same as it follows from the original Picard method after n iterations. An estimation of the error of the Runge-Kutta method was presented in [19] based on this connection.

165
In Sect. 3 below, the explicit time-integration techniques mentioned above exemplarily are applied as a first guess in the Picard-type iteration for computing the motion of the linearized mathematical pendulum, for which an analytic solution in terms of trigonometric functions exists. A single time step with arbitrary initial conditions is studied in order to compare the accuracy of the resulting five 170 approximation schemes, denoted as G1, G2, G2M, KL and RK4, respectively. The errors in the inclination angle and in the total energy, as well as the spectral radii of the provided solutions are presented at the end of the dimensionless time step after four Picard-type iterations. Moreover, a comparison of computational costs of the different approaches is provided. It can be seen that the suggested Picard-type 175 procedures lead to substantial improvements with respect to the accuracy, as well as with respect to the computational costs, since a large number of time steps is needed for the original procedures in order to reach the same accuracy as the novel Picard-type formulations.
Section 4 then delas with the motion of the non-linear mathematical pendu-180 lum, considering a periodic motion with insufficient energy for swinging over, for which an exact solution in terms of Jacoby elliptic functions is at disposal. Periodic motions with a large number of periods are considered using multiple time steps and a comparatively large time-step length. Since the mechanical system under consideration is autonomous, the Picard-type iterations for the five schemes can 185 be performed once, before the time-stepping procedure begins. This is an essential feature with respect to computational efficiency. We note that non-autonomous systems under certain conditions may be rewritten into homogeneous systems, see [5].
For the non-linear pendulum, the integrands in (1) and (2) contain non-linear functions. For a convenient integration, the non-linear terms are expanded into Taylor 190 series, and integration is performed by means of symbolic computation. The number of terms that is taken into account in the Taylor series of course does influence the accuracy at the end of the time step, likewise to the number of Picard-type iterations. An exemplary study with respect to a single time step is presented for the five novel methods G1, G2, G2M, KL and RK4, respectively. Afterwards, a multi-195 step study is performed for the pendulum with large initial angles and over a large number of periods of motion. While the single-step study demonstrates that quite good accuracy can be reached with smaller numbers, particularly for G2M, KL and RK4, the multi-step computations are exemplarily performed using fourteen terms in the Taylor series and four Picard-type iterations, since the computational 200 costs are quite small in comparison to the original versions of the five methods, which need substantially smaller time steps for reaching similar orders of accuracy.
A study on accuracy in comparison to results from the literature is also successfully performed. The paper ends with concluding remarks containing an outlook on future work. 205 2. Some explicit time integration schemes for starting the Picard-type iteration The two explicit time integration schemes introduced by Greenspan in [11] were written for the motion of a single particle, and they are shortly repeated subsequently. The first explicit scheme reads, when applied to the time instant τ : The initial resultant force is written as F (0). Note that (2.4) also follows when substituting (2.1) into (1.2). The scheme outline will be further denoted by G1orig.

215
An extended scheme has been suggested in Greenspan [11] (G2orig is its reference below), in which the mid-point approximation (2.2) is retained, but in which (2.3) is extended. When applied for the time instant τ , this second scheme by Greenspan [11] reads Here, F 1 means the resultant force at time τ /2 corresponding to the arguments for velocity and position, respectively, see formulae (4.1) and (4.2) of [11].
To approximate the position vector at time τ , the estimate originally suggested by Greenspan [11] gives We may use the mid-point approximation rule (2.2) together with (2.5) and (2.6) We may also build our numerical method on the basis of the fourth-order Runge-Kutta method and the recent explicit method of Kim and Lee [18] (the scheme with no dissipation appears to be the most accurate for the considered examples and is used in the paper). (The original versions will be further called RKorig and KLorig.)

230
The Kim-Lee formula leads to where F τ1 and F τ2 mean the resultant force at time τ 1 and τ 2 with the position and velocity given by (2.12) Following our hypothesis described above, the four explicit schemes presented in 235 (2.1)-(2.12) and the fourth-order Runge-Kutta scheme, subsequently abbreviated as RK4orig, now will be applied as starting points for developing novel schemes that are based on the integral formulations stated in (1.1) and (1.2), using a Picard-type iteration. For the formulas necessary to set up RK4orig, see [20]. The displacements and velocity approximations in (2.1)-(2.12) are substituted into the source 240 R(τ ) in (1.1), which here represents the resulting force F (τ ) that acts upon the single particle, and the integral in (1.1) is eventually solved, obtaining a first new estimate for ψ(t) = mv(t). This is facilitated because the estimates in (2.1)-(2.12) enter the specific functional form of F (τ ) via the initial value F (0) only, which is a known entity. In order to proceed further, time t is replaced by τ in (1.1), and 245 the corresponding function v(τ ) is substituted into (1.2) in order to obtain a first new estimate for the displacement r(t) by integration. Alternatively, the mid-point approximation, see (2.2), may be used. Replacing again time t by τ , we obtain r(τ ), which, together with v(τ ), gives a new estimate for R(τ ) = F (τ ) in (1.1). The procedure thus can be repeated a second time, where a second novel estimate is 250 obtained for v(t), and the corresponding v(τ ) can be substituted into (1.2) giving a second novel estimate for r(t). Proceeding further, a single time-step formulation with n iterative steps is obtained, where it is supposed that the nth iteration gives estimates that are more accurate than the (n − 1)th one. Utilizing the results of the nth repetition of r(t) and v(t) as initial values, one can proceed to a second time step, m = 2. Performing analytic integrations, the formulas for v(τ ) and r(τ ) from the nth repetition of the first time step can be used directly, so that no formula updates then are to be performed in the second and in the following time steps, up to the final time step with number m = N . This procedure will be subsequently applied to the motion of the plane mathe-

Linearized pendulum equation
In order to demonstrate advantages of our above suggested strategy, namely to involve and utilize (1.1) and (1.2), we first consider small pendulum oscillations 270 with respect to a stable equilibrium point. The restoring moment is linearized as follows: The exact solution to the linearized problem is known to be, see [2], ϕ(t) = ϕ(0) cos Ωt + ω(0) Ω sin Ωt, Ω = g L .  We use (3.4) and (3.5) to demonstrate exemplarily that this method (and the 285 following novel procedures) leads to an algorithmic satisfaction of the power law. Substituting the Greenspan estimate (3.3) into (3.1), we obtain In order to obtain (3.4), this moment has been substituted into the balance relation (1.1) and the integration has been performed. Now, in the framework of our algorithm, the power of M (τ ) becomes with the novel result ω(τ ) according to (3.4). Setting R(τ ) = P (τ ) in (1.1) and performing the integration, it turns out that with ω(t) from (3.4). This result indeed can be identified as the difference of the instantaneous and initial kinetic energy of the pendulum, see [2], when using the new estimate (3.4) for ω(t). Hence, balance of kinetic energy, or the power law, is 295 satisfied in an algorithmic sense. This is because balance of kinetic energy can be considered as a mathematical consequence of balance of momentum, see [2] and [3] for general derivations. Note that the power law remains satisfied also for the higher-order iterations, since balance of momentum is utilized in each of them. This is also true for the schemes starting with the Greenspan estimates (2.5)-(2.9). 300 We now proceed to the two schemes (G2 and G2M) based on the second scheme in Greenspan [11], see (2.5)-(2.9), to the Kim-Lee scheme (KL) (2.10), and to the fourth-order Runge-Kutta scheme (RK). For the moment in (3.1), we obtain with (2.7) that Hence, using M (0) = −mgLϕ(0), (2.5) becomes The original Greenspan approximation, G2orig, (2.8) then gives The mid-point approximation, G2Morig, (2.9) yields (3.12) (3.15) It is seen that this scheme leads to a further improvement, compare (3.4) and (3.5), which however does not coincide with the next higher terms in the Taylor expansion of the exact solution. The same is true when starting with (3.11) or 315 (3.13).
In Fig. 1 we show the relative error of the five approximate schemes (G1, G2, KL, G2M, and RK) at the end of first time step and for four iterations, n = 4. The error is computed as where T = Ωt is the dimensionless time-step length. The pendulum is assumed to 320 start from ϕ(0) = 1 at rest, ω(0) = 0. Everywhere in this paper we use the symbolic computer code Maple with 24 significant digits. An excellent degree of accuracy can be observed from Fig. 1; e. g., in case of T = π/2, the magnitude of the relative error is in the range from 4.6 × 10 −7 to 6.3 × 10 −9 . Since the mechanical system under consideration, the mathematical pendulum, 325 is conservative, we also present the error in conservation of total energy (kinetic plus potential energy) for these approaches, where we substitute the approximations into the exact formula, see [2], which yields: (3.17) In Fig. 2 the case ω(0) = 0, ϕ(0) = 1 is treated. It turns out that the error is comparatively small, and the method based on the Runge-Kutta scheme (RK) 330 leads to superior results. Note that theses methods differ also by the number of used elementary computations. The comparison is given in Table 1, where RK scheme is most advantageous.  Note: *To reach the accuracy of RK at T = 1.
The last step of investigating the linear problem will be a study on the stability. Following the representation by Hughes [21], we evaluate here the spectral radii 335 ρ of the solutions provided above. The results are shown in Fig. 3. Following [21], we should require that ρ ≤ 1.
It turns out that the stablest method with respect to initial conditions is the G2M, as long as T ≤ 1.903. The RK scheme is stable until 1.854, and the KL scheme is stable until 1.457. We also outline the possible directions of the algorithm optimization to achieve better cost/efficiency ratio: we compare the results for different numbers of Picardtype iterations, see Table 2. Concluding the present linear section, we can state that the suggested procedure leads to substantial improvements with respect to the accuracy of the approxima-345 tions. Starting from discrete mechanics schemes presented in [11], from the wellknown Runge-Kutta scheme, and from the novel Kim-Lee scheme and using four iterations, i. e., four repeated applications of the balance law (1.1), the obtained high degree of accuracy is demonstrated in Fig. 1 for a single time step, m = 1, n = 4, and t = T . All of the schemes do satisfy the power law in an algorithmic 350 sense, but the total energy is not exactly conserved, see Figs. 2 and 3. As can be seen, however, the errors in energy conservation are comparatively small. Conditional stability according to the spectral radius requirement of the G2M, KL, and RK schemes could be proved, while stability is absent for the other schemes.
In connection with the high degree of accuracy, it can be supposed, however, 355 that the solution procedures will work also after a large number of time steps. This will be demonstrated for the nonlinear pendulum in the next section.

Single time step
Consider the plane motion of a mathematical pendulum in the field of gravity, now 360 utilizing the following exact nonlinear expression for the restoring moment, For a periodic motion with insufficient energy for swinging over, the exact solution of the corresponding local relation of balance of angular momentum, is, see e. g. [2] and [22]: ϕ(t) = 2 arcsin (k sn(θ(t)|k)) (4. 3) The notation sn(...) means the Jacobi elliptic sine, and sgn(...) is the signum func-365 tion. The following abbreviations are used: For the sake of numerical approximation, we start the proposed schemes by determining the angular velocity and the angle according to the formulae (2.6) and (2.7): (ω 1 + ω(0)).
The determination of the angle are different for the methods studied here (the conditionally stable method and the most accurate one, see Sec. 3 for modeling of a linear pendulum). For the mid-point approximation method we have Similarly, we may write the formulas for the Kim-Lee (2.10) and fourth-order Runge-Kutta schemes. We now apply (1.1) by analogy in order to obtain an estimate for the angular velocity. Instead of using numerical integration, we expand expressions for sin ϕ in the moment, see (4.1), that result from (4.6), into their Taylor series, in order to simplify the integration of (1.1). The number of expansion terms taken into account Subsequently, we will present results for using 14 terms in the Taylor expansions. The computations can be easily performed symbolically, e. g. by using the Maple procedures taylor and int. The angles at the end of the time step then follow 385 either from the mid-point rule, from the Kim-Lee, or from the fourth-order Runge-Kutta schemes. As in the linear case in Sec. 3, we repeat this procedure according to (1.1) and (1.2) four times, as a basis for the subsequent time-stepping procedure. As above, the schemes are referred to as G1, G2, KL, G2M, and RK.

Multiple time steps 390
In order to compare the outcomes of our suggested explicit methodology with the exact solution presented in [2] and [22], as well as with the outcomes of powerful recent implicit numerical procedures, as presented by Kim & Reddy in [17] (KR), and by Fung in [23] (F), we apply the above derived single step formulae in the framework of a multiple time step procedure. For direct comparison with the nu-395 merical solutions in [17] and [23], we set Ω = 1 s −1 and take the initial velocity as ω(0) = 1.999999238 s −1 , which is used in [17] and [23]. The outcomes of the different computations are shown in Table 3 for the time instant t = T f , which is the quarter of the period of oscillation of the nonlinear pendulum. We note that, in our hands, the exact solution obtained by Maple yields T f = 8.430254948 s, which 400 is slightly different from the value T f = 8.430255141 s mentioned in [17]. However, the dicrepancy in results because of this difference is negligible. Having subdivided T f into various time steps m, in Table 3 we show (i) the exact values as derived by Maple using 24 digits, (ii) the results of the higher order implicit schemes presented in [17] and [23], (iii) the results of the two explicit formulations originally 405 presented by Greenspan [11] (G1orig and G2orig), as well as (iv) the results of the present explicit procedures according to formulae (4.6), Kim-Lee formula (2.10), and fourth-order Runge-Kutta scheme (G2M, KL, and RK). For the meanings of the respective orders of the implicit methods and the parameter µ, the reader is referred to [17] and [23]. It is seen that the formulations in [11] fail to give adequate 410 results, even for a larger number of time steps, such as for m = 100. For m = 25, the present explicit formulations are one order more accurate than the high-order implicit formulations reported in [17] and [23]; for larger time steps, m = 10 and m = 25, the implicit schemes in [17] and [23] are more accurate.
Using the same nonlinear example, we also compare the errors when the different 415 number of Taylor terms is considered. Table 4 presents the results for the Runge-Kutta method with the different number of iterations. It may be useful to compare the schemes where an increasing number of Taylor terms is used for the subsequent iterations. Table 5 shows the corresponding results with the number of Taylor terms different at each iteration. We denote one term 420 in the first iteration step, two in the second one, three in the third one, and four in the fourth one by by 1,2,3,4. Table 3. Comparison of nonlinear numerical solutions at T f = 8.430254948 s for various methods in highly nonlinear case (ω(0) = 1.999999238 s −1 ; ϕmax = 3.139847324). Method Exact 3.139847324 KR (7th-order, µ = 0), [17] 25 3.139833891 0.427840 × 10 −5 KR (8th-order, µ = 1), [17] 25 3.139846872 0.144169 × 10 −6 KR (10th-order, µ = 1), [17] 10 3.139848406 0.344411 × 10 −6 F (7th-order, µ = 0), [23] 25 3.139834046 0.422900 × 10 −5 F (8th-order, µ = 1), [23] 25 3.139847686 −0.115138 × 10 −6 F (10th-order, µ = 1), [  To compare the results of the proposed methods with the exact solution for a longer observation period, the errors in conservation of the total energy are presented first. For simplicity we use a zero initial condition ω(0) = 0. The relative energy error of the proposed solutions can be computed using the following formula, which becomes zero for the exact solution, compare [2]: (4.7) We show this relative error ε in Fig. 4 for the comparatively large number of 50 periods of oscillations, using m = 10000 time steps. In the first case, see Fig. 4a, we treat a large initial angle equal to ϕ(0) = 3.139847324, i. e. nearly vertical 430 initial position. In the second case, see Fig. 4b, we study an initial inclination of ϕ(0) = 3π/4. It can be seen that the error remains small in both cases and for both methods, the RK scheme being superior.
In Fig. 5 we show the phase trajectories of the two considered cases, ϕ(0) = 3.139847324 and ϕ(0) = 3π/4. It turns out that the results of the two proposed 435 methods are practically indistinguishable (they are also practically identical to the exact phase portraits). Most important, the solutions remain at the same phaseplane trajectory during the whole 50 periods of observation. After that we treat the angle over time steps in Figs. 6 and 7. Here we use formula (3.16) to obtain the results for the relative error ε. As shown in Figs. 4 and 440 6, the angle error is more pronounced than the error in energy conservation, but it remains comparatively small, even after 50 periods of nonlinear motion.
We performed the computations for a number of cases and found that in the vast majority the mid-point approximation scheme provides better accuracy with a less computational effort. The methods differ by number of elementary computations, 445 see Table 6, where the mid-point approximation formula is most advantageous. We also have compared our formulation with an explicit time-stepping method  that has been recently suggested for nonlinear problems by Kim  In concluding this section, we can state that the proposed explicit methods yield accurate results also for highly nonlinear pendulum motions and during long observation periods, the computational effort being comparatively low; the errors in energy conservation turn out to be very small; the errors in the angle being 455 also comparatively small, even after 50 periods of nonlinear motion. The solutions remain on the same phase-plane trajectories during the whole observation period and practically do coincide with the exact phase-plane curves.

Conclusions
The main conclusions of the present paper are shortly summarized as follows:

460
• The general relations of balance (the time-discrete integral form of the relations of balance) has been involved in time integration of dynamical systems. • Explicit numerical integration schemes from the literature, such as discrete mechanics schemes, have been taken as starting points for solving the in-465 tegrals, and the procedure has been improved by repeated substitution, resulting in accurate and computationally efficient explicit approximate integration schemes. • Powerful symbolic computer codes can be involved in a straightforward manner.

470
• Linear and nonlinear exact solutions for the plane pendulum motion have been exemplarily used for comparison sake, as well as with results of implicit higher-order formulations and a recent explicit numerical method from the literature. • Improved accuracy also at quite large time steps and during large observa-