## A note about high-order semi-implicit differentiation: application to a numerical integration scheme with Taylor-based compensated error

Lo¨ ıc MICHEL 1 and Jean-Pierre BARBOT 1 2 ,

1. Nantes Universit´ - Ecole Centrale de Nantes e Laboratoire des Sciences du Num´ erique de Nantes - UMR 6004 CNRS 1 rue de la No¨ e 44300 Nantes, France (corresponding author : loic.michel54@gmail.com ) 2. QUARTZ EA 7393, ENSEA, Cergy-Pontoise, France

Abstract: In this brief, we discuss the implementation of a third order semiimplicit differentiator as a complement of the recent work by the author that proposes an interconnected semi-implicit Euler double differentiators algorithm through Taylor expansion refinement. The proposed algorithm is dual to the interconnected approach since it offers alternative flexibility to be tuned and to be implemented in real-time processes. In particular, an application to a numerical integration scheme is presented as the Taylor refinement can be of interest to improve the global convergence. Numerical results are presented to support the rightness of the proposed method.

Preprint - Work in progress

## 1 Introduction

Acary et al. [1] introduced an implicit discretization technique from which differentiation techniques, for implicit differential inclusion, have emerged recently (see e.g. a survey in [2]). The principle of this method was initially to replace the classical sign function by an implicit projector, then it has been combined with explicit (homogeneous) terms to define the semi-implicit differentiation [3]; the overall goal is to reduce the effects of chattering. Based on the latest work from the author et al. [4] that proposes a semi-implicit cascade differentiator including Taylor refinement to reduce the estimation error, in this brief, an alternative version of this latter algorithm is proposed that could possibly offer more flexibility towards practical implementation.

The brief is organized as followed. The first Section defines the third order semi-implicit differentiator as a natural extension of the previous works [5] [4] including numerical results, the second Section illustrates an application towards a numerical integration scheme and the third Section gives some concluding remarks and perspectives. Simulation code can be found in the repository:

https://github.com/LoicMichelControl/SemiImplicitDifferentiation.git

## 2 Semi-implicit differentiation algorithm including Taylor refinement

## 2.1 Presentation

Considering a real signal x 1 ( ) t to differentiate, denote x 1 ( kδ ) its corresponding discretization where k is an integer and δ the sampling period. Moreover, x + 1 stands for x 1 (( k + 1) ). δ From the cascaded differentiator [4], we derive a third-order Semi-Implicit Differentiator (SIHD-3), including correction terms (in bold) that refine the Taylor expansion of the two first estimation lines:

$$\begin{cases} \ z _ { 4 } ^ { + } = z _ { 4 } + E _ { 1 } ^ { + } E _ { 2 } ^ { + } E _ { 3 } ^ { + } h \left ( \lambda _ { 4 } \mu ^ { 4 } | e _ { 1 } | ^ { 4 \alpha - 3 } \mathcal { N } _ { 4 } \right ) \\ z _ { 3 } ^ { + } = z _ { 3 } + E _ { 1 } ^ { + } E _ { 2 } ^ { + } h \left ( z _ { 4 } ^ { + } + \lambda _ { 3 } \mu ^ { 3 } | e _ { 1 } | ^ { 3 \alpha - 2 } \mathcal { N } _ { 3 } \right ) \\ z _ { 2 } ^ { + } = z _ { 2 } + E _ { 1 } ^ { + } h \left ( z _ { 3 } ^ { + } - E _ { 3 } ^ { + } h \frac { 1 } { 2 } z _ { 4 } ^ { + } + \lambda _ { 2 } \mu ^ { 2 } | e _ { 1 } | ^ { 2 \alpha - 1 } \mathcal { N } _ { 2 } \right ) \\ z _ { 1 } ^ { + } = z _ { 1 } + h \left ( z _ { 2 } ^ { + } - E _ { 2 } ^ { + } \frac { h } { 2 } z _ { 3 } ^ { + } + E _ { 3 } ^ { + } E _ { 2 } ^ { + } \frac { h ^ { 2 } } { 3! } z _ { 4 } ^ { + } + \lambda _ { 1 } \mu | e _ { 1 } | ^ { \alpha } \mathcal { N } _ { 1 } \right ) \\. \ \, \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \end{cases} \quad.$$

or, equivalently, introducing the SIHD related operator:

$$z _ { i } = \mathfrak { D } _ { f } ^ { ( i ) } ( z, z ^ { + }, x _ { 1 } ), \ \ i = 1 \dots 4$$

where z i is the estimation associated to the i th line of the differentiator (1). Denote y = x 1 the signal to differentiate; ˙ y = x 2 , ¨ = y x 3 and y (4) = x 4 are respectively the first, second and third order differentiation; z , z 2 3 and z 4 are respectively the first, second and third order related estimations; e 1 = x 1 -z 1 is the estimation error.

The projector 1 N q for q = 1 ... 4 is defined by:

$$\mathcal { N } _ { q } ( \epsilon _ { 1 } ) \coloneqq \begin{cases} \ \epsilon _ { 1 } \in S D & \to \mathcal { N } _ { q } = \frac { \lceil \epsilon _ { 1 } \rceil ^ { q ( 1 - \alpha ) } } { \lambda _ { q } ( \mu h ) ^ { q } }, \ E _ { q } ^ { + } = 1 \\ \ \epsilon _ { 1 } \notin S D & \to \mathcal { N } _ { q } = \text{sign} ( \epsilon _ { 1 } ), \ E _ { q } ^ { + } = 0 \end{cases}$$

with the convergence domain 2 SD = { ϵ 1 / ϵ | 1 | ≤ ( λ µh 1 ) 1 q (1 -α ) } .

In eq. (1), the variable E + q allows activating the upper ( q + 1)th line if the convergence domain has been reached at the q th line. This prevent from getting unexpected behavior between lines like oscillations that could destabilize the differentiator.

The corrections terms have been introduced by the author et al. in [5] with respect to the second order differentiation and a study of the convergence has been proposed in [6]. The complete derivation of the refined Taylor expansion in (1) is described in Appendix.

## 2.2 Numerical results

The differentiator (1) is simulated over 10 sec. considering a sine function signal

$$y ( t ) = \sin ( 2 \pi t )$$

under the following initial conditions (that match the derivatives): x 1 (0) = 0 , x 2 (0) = 1 , x 3 (0) = 0 and x 4 (0) = -1. The gains and parameters are set as follows: λ 1 = 10 , 3 λ 2 = 10 , 6 λ 3 = 10 , 9 λ 4 = 10 9 and α = 0 95. .

Figures 1 and 2 illustrate respectively the evolution of the estimated z , z 1 2 , z 3 and z 4 with respect to the corresponding exact differentiation orders of the sine function without extra noise on y .

Figures 3 and 4 illustrate respectively the evolution of the estimated z , z 1 2 , z 3 and z 4 with respect to the corresponding exact differentiation orders of the sine function including extra noise of very small amplitude on y (obtained by a random generator 3 that is multiplied by η 0 = 10 -6 ).

Low estimation errors are observed, thanks to the correction terms, without extra noise. In particular, the error is of order &lt; e 1 &gt; = 2 96 10 . -12 on the z 1 estimation, of order &lt; e 2 &gt; = 4 69 10 . -9 on the estimation of the first derivative z 2 and of order &lt; e 3 &gt; = 4 49 10 . -6 on the estimation of the second derivative z 3 . In presence of noise, the accuracy of the estimation remains correct for the two first orders ( &lt; e 1 &gt; = 7 55 10 . -7 and &lt; e 2 &gt; = 7 47 10 . -4 ) but the noise becomes predominant for the third order, which makes the estimation very noisy. Note that the θ projector [5], which has shown interesting performances with respect to the rejection of the noise, has not been implemented in the current version

1 The projector behaves either as a sign function, or as a nonlinear 'slope' depending on the value of the error e 1 [3].

2 Note that, according to [6], it has been shown that the error e 2 , corresponding to the second estimation line of (1), converges asymptotically to e + 2 = e + 1 /h , hence 'propagating' the errors of estimation through the differentiation orders by 'increment' of h .

3 The noise is computed under Matlab ® using the instruction eta0 * rand() where eta0 is a factor that multiplies the noise.

of the proposed SIHD3 differentiation algorithm. It constitutes an important issue to deal with noisy signals for future investigations.

Figure 1: Example of SIHD3-based differentiation of a sine function without extra noise.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00497v1-figure-fd979dea-19d5-4203-ab2c-a62a4b331a63.jpg)

Figure 2: Example of SIHD3-based differentiation of a sine function without extra noise - estimation errors.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00497v1-figure-a5df81c1-50d3-4759-854a-dc3c23e9121a.jpg)

Figure 3: Example of SIHD3-based differentiation of a sine function including extra noise.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00497v1-figure-7fc27e82-55f2-4dc4-910e-4daed7ee0512.jpg)

Figure 4: Example of SIHD3-based differentiation of a sine function including extra noise - estimation error.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00497v1-figure-6ea587c6-c435-4a19-8309-4270241e965b.jpg)

## 3 Application to numerical integration scheme

## 3.1 Observer-based numerical scheme

Consider an ordinary differential equation (ODE) such as:

$$\frac { d \, y ( t ) } { d \, t } = f ( y ( t ), u ( t ) ), \quad t > t _ { 0 }, \ \ y ( 0 ) = y _ { 0 }.$$

The quantities u and y represent respectively the input and the solution of (4). The corresponding usual discrete explicit Euler scheme reads:

$$\frac { y _ { k + 1 } - y _ { k } } { h } = f ( y _ { k }, u _ { k } ), \quad k \in \mathbb { N }, \ \ y ( 0 ) = y _ { 0 }$$

where k is the sampled time, h is the time-step, and y k , y k +1 are respectively the solution of (5) at the discrete instants t k and t k +1 . The sampled are supposed equally distributed i.e. t k -t k -1 is constant for all k .

In the context of solving ODE, an observer-based numerical integration scheme is proposed to provide a better precision of the estimation of the ODE solution. The proposed scheme is close to the Taylor-based numerical schemes [7, 8] for which we point out the refined Taylor based estimation of f ( y k , u k ) via our proposed SIHD algorithm.

Our proposed differentiation method can be associated to the nonstandard finite-difference (NSFD) methods [9, 10], described as 'powerful numerical methods that preserve significant properties of exact solutions of the corresponding differential equations' (a survey can be found e.g. in [11]) 4 .

The NSFD &amp; SIHD3 based numerical integration scheme reads:

$$\begin{cases} \begin{array} { c } z _ { 1 } ^ { + } = z _ { 1 } + \psi \mathfrak { D } _ { f } ^ { ( 1 ) } ( z, z ^ { + }, \hat { y } ) \end{array}, & k \in \mathbb { N }, y ( 0 ) = y _ { 0 } \end{array} \text{ and } \psi = h + O ( h ^ { 2 } ) \quad ( 6 ) \end{cases}$$

where the SIHD3 differentiator is used to estimate the f part of the ODE and ψ = h + O h ( 2 ) that do satisfy h → 0 according to the NSFD rules . 5

In this work, the purpose is to illustrate the convergence of the proposed NSFD &amp; SIHD3-based algorithm by comparing the evolution of the error with the classical forward Euler scheme and the Runge-Kutta (RK) scheme according to the time. Note that in this scheme, the semi-implicit differentiator (2) appears as a closed-loop observer which tends to adjust and propagate the estimation of f through the numerical integration scheme.

## 3.2 Numerical examples and discussion

Three examples are presented to illustrate some first convergence properties through the evolution of the error | y exact ( ) t -y ODEscheme ( ) t | with the classical forward Euler scheme and the Runge-Kutta (RK) scheme according to the time. y exact ( ) denotes the exact solution of the ODE (4) and t y ODEscheme ( ) denotes t the approximated solution of the ODE (5) using either the classical Euler forward method or the RK method. The internal parameters ( λ, α ) of SIHD3 are

4 Applying a particular differentiator algorithm to the ODE discretization in the context of NFSD was proposed by the author in [12]

5 For example, ψ can be chosen linear of h or exponential of h [13].

the same as in the first Section ( λ 1 = 10 , 3 λ 2 = 10 , 6 λ 3 = 10 , 9 λ 4 = 10 9 and α = 0 95). .

## Example 1 Consider the basic first order equation:

$$\frac { d \, y ( t ) } { d \, t } = - y ( t ) + 1, \ \ y ( 0 ) = 1 0 ^ { - 2 }, \ \ t > 0$$

whose exact solution is y exact = 1 -exp( -t ). Fig. 5 illustrates the comparaison of the errors between the schemes.

## Example 2 Consider the logistic equation [7]:

$$\frac { d y ( t ) } { d t } = 2 y ( 1 - y ), \quad y ( 0 ) = 1 / 5, \quad t > 1 0$$

whose exact solution is y exact = 1 1+4exp(2(10 -t )) . Fig. 6 illustrates the comparaison of the errors between the schemes.

## Example 3 Consider the simple model of a water-discharge problem [14]:

$$\frac { d \, y ( t ) } { d \, t } = - 1 0 y ^ { 2 } + 2 0, \ \ y ( 0 ) = 1 0 ^ { - 2 }, \ \ t > 0$$

whose exact solution is y exact = √ 2 tanh( √ 300 ). t Fig. 7 illustrates the comparaison of the errors between the schemes.

These preliminary numerical experiments show that the proposed NSFDSIHD3 numerical scheme provides a slightly reduced error compared with the

Figure 5: Illustration of the SIHD3 scheme convergence compared with Euler forward scheme and Runge-Kutta scheme - Example 1.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00497v1-figure-fb07227e-929c-4389-8523-0dd8310428ce.jpg)

Euler and Runge-Kutta schemes, especially when the solution is converged towards the asymptotic solution. It is not easy to claim that the transient behavior is more precise than the other scheme since SIHD3 performs a dynamic estimation that depends on the ( λ, α ) parameters, inducing delays during the estimation. Moreover, the NSFD approach appears naturally to adjust the rate of convergence of SIHD3 thanks to the ψ time-step: indeed, ψ can be considered

Figure 6: Illustration of the SIHD3 scheme convergence compared with Euler forward scheme and Runge-Kutta scheme - Example 2.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00497v1-figure-abc741c5-6bc2-4e14-93f8-2c9286b7eaeb.jpg)

Figure 7: Illustration of the SIHD3 scheme convergence compared with Euler forward scheme and Runge-Kutta scheme - Example 3.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00497v1-figure-12ab68ee-de65-4091-81e9-25dbb17cdeb0.jpg)

as an additional SIHD3 'external' gain that can 'virtually' increase the speed of convergence of the scheme, hence giving some flexibility between the SIHD3 parameters ( λ, α ) from the estimation side, and ψ from the numerical scheme side.

## 4 Conclusion and future works

We have derived a differentiator scheme of order three based on semi-implicit techniques that uses Taylor expansion in an implicit way. This differentiator allows reducing the estimation error in comparaison with exact derivatives. Future works include the implementation in real-time in the context of robotics applications. We have also proposed a numerical scheme based on NSFD scheme including the differentiator as a refined estimator of the f term in the right side of the ODE. Further investigations concern a complete study of convergence of the proposed scheme, as well as a comparison of the global error with existing algorithms and possibly an extension toward a semi-implicit scheme.

## Appendix

In this appendix, the corrections terms of the differentiator (1) are calculated from the assumption that the discretized signal x 1 can be differentiated until the third order (inducing a constant fourth order derivative).

Assuming that the discretized signal x 1 can be differentiated though the following 'multi-step' Taylor expansion:

$$x _ { 1 } ^ { + } \ & = \ x _ { 1 } + h x _ { 2 } + \frac { h ^ { 2 } } { 2! } x _ { 3 } + \frac { h ^ { 3 } } { 3! } x _ { 4 } + \frac { h ^ { 4 } } { 4! } x _ { 5 } \\ x _ { 2 } ^ { + } \ & = \ x _ { 2 } + h x _ { 3 } + \frac { h ^ { 2 } } { 2! } x _ { 4 } + \frac { h ^ { 3 } } { 3! } x _ { 5 } \\ x _ { 3 } ^ { + } \ & = \ x _ { 3 } + h x _ { 4 } + \frac { h ^ { 2 } } { 2! } x _ { 5 } \\ x _ { 4 } ^ { + } \ & = \ x _ { 4 } + h x _ { 5 } \\ x _ { 5 } ^ { + } \ & = \ x _ { 5 } \\ \text{ using that the fifth order remains constant } x _ { 5 } = x _ { 5 } ^ { + }, \, the retro-propagation$$

and assuming that the fifth order remains constant x 5 = x + 5 , the retro-propagation of the next-step predictions x + q , q = 2 ... 4 terms allows refining the first Taylor expansions by substituting each x q terms by their such corresponding predictions.

To begin with, the third line:

$$x _ { 3 } ^ { + } = x _ { 3 } + h \, x _ { 4 } + \frac { h ^ { 2 } } { 2! } x _ { 5 }$$

can be rewritten easily and substituting both x 4 = x + 4 -hx 5 and x + 5 = x 5 , allows rewriting:

$$x _ { 3 } ^ { + } \ = \ x _ { 3 } + h \, ( x _ { 4 } ^ { + } - h x _ { 5 } ^ { + } ) + \frac { h ^ { 2 } } { 2! } x _ { 5 } ^ { + }$$

that gives:

$$x _ { 3 } ^ { + } = x _ { 3 } + h \, x _ { 4 } ^ { + } - \frac { h ^ { 2 } } { 2! } x _ { 5 } ^ { + }.$$

In the same way, substituting x 3 = x + 3 -hx + 4 -h 2 2! x 5 and x 4 = x + 4 -hx 5 in the second line gives after few calculations:

$$x _ { 2 } ^ { + } = x _ { 2 } + h \, x _ { 3 } ^ { + } - \frac { h ^ { 2 } } { 2! } x _ { 4 } ^ { + } + R ( x _ { 5 } )$$

The term R x ( 5 ) denotes the 'residual' term associated to the useless fifth order.

Lastly, in the same way, substituting x , x 3 4 and x 2 into x 1 gives:

$$x _ { 1 } ^ { + } = x _ { 1 } + h \, x _ { 2 } ^ { + } - \frac { h ^ { 2 } } { 2! } \, x _ { 3 } ^ { + } + \frac { h ^ { 3 } } { 3! } x _ { 4 } ^ { + } + R ( x _ { 5 } ^ { + } )$$

The bold terms in (13) and (14) correspond indeed to the correction terms that refine the Taylor expansion of the two first lines in (1), taking into account the refinement of the precision of the Taylor terms thanks to the 'retropropagation' of the high order implicit terms.

## Acknowledgments

This work was supported partially by the ANR project DigitSlid ANR-18CE40-0008-01 (https://anr.fr/Project-ANR-18-CE40-0008). Jean-Pierre Barbot is supported with ( R´ egion Pays de la Loire ) Connect Talent GENYDROGENE project.

## References

- [1] V. Acary, B. Brogliato, and Y. Orlov. Chattering-free digital sliding-mode control with state observer and disturbance rejection. IEEE Trans. on Automatic Control , 57(5):1087-1101, 2012.
- [2] Mohammad Rasool Mojallizadeh, Bernard Brogliato, Andrey Polyakov, Subiksha Selvarajan, Lo¨ ıc Michel, Franck Plestan, Malek Ghanes, JeanPierre Barbot, and Yannick Aoustin. A survey on the discrete-time differentiators in closed-loop control systems: Experiments on an electropneumatic system. Control Engineering Practice , 136:105546, 2023.
- [3] Lo¨ ıc Michel, Malek Ghanes, Franck Plestan, Yannick Aoustin, and JeanPierre Barbot. Semi-implicit euler discretization for homogeneous observerbased control: one dimensional case. IFAC-PapersOnLine , 53(2):51355140, 2020. 21st IFAC World Congress.
- [4] Lo¨ ıc Michel, Malek Ghanes, Yannick Aoustin, and Jean-Pierre Barbot. An interconnected discrete time cascaded semi-implicit differentiation. 17th International Workshop on Variable Structure Systems (VSS2024). Accepted , April 2024. https://hal.science/hal-04564290/ .

- [5] Lo¨ ıc Michel, Malek Ghanes, Yannick Aoustin, and Jean-Pierre Barbot. A third order semi-implicit homogeneous differentiator: Experimental results. In 2022 16th International Workshop on Variable Structure Systems (VSS) , pages 77-82, 2022.
- [6] Lo¨ ıc Michel, Malek Ghanes, Franck Plestan, Yannick Aoustin, and JeanPierre Barbot. Semi-implicit homogeneous euler differentiator for a secondorder system: Validation on real data. In 2021 60th IEEE Conference on Decision and Control (CDC) , pages 5911-5917, 2021.
- [7] David Griffiths and Desmond Higham. Numerical Methods for Ordinary Differential Equations: Initial Value Problems . 01 2010.
- [8] Victor G. Ganzha, Ernst W. Mayr, and Evgenii V. Vorozhtsov, editors. Proceedings of the Second Workshop on Computer Algebra in Scientific Computing, CASC 1999, Munich, Germany, May 31 - June 4, 1999 . Springer, 1999.
- [9] Ronald E Mickens. Nonstandard Finite Difference Models of Differential Equations . WORLD SCIENTIFIC, 1993.
- [10] Ronald E Mickens. Advances in the Applications of Nonstandard Finite Difference Schemes . WORLD SCIENTIFIC, 2005.
- [11] Kailash C. Patidar. On the use of nonstandard finite difference methods . † Journal of Difference Equations and Applications , 11(8):735-758, 2005.
- [12] Lo¨ ıc Michel. Algebraic-based nonstandard time-stepping schemes. preprint arXiv 1509.04566 , 2015.
- [13] J. Vigo-Aguiar and Higinio Ramos. A numerical ode solver that preserves the fixed points and their stability. Journal of Computational and Applied Mathematics , 235(7):1856-1867, 2011. Advances in Computational and Mathematical Methods in Science and Engineering.
- [14] C. Vuik, F. Vermolen, M.B. Gijzen, and Mathea Vuik. Numerical Methods for Ordinary Differential Equations . 01 2023.
Loïc Michel Jean-Pierre Barbot
A note about high-order semi-implicit differentiation: application to a numerical integration scheme with Taylor-based compensated error
In this brief, we discuss the implementation of a third order semi-implicit
differentiator as a complement of the recent work by the author that proposes
an interconnected semi-implicit Euler double differentiators algorithm through
Taylor expansion refinement. The proposed algorithm is dual to the
interconnected approach since it offers alternative flexibility to be tuned and
to be implemented in real-time processes. In particular, an application to a
numerical integration scheme is presented as the Taylor refinement can be of
interest to improve the global convergence. Numerical results are presented to
support the rightness of the proposed method.