AAlgebraic-based nonstandard time-stepping schemes
Lo¨ıc MICHEL
Abstract
In this preliminary work , we present nonstandard time-stepping strate-gies to solve differential equations based on the algebraic estimation methodapplied to the estimation of time-derivative, which provides interestingproperties of ”internal” filtering. We consider firstly a classical finite dif-ference method, like the explicit Euler method for which we study thepossibility of using the algebraic estimation of derivatives instead of theusual finite difference to compute the numerical derivation. Then, we in-vestigate how to use the algebraic estimation of derivatives in order toimprove the slope predictions in RK-based schemes. Contents This work is distributed under CC license http://creativecommons.org/licenses/by-nc-sa/4.0/ . Email of the corresponding author: [email protected] a r X i v : . [ m a t h . NA ] S e p Introduction
Ordinary differential equations and stiff differential equations [1, 2] have been studiedextensively and successful methods have been proposed (e.g. [3, 4, 5, 6, 7, 8, 9, 10]),including the solver routines ’ode45’ and ’ode23s’ [11, 12, 13]. Recent alternativemethods propose for instance to use specific series, like power series [14] or Haarwavelets [15, 16] to describe the solution of ODEs; the purpose is to substitute thestandard finite difference by the numerical properties of the series. In this paper, wepropose to extend the finite difference technique by an estimation of the derivativesusing the algebraic estimation approach. Introduced in [17], the algebraic estimationmethod [18, 19] has been widely applied to many different problems of estimation thatoccur in dynamical systems. Some of these applications aim e.g. to reconstruct thestates of dynamical systems [20, 18], and a computational toolbox has been releasedto help processing efficiently the algebraic estimation of derivatives for particularproblems [21]. The algebraic estimation of derivatives could be seen as similar to thedifferentiation by integration technique, investigated especially by [22] and [23], andprobabilistic Kriging-based method [24].Our proposed method belongs a priori to the class of nonstandard finite-difference(NSFD) methods [25, 26, 27, 28, 29], described as ”powerful numerical methods thatpreserve significant properties of exact solutions of the corresponding differential equa-tions” (an interesting survey can be found in [30]); explicit rules to ”design” NSFDschemes have been proposed in [25] [26]. Among the derivations that have been pro-posed (e.g. [31, 32, 33, 34]), one emphasizes the contribution of [34] that describesa nonstandard finite-difference scheme for fractional systems, that uses a discreteversion of the Caputo fractional derivation.In this paper, we attempt to design simple NSFD schemes, that are derived fromthe algebraic estimation technique in order to evaluate the discrete derivatives of firstorder ODEs considering an ”Euler framework” and a ”RK framework”. In an ”Euler-like scheme”, the algebraic estimation of the derivatives is mainly used to build amulti-step scheme, where the evaluation of the derivative depends also on the pastvalues of the solution. In a RK-like scheme, the estimation of the local slope, thatis usually ”measured” as an average of several local slopes, is performed using theproperties of filtering of the algebraic technique.The paper is structured as follow. Section 2 described the proposed methods,including a brief review of the algebraic estimation method. Some concluding remarkscan be found in Section 3. 1
Outline of the method
Consider an ordinary differential equation (ODE) such as: d y ( t ) d t = f ( y ( t ) , u ( t )) , t ∈ [0 , t f ] , y (0) = y . (1)The quantities u and y represent respectively the input and the solution of (1). Thecorresponding usual discrete explicit Euler scheme reads: d y ( t ) d t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k +1 = y k +1 − y k h ≈ f ( y k , u k ) , k ∈ N , y (0) = y (2)where k is the sampled time, h is the time-step, and y k , y k +1 are respectively thesolution of (2) at the discrete instants t k and t k +1 . The sampled are supposed equallydistributed i.e. t k − t k − is constant for all k .We identified two major strategies to build time-stepping schemes, based on therules that help designing NSFD schemes and the algebraic estimation technique. The”Euler-like”(or ”multi-step”) strategy follows the principle of the Euler method, andthe ”RK-like” strategy is based on the general Runge-Kutta method [12]. Algebraic estimator of derivatives
Using the previous notations, consider a function g ( x ) and define the algebraic esti-mator of derivative D ( g ) n of order n such as: d g ( x ) d x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k +1 ≈ D ( g ) n = α g k +1 + α g k + · · · + α n g k − n φ ( h ) ,φ ( h ) = h + O ( h ) as h → h = x k − x k − is constant (3)The coefficients α , α , · · · , α n are real coefficients. These coefficients and the function φ ( h ) are defined by (17) in the proposition (2.1) as calculating steps of the proposedtime-stepping schemes . We have: α = Tα i = 2( T − hi ) , i ∈ { · · · η − } α n = ( T − hη )1 φ ( h ) = 3 KhT = 3 Kh ( ηh ) (4) The algebraic estimator of derivative is defined with the highest degree of g that is equal to k + 1 . Depending on the application, this degree can be decreased to k . .1 Euler-like Algebraic-NSFD scheme The proposed ”Euler-like” Algebraic-NSFD scheme aims to extend the finite differencein (1) by the algebraic estimator D ( f ) n . First, one proposes a ”symbolic” nonstandardscheme, that is of the form: D ( f ) n ≈ f ( y k , u k )i . e . α y k +1 + α y k + · · · + α n y k − n +1 φ ( h ) ≈ f ( y k , u k ) , k ∈ N ∗ + , y = y (0) . (5)A forward scheme can be easily deduced from the proposed general scheme (5): y k +1 = − α y k + · · · + α n y k − n +1 α + φ ( h ) α f ( y k , u k ) , k ∈ N ∗ + , y = y (0) . (6)As presented in [34] regarding the Caputo fractional derivative, due to the nonlocalnature of the algebraic-based derivative operator, the discrete representation of thederivative must take into account a part of the past history of the solution. Thenumber n of the involved sampled solutions y k , · · · , y k − n defines a window that char-acterizes the ”precision” of the derivative estimation .Then, in the following proposition, we formalize the proposed Euler-like nonstan-dard time-stepping scheme based on the algebraic estimation framework. Proposition 2.1.
Consider the following nonstandard numerical scheme associatedto the ODE (1), that verifies the general scheme (5): K hT σ + η − X j =1 y k − η + j +1 ( T − jh ) ≈ f ( y k , u k ) , k ∈ N ∗ + , y = y (0) with T = ηh > η ∈ N ∗ + and σ = y k − η +1 T + y k +1 ( T − ηh ) (7) where k is the sampled time, h is the time-step, K is a real constant, and T is amultiple of h that characterizes the ”low filtering” property of the algebraic derivative(see § T . This estimation window implies that the n initial sampled solutions are obviously not known atthe beginning of the algorithm. To initialize the estimation window, one may consider e.g. usingthe classical finite difference scheme. roof. Hypothesis We consider solving the ODE (1), for which the solution y ( t ) isassumed, in the time domain, to be locally represented by a linear function of thetime i.e.: y ( t ) = a + a t (8)where, in particular, the coefficient a is calculated from the algebraic estimationtechnique. The lowest degree of ”time-filtering” is considered.Step 1 - Algebraic Derivation The general technique to perform the estimation ofderivatives using the algebraic estimation strategy is mainly described in [17, 18].It allows estimating the coefficient a in (8). We consider only the first algebraicderivative of f according to the definition of the initial-value problem (1). Transformfirst (8) in the Laplace domain, then take the derivative d/ds : y ( s ) = a s + a s ⇐⇒ sy ( s ) = a + a s (9) y ( s ) + s d yd s = − a s (10)Step 2 - Back to the time domain Multiply first (10) by s − p : s − p y ( s ) + s − p +1 d yd s = − s − p − a (11)then using the Cauchy formula applied to each term of (11): Z ( β ) s α d β x ( τ ) d τ β d β τ = 1( α − Z t ( t − τ ) α − ( − β τ β x ( τ ) dτ (12)we deduce the expression of a that corresponds to the derivative of y through (8): a = 1( p − Z tt − T ( t − τ ) p − y ( τ ) d τ − p − Z tt − T ( t − τ ) p − τ x ( τ ) d τ p + 1)! Z tt − T ( t − τ ) p +1 d τ (13)Since we consider the lowest degree of ”time-filtering”, we take p = 2 . Some elements of proof can be found in [35] (pp.22-23). The expression of a corresponds to the direct application of (12), whose integral is evaluatedover a moving windows T . a = − T Z tt − T ( t − τ ) y ( τ ) d τ (14)Step 3 - Integral expression of the numerical scheme Replacing in (1) the time-derivative operator by (14) gives an integral expression of the proposed E-A-NSFDscheme: − T Z tt − T ( t − τ ) y ( τ ) d τ ≈ f ( y ( t ) , u ( t )) (15)Step 4 - Finally, the discretized version... The usual Trapezoidal scheme allows in-tegrating numerically and thus providing a discretized expression of (15). We have: a = − K hT y k − η +1 T + y k +1 ( T − ηh ) + η − X j =1 y k − η + j +1 ( T − jh ) (16)which, as a result, gives the proposed NSFD scheme (7) when substituted in (2) ac-cording to the rule [30]. Dependingon the value of T , the total gain of | a | could be different, that is why we multiply a by a positive constant K that should be estimated at the initial instants . Note thatif η = 1 , then T = h which cancel the summation term P η − j =1 y k − η + j +1 ( T − jh ) .Step 5 - In addition, identification with the general scheme The coefficients α i , i ∈{ · · · n } can be identified to the general form (5). We have: α = Tα i = 2( T − hi ) , i ∈ { · · · η − } α n = ( T − hη )1 φ ( h ) = 3 KhT = 3 Kh ( ηh ) (17)As an example of the Prop. 2.1, the simplest version of the proposed E-A-NSFDis a particular case where T = h that verifies (2). Since the complete finite difference scheme in (2) is a priori replaced by a weighted sum of thediscrete values y k , we deduce that the rule ” h ” term has to be changed). A future work should investigatethis ambiguity. In simulation, it is observed that K → when T increases.
5o illustrate how the E-A-NSFD scheme is written, consider the simple ODE: d y ( t ) d t = 5 y ( t ) + u ( t ) , t ∈ [0 , t f ] , y (0) = y (18) • Considering η = 3; h = 10 − ( Ψ is a real constant), the E-A-NSFD scheme of(18) reads: Ψ (3 u k − + 2 y k − − y k − y k +1 ) ≈ y ( t ) + 50 u ( t ) (19) • Considering η = 5; h = 10 − ( Ψ is a real constant), the E-A-NSFD scheme of(18) reads: Ψ (5 y k − + 6 y k − + 2 y k − − y k − − y k − y k +1 ) ≈ y ( t ) + 50 u ( t ) (20)As stated in Prop. 2.2, the distribution of the signs in these two cases is equallydistributed. This strategy could been seen as a receding horizon approach because the interval ofintegration T corresponds to the window from which the estimation of the derivativeis performed. Therefore, the constant T should be chosen small enough in order toevaluate the estimate within an acceptable short delay, but large enough, in order toensure a good low pass filtering [21].Since the whole set of coefficients α i does not depend of the discrete solution y k ,they need only to be evaluated at the initialization step of the scheme.The following simple proposition gives an interesting property of the generalscheme (5) concerning the distribution of the signs of the set of α i . Proposition 2.2.
The distribution of the signs (over the whole set of the coefficients α i , i = 0 · · · η ) verifies: sign( α , · · · , α b η/ − c ) = − sign( α b η/ c , · · · , α η ) .Proof. From Prop. 2.1, we have α i = γ ( T − hi ) , i = 0 · · · η, γ ∈ N ∗ + . To get thenumber of α i that are e.g. positive, check the inequality ( T − hi ) > . This gives i < b η c .The E-A-NSFD scheme may become unstable for high η . The window T of the algebraic estimation may be seen as an ”internal” time-delay inside the ODE(2) that may ”deform” the value of the time-derivative estimation. Further works may investigatethe conditions of stability according to the length of the window T and the coefficient K . .1.3 Toward possible generalizations One may consider a possible generalization of the scheme (7) using a more generalexpression φ ( h ) instead of the constant time-step h . The scheme can be rewritten: KφT σ + η − X j =1 y k − η + j +1 ( T − jφ ) ≈ f ( y ( t ) , u ( t )) (21)with φ ( h ) = h + O ( h ) , as h → . The family of RK methods is described generally as an Euler scheme for which mul-tiple estimations of the local slopes around f ( y k ) are computed in such manner thatthe solution increment y k +1 is based upon a weighted average of these multiple esti-mations. The general RK scheme reads: y k +1 = y k + h s X i =1 b i q i (22)where q , q , · · · , q s are the estimations of the different slopes through f . The number s characterizes the choice of the method in the RK family. The proposed ”RK-like” Algebraic-NSFD scheme aims to substitute the weighted sum P si =1 b i q i in (22) by the algebraic estimator D ( f ) n .First, one proposes a ”symbolic” nonstandard scheme, that is of the form: y k +1 = y k + φ ( h ) f ( y k + φ ( h )[ f ( y k ) + φ ( h ) D ( f ) n ]) i . e .y k +1 − y k φ ( h ) = f ( y k + φ ( h )[ f ( y k ) + φ ( h )( α f ( y k + δ ) + · · · + α n f ( y k + δ n ))]) k ∈ N ∗ + , y = y (0) (23)where the coefficients δ , δ , · · · , δ n are real coefficients such as for all i ∈ [0 , n ] , δ i ∈ [ δ min , δ max ] ; α · · · α n are real coefficients and φ ( h ) = h + O ( h ) , as h → .Then, in the following proposition, we formalize the proposed RK-like nonstandardtime-stepping scheme based on the algebraic estimation framework.7 roposition 2.3. Consider the following nonstandard numerical scheme associatedto the ODE (1), that verifies the general scheme (23): y k +1 − y k φ ( h ) = f y k + h f ( y k ) + 3 Kh T σ + η − X j =1 f ( y k + δ j )( T − jh ) , with T = ηh > η ∈ N ∗ + and σ = f ( y k + δ ) T + f ( y k + δ η )( T − ηh ) ,k ∈ N ∗ + , y = y (0) (24) where k is the sampled time, h is the time-step, K is a real constant, and T is amultiple of h that characterizes the ”low filtering” property of the algebraic derivative(see § T .Proof. Hypothesis We consider solving the ODE (1), for which the function f isassumed, to be locally represented by a linear function i.e.: f ( x ) = a + a x (25)where, in particular, the coefficient a is calculated from the algebraic estimationtechnique. The lowest degree of ”time-filtering” is considered. The purpose is toevaluate the local value of a for x ∈ [ x k − , x k ] , for all k ≥ .Step 1 - From the discretized version of the derivative estimation From Prop. 2.1,the expression of a in the discrete domain is given by: a = − K hT f ( x k − η ) T + f ( x k )( T − ηh ) + η − X j =1 f ( x k − η + j )( T − jh ) (26)The set of values x k − η , x k − η +1 , · · · , x k is, according to (3), theoretically a regular grid,but since the interval [ x k − , x k ] , k ≥ is very small, we can consider evaluating η + 1 values of f ( y k + δ j ) where for all i ∈ [0 , η + 1] , δ i ∈ [ δ min , δ max ] . In particular, the δ i coefficients can be random numbers. We have: a = − K hT f ( x k + δ ) T + f ( x k + δ η )( T − ηh ) + η − X j =1 f ( x k + δ j )( T − jh ) . (27)8tep 2 - Definition of an RK2 scheme with algebraic estimation of f To build theproposed NSFD scheme, we start from the standard RK2 scheme.Considering the evolution of the solution y between the steps y k and y k +1 , for any k ≥ . The first step of the RK2 scheme, called the prediction , is to evaluate thesolution y k +1 / in the middle of the step i.e.: y k +1 / = y k + h f ( y k ) (28)that gives, via (2), the derivative of y at the middle of the time-step: d yd t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k +1 / = f ( y k +1 / ) = f y k + h f ( y k ) ! . (29)Then, the second step, called the correction , is to evaluate the solution at the end ofthe step y k taking into account the derivative d yd t (cid:12)(cid:12)(cid:12) k +1 / i.e.: y k +1 = y k + hf ( y k +1 / ) . (30)Instead of considering a single evaluation of the derivative at the middle of the time-step (or multiple evaluations of the derivatives like e.g. in the RK4 case), the proposedmodification of the RK2 strategy is to estimate multiple local values of the slopes via f between the steps y k and y k +1 using the algebraic estimator D ( f ) n . First, onegeneralizes the ”prediction” step (28) that is rewritten: y k +1 / = y k + h < f ( y k ) > h (31)where we denote < f ( y k ) > h the averaged value of f ( y ) between the steps k and k + 1 (in other words, over the time-step h , regarding the notations) .Before rewriting the ”correction” step, one needs the formal definition of the deriva-tive of a function g ( x ) taken at the point x = a : d gd x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = a = lim h → g ( a + h ) − g ( a ) h ≈ g ( a + h ) − g ( a ) h for h very small (32)that allows defining the quantity < f ( y k ) > h , by substituting (27) in (32): < f ( y k ) > h ≈ f ( y k ) + d f ( y ) d y (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k h = f ( y k ) + ha (33)Finally, from (31) and (33), the correction step is thus rewritten : Regarding the RK4 algorithm (with no time dependence), we have: < f ( y k ) > h = ( f ( y k ) +2( f ( y k + k ) + 2 f ( y k + k ) + f ( y k + k )) where k , k and k are the estimated local slopes. A simple approximation in (30) can be considered: since we assumed that f is locally describedas a first order polynomial (25), it could be possible to approximate f ( y k ) by γy k where γ is a realconstant number that may be estimated at each time-step. k +1 = y k + hf ( y k +1 / ) = y k + hf y k + h < f ( y k ) > h ! (34) ⇐⇒ y k +1 = y k + hf y k + h f ( y k ) + ha ) ! . (35)As a result, (35) gives the proposed NSFD scheme (24) according to the rule T , the total gain of | a | could be different, thatis why we multiply a by a positive constant K that should be estimated at the initialinstants.Since f ( y k +1 ) − f ( y k ) represents the slope that locally ”drives” the discrete differen-tial equation (2), we are looking for a good estimation of the predicted slope through f at a point y , that is located between f ( y k ) and f ( y k +1 ) . Considering smooth enoughthe function f ( x ) where x ∈ [ y k , y k +1 ] , the purpose of the operator D ( f ) n is to esti-mate the derivative of f in this particular interval in such manner that an ”averaged”value of f ( x ) for x ∈ [ y k , y k +1 ] can be deduced. Remark 2.4.
This scheme (23) is the dual form of the Euler-like scheme (5), wherethe algebraic estimator is located on the left side of (1) and is directly ”connected” tothe discrete values of y . In this RK-like scheme, the algebraic estimator is located onthe right side of (1) and is connected to the discrete values of y through f . Remark 2.5.
The lowest degree of ”time-filtering” has been considered in the proof ofthe Prop. 2.1 2.3. More investigations would clarify the impact of higher degrees onthe precision of the solution and the stability of the proposed schemes. Moreover, thiscould be an essential factor regarding the possibility of simulating dynamical modelswith noisy signals (see Rem. 2.6).
The main difference in the utilization of the algebraic operator between the two pro-posed schemes is the ”management” of the sampled y solution. In the Euler-likescheme, the window T is a moving window that aims to performs the derivative es-timation while the scheme runs; the window is thus initialized only at the beginningof the scheme. In the RK-like scheme, the window T is a static window, that aimsto perform the derivative estimation between two steps [ y k , y k − ] ; the window is thusinitialized at each time-step and the number of evaluations of f is proportional to thelength of the window T (especially, to compute 27).The Euler-like scheme may require a priori few samples of y in order to compute D ( f ) n while preserving a priori the global stability of the scheme. At the opposite,10he RK-like scheme requests many samples due to the fact that a random process isinvolved in the evaluation of f / estimation of the slopes. In such case, the average ofthe estimated slopes is performed naturally by the filtering property of the algebraicestimator.As presented in § h by a function φ that follows the rule Remark 2.6.
Taking into account the properties of filtering that provide the algebraicestimation framework, we assume that it could be possible to simulate dynamical mod-els that involve noisy signals, like e.g. models in electronics that sometimes can notbe dissociated from the physical noise [36].
We described in this paper, nonstandard finite-difference schemes that use the alge-braic estimation framework in order to compute an estimate of the derivatives. Wederived two algebraic-based nonstandard schemes. The approach used for the firstscheme is similar to [34], which derives a nonstandard scheme from the Caputo deriva-tive definition in order to solve fractional dynamical systems, and the second schemeaims to extend the Runge-Kutta method. Further work include: • extensive tests of practical problems (e.g. using the set of practical problemspresented in [37]); • a complete stability study in order to characterize the stability domains accord-ing to the stiffness of the ODE; • the application to higher order ODEs and utilization of the dedicated toolbox[21] to systematize the proposed A-NSFD procedure; • the application to fractional differential equations. Acknowledgement
The author is sincerely grateful to Dr. Edouard Thomas for his strong guidance andhis valuable comments that improved this paper.11 eferences [1] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff andDifferential-Algebraic Problems, Springer-Verlag, 1991.[2] J. C. Butcher, Numerical Methods for Ordinary Differential Equations, WileyBook, 2003.[3] N. Guzel and M. Bayram, On the numerical solution of stiff systems, Appl. Math.Comp., 170(1):230-236, 2005.[4] S. Abelman and K.C. Patidar, Comparison of some recent numerical methods forinitial-value problems for stiff ordinary differential equations, Comput. Math. withApplications, 55(4):733-744, 2008.[5] W.H. Enright, T.E. Hull and B. Lindberg, Comparing numerical methods for stiffsystems of ODEs, BIT, 15(1):10-48, 1975.[6] A. Jannelli and R. Fazio, Adaptive stiff solvers at low accuracy and complexity,Journal of Computational and Applied Mathematics, 191(2):246-258, 2006.[7] G. Hojjati, M.Y. Rahimi Ardabili and S.M. Hosseini, A-EBDF: an adaptivemethod for numerical solution of stiff systems of ODEs, Math. Computers in Sim-ulation, 66(1):33-41, 2004.[8] M.T. Darvishia, F. Khanib and A.A. Soliman, The numerical simulation for stiffsystems of ordinary differential equations, Computers & Mathematics with Appli-cations, 54(7-8):1055-1063, 2007.[9] J. R. Cash, On the integration of stiff systems of O.D.E.s using extended backwarddifferentiation formulae, Numerische Mathematik, 34(3):235-246, 1980.[10] D. I. Ketcheson and U. bin Waheed, A comparison of high order explicitRunge-Kutta, extrapolation, and deferred correction methods in serial and par-allel, preprint arXiv, Mar. 2014.[11] G. D. Byrne and A. C. Hindmarsh, Stiff ODE solvers : A review of current andcoming attractions, Journal of Computational Physics archive, 70(1):1-62, 1987.[12] L. F. Shampine and M. W. Reichelt, The MATLAB ODE Suite, SIAM J. Sci.Comput., 18(1):1-22. 1997.[13] R. Ashino, M. Nagase and R. Vaillancourt, Behind and beyond the Matlab ODEsuite, Computers & Mathematics with Applications, 40(4-5):491-512, 2000.1214] B. Benhammouda et al. , Procedure for Exact Solutions of Stiff Ordinary Difer-ential Equations Systems, British Journal of Mathematics & Computer Science4(23):3252-3270, 2014.[15] C.H. Hsiao, Haar wavelet approach to linear stiff systems, Mathematics andComputers in Simulation, 64(5):561-567, 2004.[16] ¨U. Lepik, Haar wavelet method for solving stiff differential equations, Mathe-matical Modelling and Analysis, 14(4):467-481, 2009.[17] M. Fliess and H. Sira-Ramirez, Closed-loop parametric identification forcontinuous-time linear systems via new algebraic techniques, in eds. H. Garnier andL. Wang, Identifcation of Continuous-time Models from Sampled Data , Springer,pp. 362-391, 2008.[18] H. Sira-Ram´ırez, C. Garc´ıa Rodr´ıguez, J. C. Romero and A. L. Ju´arez, AlgebraicIdentification and Estimation Methods in Feedback Control Systems, Wiley, 2014.[19] M. Fliess, C. Join, and H. Sira-Ramirez, Non-linear estimation is easy, Int. Jour-nal of Modelling, Identification and Control, vol. 3, 2008.[20] J. Reger and J. Jouffroy, On Algebraic Time-Derivative Estimation and DeadbeatState Reconstruction. 48th IEEE Conference on Decision and Control, Shanghai,China, 2009[21] J. Zehetner and J. Reger and M. Horn, A Derivative Estimation Toolbox basedon Algebraic Methods - Theory and Practice, IEEE International Conference onControl Applications, pp.331,336, 2007.[22] X. Huang, C. Wu and J. Zhou, Numerical differentiation by integration, Journal:Math. Comp., 83(286):789-807. 2014.[23] Z. Wang and R. Wen, Numerical differentiation for high orders by an integra-tion method, Journal of Computational and Applied Mathematics, 234(3):941-948,2010.[24] E. Vazquez and E. Walter, Estimating derivatives and integrals with Kriging,44th IEEE Conference on European Control Conference Decision and Control(CDC-ECC ’05), pp.8156-8161, 2005.[25] R.E. Mickens, ”Nonstandard Finite Difference Models of Differential Equations”,World Scientific (book), 2000.[26] R.E. Mickens, ”Advances in the Applications of Nonstandard Finite DifferenceSchemes”, World Scientific (book), 2005.1327] P. M. Manning and G. F. Margrave, Introduction to non-standard finite-difference modelling, CREWES Research Report - Volume 18 (2006).[28] R. Anguelov and J.M.-S. Lubuma, Contributions to the mathematics of the non-standard finite difference method and applications, Num. Methods Partial Differ-ential Equations, 17(5):518-543, 2001.[29] S. Abelmana and K. C. Patidar, Comparison of some recent numerical methodsfor initial-value problems for stiff ordinary differential equations, Computers & Mathematics with Applications, 55(4):733-744, 2008.[30] K. C. Patidar, On the use of nonstandard finite difference methods, Journal ofDifference Equations and Applications, 11(8):735-758, 2005.[31] U. Erdogana and T. Ozis, A smart nonstandard finite difference scheme for sec-ond order nonlinear boundary value problems, Journal of Computational Physics,230(1):6464-6474, 2011.[32] L.-I. W. Roeger, ”General nonstandard finite-difference schemes for differentialequations with three fixed-points”, Computers & Mathematics with Applications,57(3):379-383, 2009.[33] R. Anguelov and J. M.-S. Lubuma, Nonstandard finite difference method bynonlocal approximation, Mathematics and Computers in Simulation, 61(3-6):465-475, 2003.[34] M. Y. Ongun, D. Arslan, and R. Garrappa, Nonstandard finite difference schemesfor a fractional-order Brusselator system, Advances in Difference Equations, vol.2013, article 102, 2013.[35] R. Delpoux, Contribution `a l’identification, l’estimation et la commande de Mo-teurs Synchrones `a Aimants Permanents (MSAP), Ph.D. dissertation thesis, EcoleCentrale de Lille, France, 2012.[36] T. Sickenberger and R. Winkler, Efficient transient noise analysis in circuit sim-ulation, PAMM, 6(1):55-58, 2006.[37] Test set for IVP solvers