Qinghua Zhang
email: qinghua.zhang@inria.fr

Adaptive Kalman Filter for Actuator Fault Diagnosis

Keywords: FDI, adaptive observer, joint state-parameter estimation, LTV/LPV systems

An adaptive Kalman filter is proposed in this paper for actuator fault diagnosis in discrete time stochastic time varying systems. By modeling actuator faults as parameter changes, fault diagnosis is performed through joint state-parameter estimation in the considered stochastic framework. Under the classical uniform complete observability-controllability conditions and a persistent excitation condition, the exponential stability of the proposed adaptive Kalman filter is rigorously analyzed. The minimum variance property of the combined state and parameter estimation errors is also demonstrated. Numerical examples are presented to illustrate the performance of the proposed algorithm.

INTRODUCTION

In order to improve the performance and the reliability of industrial systems, and to satisfy safety and environmental requirements, researches and developments in the field of fault detection and isolation (FDI) have been continuously progressing during the last decades [START_REF] Hwang | A survey of fault detection, isolation, and reconfiguration methods[END_REF]. Model-based FDI have been mostly studied for linear time invariant (LTI) systems [START_REF] Gertler | Fault Detection and Diagnosis in Engineering Systems[END_REF][START_REF] Chen | Robust Model-Based Fault Diagnosis for Dynamic Systems[END_REF][START_REF] Isermann | Fault Diagnosis Systems: An Introduction From Fault Detection To Fault Tolerance[END_REF][START_REF] Ding | Model-based Fault Diagnosis Techniques: Design Schemes, Algorithms, and Tools[END_REF], whereas nonlinear systems have been studied to a lesser extent and limited to some particular classes of systems [START_REF] De Persis | A geometric approach to nonlinear fault detection and isolation[END_REF][START_REF] Xu | Nonlinear system fault diagnosis based on adaptive estimation[END_REF][START_REF] Berdjag | Nonlinear model decomposition for robust fault detection and isolation using algebraic tools[END_REF]. This paper is focused on actuator fault diagnosis for linear time-varying (LTV) systems, including the particular case of linear parameter varying (LPV) systems. The problem of fault diagnosis for a large class of nonlinear systems can be addressed through LTV/LPV reformulation and approximations (Lopes dos [START_REF] Lopes Dos Santos | An LPV modeling and identification approach to leakage detection in high pressure natural gas transportation networks[END_REF][START_REF] Tóth | The behavioral approach to linear parametervarying systems[END_REF]. It is thus an important advance in FDI by moving from LTI to LTV/LPV systems. In this paper, actuator faults are modeled as parameter changes, and their diagnosis is achieved through joint estimation of states and parameters of the considered LTV/LPV systems. Usually the problem of joint stateparameter estimation is solved by recursive algorithms known as adaptive observers, which are most often studied in deterministic frameworks for continuous time systems [START_REF] Marino | Nonlinear control design. Information and system sciences[END_REF][START_REF] Zhang | Adaptive observer for multiple-inputmultiple-output (MIMO) linear time varying systems[END_REF][START_REF] Besançon | On adaptive observers for state affine systems[END_REF][START_REF] Farza | Adaptive observers for a class of uniformly observable systems with nonlinear parametrization and sampled outputs[END_REF]. Discrete time systems have been considered in [START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF][START_REF] Besançon | Adaptive observer design for discrete time LTV systems[END_REF], also in deterministic frameworks. In order to take into account random uncertainties with a numerically efficient algorithm, this paper considers stochastic systems in discrete time, with an adaptive Kalman filter, which is structurally inspired by adaptive observers [START_REF] Zhang | Adaptive observer for multiple-inputmultiple-output (MIMO) linear time varying systems[END_REF][START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF], but with well-established stochastic properties.

The main contribution of this paper is an adaptive Kalman filter for discrete time LTV/LPV system joint state-parameter estimation in a stochastic framework, with rigorously proved stability and minimum variance properties. Different adaptive Kalman filters have been studied in the literature for state estimation based on inaccurate statespace models. Most of these algorithms address the problem of unknown (or partly known) state noise covariance matrix or output noise covariance matrix [START_REF] Mehra | On the identification of variances and adaptive Kalman filtering[END_REF][START_REF] Brown | Adaptive Kalman filtering[END_REF], whereas the case of incorrect state dynamics model is treated as incorrect state covariance matrix. In contrast, in the present paper, the new adaptive Kalman filter is designed for actuator fault diagnosis, by jointly estimating states and parameter changes caused by actuator faults.

PROBLEM STATEMENT

The discrete time LTV system subject to actuator faults considered in this paper is generally in the form of1

x(k) = A(k)x(k-1) + B(k)u(k) + Φ(k)θ + w(k) (1a) y(k) = C(k)x(k) + v(k), (1b)
where k = 0, 1, 2, . . . is the discrete time instant index,

x(k) ∈ R n is the state, u(k) ∈ R l the input, y(k) ∈ R m the output, A(k), B(k), C(k) are time-varying matrices of appropriate sizes characterizing the nominal state-space model, w(k) ∈ R n , v(k) ∈ R m are mutually independent centered white Gaussian noises of covariance matrices Q(k) ∈ R n×n and R(k) ∈ R m×m ,
and the term Φ(k)θ represents actuator faults with a known matrix sequence Φ(k) ∈ R n×p and a constant (or piecewise constant with rare jumps) parameter vector θ ∈ R p . A typical example of actuator faults represented by the term Φ(k)θ is actuator gain losses. When affected by such faults, the nominal control term

B(k)u(k) becomes B(k)(I l -diag(θ))u(k) = B(k)u(k) -B(k)diag(u(k))θ
where I l is the l × l identity matrix, the diagonal matrix diag(θ) contains gain loss coefficients within the interval [0, 1], and Φ(k) ∈ R n×l (p = l) is, in this particular case, Φ(k) = -B(k)diag(u(k)).

(2)

A straightforward solution for the joint estimation of x(k) and θ is to apply the Kalman filter to the augmented system

x(k) θ(k) = A(k) Φ(k) 0 I p x(k-1) θ(k-1) + B(k) 0 u(k) + w(k) 0 y(k) = [C(k) 0] x(k) θ(k) + v(k).
However, to ensure the stability of the Kalman filter, this augmented system should be uniformly completely observable and uniformly completely controllable regarding the state noise [START_REF] Kalman | New methods in Wiener filtering theory[END_REF][START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF]. Notice that, even in the case of time invariant matrices A and C, the augmented system is time varying because of Φ(k), which is typically time varying. The uniform complete observability of an LTV system is defined as the uniform positive definiteness of its observability Gramian [START_REF] Kalman | New methods in Wiener filtering theory[END_REF][START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF]. In practice, it is not natural to directly assume properties (observability and controllability) of the augmented system (a more exaggerated way would be directly assuming the stability its Kalman filter, or anything else that should be proved!). Moreover, the augmented system is clearly not uniformly completely controllable regarding the state noise w(k), since the augmented states θ(k) are not controlled at all by w(k).

In contrast, in the present paper, the classical uniform complete observability and uniform complete controllability are assumed for the original system (1), in terms of the Gramian matrices defined for the [A(k), C(k)] pair and the [A(k), Q 1 2 (k)] pair. These conditions, together with a persistent excitation condition (see Assumption 3 formulated later), ensure the stability of the adaptive Kalman filter presented in this paper.

THE ADAPTIVE KALMAN FILTER

In the adaptive Kalman filter, the state estimate x(k|k) ∈ R n and the parameter estimate θ(k) ∈ R p are recursively updated at every time instant k. This algorithm involves also a few other recursively updated auxiliary variables:

P (k|k) ∈ R n×n , Υ(k) ∈ R n×p , S(k) ∈ R p×p
and a forgetting factor λ ∈ (0, 1). At the initial time instant k = 0, the initial state x(0) is assumed to be a Gaussian random vector

x(0) ∼ N (x 0 , P 0 ).

(3) Let θ 0 ∈ R p be the initial guess of θ, λ ∈ (0, 1) be a chosen forgetting factor, and ω be a chosen positive value for initializing S(k), then the adaptive Kalman filter consists of the initialization step and the recursion steps described below. Each part of this algorithm separated by horizontal lines will be commented after the algorithm description.

Initialization

P (0|0) = P 0 Υ(0) = 0 S(0) = ωI p (4a) θ(0) = θ 0 x(0|0) = x 0 (4b)
Recursions for k = 1, 2, 3, . . . θ(k-1)] is for the purpose of compensating the error caused by θ(k-1) = θ. This term is essential for the analysis of the properties of the adaptive Kalman filter in the following sections. It has also been introduced in the deterministic adaptive observer in [START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF], and its continuous time counterpart in [START_REF] Zhang | Adaptive observer for multiple-inputmultiple-output (MIMO) linear time varying systems[END_REF].

P (k|k-1) = A(k)P (k-1|k-1)A T (k) + Q(k) (5a) Σ(k) = C(k)P (k|k -1)C T (k) + R(k) (5b) K(k) = P (k|k -1)C T (k)Σ -1 (k) (5c) P (k|k) = [I n -K(k)C(k)]P (k|k -1) (5d) Υ(k) = [I n -K(k)C(k)]A(k)Υ(k-1) + [I n -K(k)C(k)]Φ(k) (5e) Ω(k) = C(k)A(k)Υ(k-1) + C(k)Φ(k) (5f) Λ(k) = [λΣ(k) + Ω(k)S(k-1)Ω T (k)] -1 (5g) Γ(k) = S(k-1)Ω T (k)Λ(k) (5h) S(k) = 1 λ S(k-1) - 1 λ S(k-1)Ω T (k)Λ(k)Ω(k)S(k-1) (5i) ỹ(k) = y(k) -C(k) A(k)x(k-1|k-1) + B(k)u(k) + Φ(k) θ(k-1) (5j) θ(k) = θ(k-1) + Γ(k)ỹ(k) (5k) x(k|k) = A(k)x(k-1|k-1) + B(k)u(k) + Φ(k) θ(k-1) + K(k)ỹ(k) + Υ(k)[ θ(k) -θ(k-1)]. (5l)

STABILITY OF THE ADAPTIVE KALMAN FILTER

Assumption 1. The matrices

A(k), B(k), C(k), Φ(k), Q(k), R (
(k), Q 1 2 (k)
] pair is uniformly completely controllable, in the sense of the uniform positive definiteness of the corresponding Gramian matrices. [START_REF] Kalman | New methods in Wiener filtering theory[END_REF][START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF].

2 Assumption 3. The signals contained in the matrix Φ(k) are persistently exciting in the sense that there exit an integer h > 0 and a real constant α > 0 such that, for all integer k ≥ 0, the matrix sequence Ω(k) driven by Φ(k) through the linear system (5e)-(5f), satisfies

h-1 s=0 Ω T (k + s)Σ -1 (k + s)Ω(k + s) ≥ αI p . (6)
2 Proposition 1. The matrices P (k|k), Υ(k), Σ(k), K(k), Ω(k) all have a finite upper bound, and Σ(k) has a strictly positive lower bound.

2

Proof. The proof based on classical results is quite straightforward. The recursive computations (5a)-(5d) for Σ(k), K(k), P (k|k) are identical to the corresponding part in the classical Kalman filter, hence like in the classical Kalman filter theory [START_REF] Kalman | New methods in Wiener filtering theory[END_REF][START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF], the boundedness of P (k|k) is ensured by the complete uniform observability of the [A(k), C(k)] pair and the complete uniform controllability of the [A(k), Q 1 2 (k)] pair stated in Assumption 2. As simple corollaries of this result and Assumption 1, the matrices Σ(k) and K(k) are also bounded. Equation (5b) implies Σ(k) > R(k). The covariance matrix R(k) is assumed to have a strictly positive lower bound, that is also a lower bound of Σ(k). Again under the complete uniform observability and the complete uniform controllability conditions, the homogenous part of the Kalman filter, corresponding to the homogenous system

ζ(k) = [I n -K(k)C(k)]A(k)ζ(k-1), (7)
is exponentially stable [START_REF] Kalman | New methods in Wiener filtering theory[END_REF][START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF]. This result implies that Υ(k) driven by bounded Φ(k) through (5e) is bounded. It is then follows from (5f) that Ω(k) is also bounded.

2 Notice that S(k) was missing in Proposition 1. It is the object of the following proposition. Proposition 2. The matrix S(k) recursively computed through (5i) has a finite upper bound and a strictly positive lower bound for all k ≥ 0.

2

Proof. In order to show the upper and lower bounds of S(k), let us first study another recursively defined matrix sequence

M (0) = S -1 (0) = ω -1 I p > 0 (8) M (k) = λM (k-1) + Ω T (k)Σ -1 (k)Ω(k). ( 9
)
It will be shown later that S(k) = M -1 (k), but for the moment this relationship is not used. According to Proposition 1, Ω(k) is upper bounded and Σ(k) has a strictly positive lower bound, hence Ω T (k)Σ -1 (k)Ω(k) is upper bounded. Then M (k) recursively generated from ( 9) with λ ∈ (0, 1) is also upper bounded. The lower bound of M (k) is investigated in the following.

Repeating the recursion of M (k) in ( 9) yields

M (k) = λ k M (0) + k-1 j=0 λ j Ω T (k -j)Σ -1 (k -j)Ω(k -j) For 0 ≤ k ≤ h, M (k) ≥ λ h M (0).
For k > h, let [k/h] denote the largest integer smaller than or equal to k/h, and break the sum k-1 j=0 into sub-sums of h terms, then

M (k) ≥ [k/h] i=1 (i-1)h+(h-1) j=(i-1)h+0 λ j Ω T (k -j)Σ -1 (k -j)Ω(k -j) ≥ [k/h] i=1 λ ih-1 (i-1)h+(h-1) j=(i-1)h+0 Ω T (k -j)Σ -1 (k -j)Ω(k -j) ≥ [k/h] i=1 λ ih-1 αI p
where the last inequality is based on Assumption 3 (persistent excitation). The geometric sequence sum in this last result can be explicitly computed, so that

M (k) ≥ αλ h-1 (1 -λ h[k/h] ) 1 -λ h I p ≥ αλ h-1 I p .
Therefore, M (k) has a strictly positive lower bound for either 0 ≤ k ≤ h or k > h. Now take the matrix inverse of both sides of equation ( 9), and apply the matrix inversion formula

(A+V T BV ) -1 = A -1 -A -1 V T (B -1 +V A -1 V T ) -1 V A -1 with A = λM (k-1), B = Σ -1 (k) and V = Ω(k), then M -1 (k) = 1 λ M -1 (k-1) - 1 λ M -1 (k-1)Ω T (k) λΣ(k) + Ω(k)M -1 (k-1)Ω T (k) -1 Ω(k)M -1 (k-1)
This recursion in M -1 (k) coincides exactly with that of S(k) as formulated in (5i). Moreover, as defined in (8), M (0) = S -1 (0), hence M -1 (k) = S(k) for all k = 0, 1, 2, . . . . It then follows from the already proved upper and lower bounds of M (k) that S(k) has also a finite upper bound and a strictly positive lower bound. 2 Define the state and parameter estimation errors

x(k|k) x(k) -x(k|k) (10) θ(k) θ -θ(k) (11)
The main result of this section is stated in the following proposition. Proposition 3. The mathematical expectations Ex(k|k) and E θ(k) tend to zero exponentially when k → ∞.

2

In other words, the state and parameter estimates of the adaptive Kalman filter converge respectively to the true state x(x) and to the true parameter θ in mean. It also means that the deterministic part of the error dynamic system, ignoring the random noise terms, is exponentially stable.

Proof. It is straightforward to compute from (1), ( 5j) and (5l) that andfrom (5j) and(5k) that θ

x(k|k) = A(k)x(k-1|k-1) + Φ(k) θ(k-1) + w(k) -K(k)ỹ(k) -Υ(k)[ θ(k) -θ(k-1)] = [I n -K(k)C(k)][A(k)x(k-1|k-1) + Φ(k) θ(k-1)] + Υ(k)[ θ(k) -θ(k-1)] + [I n -K(k)C(k)]w(k) -K(k)v(k),
(k) = θ(k-1) -Γ(k)ỹ(k) (12) = θ(k-1) -Γ(k)C(k)[A(k)x(k-1|k-1) + Φ(k) θ(k-1)] -Γ(k)C(k)w(k) -Γ(k)v(k). ( 13
)
Like in [START_REF] Zhang | Adaptive observer for multiple-inputmultiple-output (MIMO) linear time varying systems[END_REF], define

ξ(k) x(k|k) -Υ(k) θ(k). ( 14
)
Simple substitutions lead to

ξ(k) = [I n -K(k)C(k)][A(k)x(k-1|k-1) + Φ(k) θ(k-1)] + Υ(k)[ θ(k) -θ(k-1)] + [I n -K(k)C(k)]w(k) -K(k)v(k) -Υ(k) θ(k).
In this last result, according to (14), replace x(k-1|k-1) with x(k-1|k-1) = ξ(k-1)

+ Υ(k-1) θ(k-1), ( 15
) then ξ(k) = [I n -K(k)C(k)]A(k)ξ(k-1) + [I n -K(k)C(k)]A(k)Υ(k-1) + [I n -K(k)C(k)]Φ(k) -Υ(k) θ(k-1)] + [I n -K(k)C(k)]w(k) -K(k)v(k). The content of the curly braces {• • • } is zero, because Υ(k) satisfies (5e). Then ξ(k) = [I n -K(k)C(k)]A(k)ξ(k-1) + [I n -K(k)C(k)]w(k) -K(k)v(k). (16)
The noises w(k) and v(k) are assumed to have zero mean values (centered noises), then Eξ(k) = [I n -K(k)C(k)]A(k)Eξ(k-1). Like (7), this recurrent equation is exponentially stable.

Starting from

Eξ(0) = Ex(0) -Υ(0)E θ(0) = 0, it is recursively shown that Eξ(k) = 0 for all k ≥ 0. In (13) replace x(k-1|k-1) with (15), then θ(k

) = I p -Γ(k)C(k)A(k)Υ(k-1) -Γ(k)C(k)Φ(k) θ(k-1) -Γ(k)C(k)A(k)ξ(k-1) -Γ(k)C(k)w(k) -Γ(k)v(k) = [I p -Γ(k)Ω(k)] θ(k-1) -Γ(k)e(k) (17)
where Ω(k) is as defined in (5f), and

e(k) C(k)A(k)ξ(k-1) + C(k)w(k) + v(k).
(18) It was already shown Eξ(k) = 0 for all k ≥ 0, therefore Ee(k) = 0. Take the mathematical expectation at both sides of (17) and denote θ

(k) E θ(k), ( 19
) then θ(k) = [I p -Γ(k)Ω(k)] θ(k-1). ( 20
)
Before analyzing the convergence of θ(k) governed by ( 20), let us combine the two equations ( 5h) and (5i) into

S(k) = 1 λ [I p -Γ(k)Ω(k)]S(k-1). (21) Accordingly, M (k) = S -1 (k) satisfies 2 M (k) = λM (k-1)[I p -Γ(k)Ω(k)] -1 . ( 22
)
Let us define the Lyapunov function candidate

V (k) θ(k) T M (k) θ(k), (23)
it then follows from ( 20) and ( 22) that

V (k) = θ(k-1) T [I p -Γ(k)Ω(k)] T λM (k-1) θ(k-1) = λ θ(k-1) T M (k-1) θ(k-1) -λ θ(k-1) T Ω T (k)Γ T (k)M (k-1) θ(k-1). ( 24
)
By recalling (5h), (5g) and M (k-1) = S -1 (k-1), it yields

Ξ(k) Ω T (k)Γ T (k)M (k-1) (25) = Ω T (k)Λ(k)Ω(k) (26) = Ω T (k) λΣ(k) + Ω(k)S(k-1)Ω T (k) -1 Ω(k), (27)
which is a symmetric positive semidefinite matrix. Then (24) becomes simulated trials is generated, and all the histograms are depicted as a 3D illustration in Fig, 3. The histograms are normalized so that they are similar to probability density functions.

V (k) = λV (k-1) -λ θ(k-1) T Ξ(k) θ(k-1) ≤ λV (k-1), ( 28
) implying that V (k) = θ(k) T M (

  Recursions (5a)-(5d) compute the covariance matrix P (k|k) ∈ R n×n of the state estimate, the innovation covariance matrix Σ(k) ∈ R m×m and the state estimation gain matrix K(k) ∈ R n×m . These formulas are identical to those of the classical Kalman filter. Inspired by the recursive least square (RLS) estimator with an exponential forgetting factor, recursions (5e)-(5i) compute the parameter estimate gain matrix Γ(k) ∈ R p×m through the auxiliary variables Υ(k) ∈ R n×p , Ω(k) ∈ R m×p , S(k) ∈ R p×p . Equation (5j) computes the innovation ỹ(k) ∈ R m . Finally, recursions (5k)-(5l) compute the parameter estimate and the state estimate. Part of equation (5l), namely,x(k|k) ∼ A(k)x(k-1|k-1) + B(k)u(k) + K(k)ỹ(k) can be easily recognized as part of the classical Kalman filter, with the traditional prediction step and update step combined into a single step. The term Φ(k) θ(k -1) corresponds to the actuator fault term Φ(k)θ in (1a), with θ replaced by its estimate θ(k -1). The extra term Υ(k)[ θ(k)-

  k) and the input u(k) are upper bounded, Q(k) is symmetric positive semidefinite, and R(k) is symmetric positive-definite with a strictly positive lower bound, for all k ≥ 0. 2 Assumption 2. The [A(k), C(k)] pair is uniformly completely observable, and the [A

Fig. 2 .Fig. 3 .

 23 Fig. 1. Random switching index among 4 state-space models.

There exists a "forward" variant form of the state-space model, typically with x(k + 1) = A(k)x(k) + B(k)u(k) + w(k). While this difference is important for control problems, it is not essential for estimation problems, like the one considered in this paper. The form chosen in this paper corresponds to the convention that data are collected at k = 1,

2, 3, . . . and the initial state refers to x(0).

According to Proposition 2, S(k) is positive definite for all k ≥ 0, hence (21) implies that the matrix [Ip -Γ(k)Ω(k)] is invertible.

CONCLUSION

Unlike classical adaptive Kalman filters, which have been designed for state estimation in case of uncertainties about noise covariances, the adaptive Kalman filter proposed in this paper is for the purpose actuator fault diagnosis, through joint state-parameter estimation. The stability and minimum variance properties of the adaptive Kalman filter have been rigorously analyzed. Through LTV/LPV reformulation and approximations, this method for actuator fault diagnosis is also applicable to a large class of nonlinear systems.

is the part of the state estimation error x(k|k) due to the parameter estimation error θ(k) (this part would be zero if the parameter estimate θ(k) was replaced by the true parameter value θ, and then the adaptive Kalman filter (5) would be reduced to the standard Kalman filter). Therefore, the meaning of this proposition is that the remaining part of the state estimation error reaches its minimum variance if the Kalman gain K(k) is used. Proof. Following (16) where K(k) is replaced by L(k), and noticing that ξ(k-1), w(k) and v(k) are pairwise independent, compute the covariance matrix of ξ(k):

For shorter notations, let us denote

Define also

Because C(k)Π(k)C T (k) ≥ 0 and R(k) > 0 (R(k) is a positive definite matrix, see Assumption 1), H(k) is also positive definite, and thus invertible. Rearrange (33) as

The equivalence between ( 33) and ( 35) can be shown by first developing (35) and then by incorporating (34).

The matrix H(k) is positive definite, hence the first term in (35) (a symmetric matrix product) is positive semidefinite, that is,

When L(k) is chosen such that this inequality becomes equality, i.e., the first term of ( 35) is zero, cov[ξ(k)] reaches its minimum, since the other terms in (35) are independent of L(k). This optimal choice of L(k), denoted by L * (k), is

It remains to show that L * (k) is identical to the Kalman gain K(k) in order to complete the proof. Rewrite (33) while incorporating (34):

In the particular case

Assemble ( 32), ( 34), ( 36) and (39) together, for

(40d) These equations allow recursive computation of cov[ξ(k)|L * ] and the other involved matrices. It turns out that these recursive computations are exactly the same as those in (5a)-(5d), with Π(k), H(k), L * (k) and cov[ξ(k)|L * ] corresponding respectively to P (k|k -1), Σ(k), K(k) and P (k|k). It remains to show that these two recursive computations have the same initial condition, i.e. cov[ξ(0)|L * ] = P (0|0), in order to conclude L * (k) = K(k) for all k ≥ 0. Because Υ(0) = 0 as specified in (4a), and according to the definition of ξ(k) in ( 14),

6. NUMERICAL EXEMPLE Consider a piecewise constant system randomly switching within 4 third order (n = 3) state-space models with one input (l = 1) and 2 outputs (m = 2). Each of the 4 statespace models is randomly generated with the Matlab code (requiring the System Identification Toolbox):

zreal = rand(1,1)*1.2-0.6; pmodul = rand(1,1)*0.1+0.4; pphase = rand(1,1)*2*pi; preal = rand(1,1)-0.5; ssk = idss(zpk (zreal,[pmodul.*exp(1i*pphase) ... pmodul.*exp(-1i*pphase) preal],rand(1,1)+0.5,1));

In the simulation for k = 0, 1, 2, . . . , 1000, the the actual model at each instant k randomly switches among the 4 randomly generated state-space models, which are kept unchanged. The random switching sequence among the 4 state-space models is plotted in Figure 1. The input u(k) is randomly generated with a Gaussian distribution and the standard deviation equal to 2. The noise covariance matrices are chosen as Q(k) = 0.1I 3 and R(k) = 0.05I 2 for k ≥ 0. The matrix Φ(k) is as in (2) so that θ represents actuator gain loss. During the numerical simulation running from k = 0 to k = 1000, a gain loss of 50% at the time instant k = 500 is simulated, corresponding to a jump of θ (a scalar parameter) from 0 to 0.5. The result of parameter estimation by the adaptive Kalman filter is presented in Fig. 2. The above results are based on a single numerical simulation trial. In order to statistically evaluate the performance of the proposed method, 1000 simulated trials are performed, each corresponding to a different set of 4 statespace models and to a different random realization of the input and noises. At each time instant k, the histogram of the parameter estimation error based on the 1000


Qinghua Zhang
email: qinghua.zhang@inria.fr

Inverse Scattering for Electrical Cable Soft Fault Diagnosis Demonstrator Extended Abstract

Keywords: Fault diagnosis, electrical cable, inverse scattering, reflectometry

The reliability of electrical connections becomes a crucial issue, because more and more connecting cables are involved in modern engineering systems. In this demo, the tested cable is a twisted pair of 10 meters, with faults created during the demo by locally untwisting the cable. Easily reparable for the purpose of demo, such geometrical deformations are soft faults, since they do not cause impedance discontinuity in the cable. A network analyzer connected to one end of the cable collects spectral data, which are processed by an algorithm in real time to estimate the characteristic impedance profile distributed along the cable. The faults created during the demo are detected, localized and quantified in real time through the estimated characteristic impedance profile. As the impedance estimation algorithm does not make any assumption about the shape of the characteristic impedance profile, this method is particularly suitable for the diagnosis of multiple and distributed faults. The demo will be shown in real time with a true cable connected to a network analyzer. A video is available online.

INTRODUCTION

The fast development of electronic devices is accompanied by more and more connecting cables. Though the probability of failure is small for each connecting cable, the reliability of electrical connections becomes a crucial issue because of their large numbers [START_REF] Auzanneau | Wire troubleshooting and diagnosis: Review and perspectives[END_REF]. In this context, a promising technology for cable fault diagnosis is reflectometry, which consists in analyzing the reflection and the transmission of electrical waves observed at the ends of a cable. It has been reported that this technology is able to e ciently detect and to locate hard faults (severe impedance discontinuities, typically open circuit or short circuit). For soft faults, the problem is much more di cult, in particular when no impedance discontinuity is caused by such faults, as most methods for soft fault diagnosis are based on the detection of (relatively weak) impedance discontinuities [START_REF] Kafal | Locating faults with high resolution using single-frequency TR-MUSIC processing[END_REF]. Recently, the inverse scattering theory has been applied to the reflectometry technology for cable soft fault diagnosis [START_REF] Zhang | Inverse scattering for soft fault diagnosis in electric transmission lines[END_REF][START_REF] Tang | An inverse scattering approach to soft fault diagnosis in lossy electric transmission lines[END_REF][START_REF] Loete | Experimental validation of the inverse scattering method for distributed characteristic impedance estimation[END_REF]. By estimating the profile of the characteristic impedance distributed along a cable, this method is particularly e cient for the diagnosis of soft faults that do not imply impedance discontinuity. Without assuming any particular shape of the impedance profile of the monitored cable, it is well adapted to the case of multiple or spatially distributed faults. Its fast numerical computation is suitable for real time applications.

CABLE MODEL AND INVERSE PROBLEM

The telegrapher's equations are widely used as a mathematical model of electrical cables in engineering practice [START_REF] Paul | Analysis of multiconductor transmission lines[END_REF] :

@ @z v(t, z) + L(z) @ @t I(t, z) + R(z)i(t, z) = 0 (1a) @ @z i(t, z) + C(z) @ @t v(t, z) + G(z)v(t, z) = 0, ( 1b
)
where t represents the time, z is the longitudinal coordinate along the cable, v(t, z) and i(t, z) are respectively the voltage and the current in the cable at the time instant t and at the position z, R(z), L(z), C(z) and G(z) denote respectively the series resistance, inductance, capacitance and shunt conductance per unit length of the cable at the position z. As in typical cables R(z) and G(z) are very small (respectively at the order of 10 2 ohm and 10 9 siemens per meter), they are neglected in most studies for fault diagnosis.

In reflectometry, the geometrical coordinate z is often replaced by the wave propagation time from the position 0 to the position z, which is denoted by x and usually referred to as the coordinate in electrical distance. The transformation z ! x(z) is invertible, so is the inverse transformation x ! z(x). After this coordinate transformation, by neglecting the series resistance R(z) and the shunt conductance G(z), the telegrapher's equations (1) characterized by the 4 distributed parameters R(z), L(z), C(z), G(z) are then reduced to a simpler form characterized by a single distributed parameter, the characteristic impedance

Z 0 (x) , s L(z(x)) C(z(x)) . (2)
A common reflectometry instrument is the network analyzer that, when connected to one end of a cable, can measure its reflection coe cient r(f ) defined as

r(f ) , V (f, 0) Z ref I(f, 0) V (f, 0) + Z ref I(f, 0) (3)
where f is the frequency, V (f, 0) and I(f, 0) are respectively the Fourier transforms of v(t, 0) and i(t, 0) (the voltage v(t, z) and the current i(t, z) at the position z = 0), and Z ref is the reference impedance of the network analyzer (typically 50 ohm).

In the field of reflectometry, an important cable inverse scattering problem is the determination of the characteristic impedance profile Z 0 (x) along a cable from its reflection coe cient r(f ) measured at one end of the cable.

The word "inverse" refers to the "direct problem" that is the simulation of reflectometry measurements for a given profile of Z 0 (x). Quite often "inverse problems" are much more di cult than "direct problems".

A method for solving this cable inverse problem, based on the inverse scattering theory, has been presented in [START_REF] Zhang | Inverse scattering for soft fault diagnosis in electric transmission lines[END_REF] and experimentally validated in [START_REF] Loete | Experimental validation of the inverse scattering method for distributed characteristic impedance estimation[END_REF]. It is implemented in a software, ISTL (Inverse Scattering for Transmission Lines).

THE DEMO FOR REAL TIME FAULT DIAGNOSIS

The demo setup is illustrated in Figure 1. A two-wire cable (with parallel or twisted wires) is connected at one end to a network analyzer, which measures the reflection coe cient r(f ) of the cable. The reflection coe cient is processed by ISTL in real time to compute the characteristic impedance profile Z 0 (x), which is displayed on a screen. The tested cable of 10m long (a twisted pair) is shown in Figure 2. The values of r(f ) are measured for f =1,2, . . . , 500MHz. The characteristic impedance Z 0 (x) is computed for the tested cable discretized at 1000 points, without any assumption about the shape of the profile of Z 0 (x).

The demo repeats continuously a cycle of operations involving data acquisition, computation of Z 0 (x), and display of the curve of Z 0 (x). Each cycle takes about 1 second, which is fast enough for real time fault diagnosis. Soft faults are typically caused by chemical corrosion, overheating, mechanical wearing and geometrical deformation. In this demo, soft faults are artificially created by geometrical deformations (locally untwisting the twisted pair). Such deformations lead to local smooth variations of the cable characteristic impedance, and have the advantages of being reversible and easily operational in a real time demo.

The e↵ects of such created faults are displayed on the impedance curve with a delay of about 1 second. Figure 3 shows an example of locally untwisted cable, whose e↵ect on the impedance curve is indicated by the solid red circle. The demo will be shown in real time with a true cable connected to a network analyzer. A video is available online [START_REF] Zhang | ISTL video demo[END_REF].

Fig. 1 .

 1 Fig. 1. Illustration of the demo setup.

Fig. 2 .

 2 Fig. 2. A twisted pair cable of 10 meters.

Fig. 3 .

 3 Fig. 3. Soft fault created by untwisting the cable and its e↵ect on the impedance curve (in the solid red circle). The abscissa indicates the position along the cable, and the ordinate the characteristic impedance. The other circles indicate the e↵ects of two other similar faults.The demo will be shown in real time with a true cable connected to a network analyzer. A video is available online[START_REF] Zhang | ISTL video demo[END_REF].


Gary A Soo

A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws

The finite difference methods of Godunov, Hyman, Lax and Wendroff (two-step), MacCormack, Rusanov, the upwind scheme, the hybrid scheme of Harten and Zwas, the antidiffusion method of Boris and Book, the artificial compression method of Harten, and Glimm's method, a random choice method, are discussed. The methods are used to integrate the one-dimensional Eulerian form of the equations of gas dynamics in Cartesian coordinates for an inviscid, nonheat-conducting fluid. The test problem was a typical shock tube problem. The results are compared and demonstrate that Glimm's method has several advantages.

INTRODUCTION

In the past few years numerous methods have been developed for solving systems of nonlinear hyperbolic conservation laws ut + F(U)., = o, (1) where F is a smooth function on Rn of U and the derivative matrix DF(U) possesses strictly nonlinear or linearly degenerate eigenfunctions and eigenvalues in the sense of Lax [16]. In the terminology of gas dynamic this corresponds to systems having modes which admit respectively, discontinuities of shock type, or contact discontinuities, only. For example, a single equation Ut + g(u)., = 0, is strictly nonlinear if g(u).,., =I= 0 and is linearly degenerate if g(u) = Cu, where C is a constant.

There are several difficulties in solving such systems numerically. Inherent in any finite difference scheme is an assumption on the regularity of the solution. Typically such schemes produce oscillations behind a shock. All finite difference schemes have numerical diffusion, dispersion, or both due to the truncation error. As a result there will be diffusion across the contact discontinuity with each time step, causing the contact discontinuity to be smeared in the course of the calculation.

A shock is also smeared, though less dramatically. Finally, there is often difficulty in obtaining an accurate approximation to the smooth parts of the solution, particularly at points where certain higher order derivatives fail to exist.

The finite difference methods of Godunov [6], Hyman [START_REF] Hyman | On robust and accurate methods for the calculation of compressible fluid flows, I[END_REF], Lax and Wendroff [START_REF] Richtmyer | A Survey of Difference Methods for Non-Steady Fluid Dynamics[END_REF], MacCormack [START_REF] Maccormack | Proceedings of the Second International Conference on Numerical Methods in Fluid Dynamics[END_REF], Rusanov [START_REF] Rusanov | [END_REF], the upwind scheme [START_REF] Morton | Difference Methods for Initial-Value Problems[END_REF], and the hybrid scheme of Harten and Zwas [12] will be discussed. Glimm's method, a random choice method, will also be considered. Glimm's method at first appearance may not seem to fit in the category of the other methods. However, it shares several features with the above difference schemes. All of these methods are fixed grid methods having the same grid structure. These methods require about the same amount of storage. All of the above finite difference schemes are finite difference approximations to the derivative arising in the conservation laws and can treat in principle an arbitrary system of conservation laws. Similarly Glimm's method can treat an arbitrary system of conservation laws. See Harten and Sod [START_REF] Harten | A generalized version of Glimm's method, to appear[END_REF].

The finite difference schemes as well as Glimm's method share one other feature, namely, a family of waves interacting can be handled automatically. However, the degree of resolution to which such interactions are represented is not the same for all methods.

Most finite difference methods when applied to problems with discontinuities produce oscillations behind the shock. Von Neumann and Richtmyer [START_REF] Neumann | [END_REF] developed an artificial viscosity term which was introduced into the Lagrangian form of the equations of gas dynamics. The goal of the artificial viscosity was to reduce the oscillations while allowing the shock transition to occupy only a few mesh points and having negligible effect in the smooth regions. All other forms of artificial viscosity are variations of the one introduced by von Neumann and Richtmyer. The one used in this survey due to Lapidus [START_REF] Lapidus | [END_REF] is no exception. However, it is suited to the conservation laws and can be added to an existing method as a fractional step. One important feature of Lapidus' form of artificial viscosity is that a variant of it can be applied to general conservation laws. Artificial viscosity will not spread a contact discontinuity (see section on artificial viscosity or Sod [25]). The smearing of a contact discontinuity is due to the truncation error of the scheme and the smearing of a shock is due to the truncation error of the scheme as well as the artificial viscosity. The spreading of the contact discontinuity being much more severe, O(n 1 /(Hll) for a qth order accurate scheme, where n is the number of time steps. See Harten [START_REF] Harten | The Method of Artificial Compression[END_REF].

Recently methods have been developed for correcting the smearing of shocks and contact discontinuities. The "flux corrected transport" method of Boris and Book [I] is used to obtain a high degree of resolution without oscillations. The artificial compression method of Harten [START_REF] Harten | The Method of Artificial Compression[END_REF][START_REF] Harten | [END_REF][START_REF] Harten | The artificial compression method for shocks and contact discontinuities[END_REF] is designed to sharpen results in regions containing shocks and contact discontinuities.

The inviscid, non-heat-conducting equations of gas dynamics can be written in conservation form with the flux F having the form discussed above. However, there is a large class of problems for which little is known. This is the case where F is either nonsmooth or not strictly nonlinear. Examples are the Buckley-Leverett equation, in the case of a single equation (see Buckley and Leverett [27]) and the problem of an exothermic reacting gas, in the case of systems. In these cases special methods have to be considered. For the case of an exothermic reacting gas see Chorin [4] and Sod [START_REF] Son | A numerical model of unsteady combustion phenomena[END_REF].

The only "high" -order method considered in this survey is the hybrid scheme of Harten and Zwas. Fourth-and sixth-order methods were originally considered for the equations of gas dynamics. However, with all of the high-order methods using the appropriate high-order artificial viscosity term, the coefficient of numerical diffusion was so large that the time step became prohibitively small. The only way in which a solution could be obtained with these high-order methods was to use the artificial viscosity due to Lapidus, which is third order. With the use of a low-order artificial viscosity term, all benefits obtained by using a high-order method are lost.

In the following sections a brief discussion of the methods is given, their solution to a sample one-dimensional problem is compared, and the merits of the methods are discussed. Due to the nonstandardness of Glimm's method, as well as the difficulty in programming, its acceptance as an effective and efficient numerical tool may be restricted. For this reason the equations used by Glimm's method are derived and a flow chart for its implementation is given. For further details on these methods see Sod [25].

Basic Equations

The one-dimensional equations of gas dynamics may be written in the (conservation) form:

o 1 p + o.,(pu) = o, o 1 m + ox((m 2 /p) + p) = 0, o 1 e + 8,((m/p)(e + p)) = 0, (2a) (2b) (2c)
where p is the density, u is the velocity, m = pu is momentum, p is pressure, and e is energy per unit volume. We may write e = pE + !pu 2 , where E is the internal energy per unit mass. Assume the gas is polytropic, in which case

E = pj(y -1) p,
where y is a constant greater than one. Furthermore, from (3) we have p = A(S)pY, where S denotes entropy.

Equations (I}-(3) may be written in vector form where

U 1 + F(U), = 0, and
F(U) = [ (m 2 /~ + p ]• (m/p)(e + p) (3) (4) (5)
In order to deal with solutions containing shocks we write the equations in integral form, which is obtained by integrating Eqs. (2a)-(2c) (or Eq. ( 5)) over any region in the upper half of the (x, t) plane and applying Green's theorem

u~+l = u." -(Llt/Llx)(F~+' -FM') t t t-t-t t-!- u~+} = 1(un + u ") -(Lltj2Llx)(F~' -F n) z+ 2 2 r-1 t H 1 t u;+l = ut -(Lltj2Llx)(F;H -F;_ 1 )
+ i((a;+l + a,n)(u;+l -ut)

e<;n = w(Llt/Llx)(u + c),n u;+l = u,n -sgn(u,n)(Llt(Llx)(Gt-G;+s<ul) - (

DESCRIPTION OF THE METHODS

The methods of Godunov [6], Lax and Wendroff (two-step) [START_REF] Richtmyer | A Survey of Difference Methods for Non-Steady Fluid Dynamics[END_REF], MacCormack [START_REF] Maccormack | Proceedings of the Second International Conference on Numerical Methods in Fluid Dynamics[END_REF], Rusanov [START_REF] Rusanov | [END_REF], and the upwind difference scheme [START_REF] Morton | Difference Methods for Initial-Value Problems[END_REF] have been widely used and no benefit can be obtained by describing them here. Hence, these schemes will merely be listed in Table I. The remaining methods under consideration will be briefly discussed.

Consider the nonlinear system of Eqs. (5). Divide time into intervals of length Lit and let Llx be the spatial increment. The solution is to be evaluated at time n Llx, where n is a nonnegative integer at the spatial increments i Llx, i = 0, ±1, ±2, ... , and at time (n + -t) Lit at (i + t) Llx. Let ut approximate U(i Llx, n Lit) and uf;-l approximate U((i + t) Llx, (n + t) Llt).

-

Glimm's Method

Recently the random choice method introduced by Glimm [7] has been developed for hydrodynamics by Chorin [3]. The method is a two-step method. To find u.":l and thus define the method, consider system (5) along with the piecewise constant initial data

V(x, n Lit) = uf+l , x ~ (i + t) Llx, x < (i + t) Llx. ( 9
)
This defines a sequence of Riemann problems. If Lit < Llx/2(1 u I + c), where c is the local sound speed, the waves generated by the different Riemann problems will not interact. Hence the solution v(x, t) to the Riemann problem can be combined into a single exact solution. Let fn be an equidistributed random variable which is given by the Lebesgue measure on the interval [ -t, H Define [START_REF] Harten | The artificial compression method for shocks and contact discontinuities[END_REF] At each time step, the solution is approximated by a piecewise constant function. The solution is then advanced in time exactly and the new values are sampled. The method depends on the possibility of solving the Riemann problem exactly and inexpensively. Chorin [3] (see also Sod [START_REF] Son | The Computer Implementation of Glimm's Method[END_REF]) modified an iterative method due to Godunov [5] which will be described below.

Consider system (5) with the initial data = Sr = (pr, Ur ,pr), X< 0,

X~ 0. ( 11
)
The solution at later times looks like (see [15]) Using this iterative method we first evaluate p* in the stateS*. Define the quantity (12)

If the left wave is a shock, using the jump condition U 1 [p] = [pu ], we obtain (13

)
where U 1 is the velocity of the left shock and p* is the density in the portion of S* adjoining the left shock. Similarly, define the quantity (

If the right wave is a shock, using the jump conditions Ur[p] = [pu], we obtain (15) where Ur is the velocity of the right shock and p* is the density in the portion of S* adjoining the right shock.

In either of the two cases ((12) or ( 13) for M 1 and ( 14) and (15) for Mr) we obtain where Mr = (PrPr) 112 4>(P*/Pr), Mz = (pzpz) 112 4>(P*/pz),

y -I I -x = 2ylf2 1 _ xt-tY ' X~ I, x~I.
Upon elimination of u* from ( 12) and ( 14) we obtain ( Pr Pz )/( I 1 )

p* = Uz-Ur +-M +-M -M + M . r l l r (16a) (16b) (17) (18)
Equations ( 16a), (16b), and ( 18) represent three equations in three unknowns for which it can be seen that there exists a real solution. Upon choosing a starting value p* 0 (or M 1 ° and M/), we iterate using Eqs. (16a), (16b), and [START_REF] Maccormack | Proceedings of the Second International Conference on Numerical Methods in Fluid Dynamics[END_REF]. For details of the starting values see Chorin [3] and Sod [25]. After p*, M 1 , and Mr have been determined we may obtain u* by eliminating p* from Eqs. ( 12) and ( 14), [START_REF] Richtmyer | A Survey of Difference Methods for Non-Steady Fluid Dynamics[END_REF] from which the complete solution of the Riemann problem can be determined using the jump conditions for shock waves and the isentropic law and constancy of the Riemann invariants for the rarefaction waves. See the Appendix.

The finite difference method due to Godunov [6] in Table I is for the Eulerian form of the equations of gas dynamics. The method developed by Godunov [5] for the Lagrangian form is also a two-step method where the second step is the second half step in Table I. However, the values of u"/.N and p"/N are replaced by u* [START_REF] Richtmyer | A Survey of Difference Methods for Non-Steady Fluid Dynamics[END_REF] and

p* (18) from the Riemann problem at i + !-• - Artificial Viscosity
In the methods of Godunov, MacCormack, and Lax and Wendroff (two-step), an artificial viscosity term was added. The artificial viscosity term used was introduced by Lapidus [START_REF] Lapidus | [END_REF]. It has the advantage that it is very easy to add to an existing scheme and it retains the high-order accuracy of the scheme. Let fi"/+ 1 be the approximation at time (n + 1) Lit obtained by any one of the above schemes. This value is replaced by the new approximation where Ll'ii;n = ii;n -fi~1 and v is an adjustable constant.

Equation ( 20) is a fractional step for the numerical solution of the following diffusion equation

It is shown (see Ladidus [START_REF] Lapidus | [END_REF]) that this new difference scheme (obtained by adding the artificial viscosity) satisfies the same conservation law that the previous equation did. The values of the constant v used varied from method to method. This is discussed in the section on numerical results. This artificial viscosity was not added in the smooth regions and not applied to the first component (mass) equation of (5). The artificial viscosity is not applied to the mass equation so as to minimize the amount of mass diffused across the contact discontinuity.

Harten's Corrective Method of Artificial Compression

In this section we discuss the Artificial Compression Method (ACM) developed by Harten [9]. This method is designed to be used in conjunction with an already existing finite difference scheme. The purpose of this method is to sharpen the regions which contain discontinuities whether shocks or contact discontinuities.

Only the basic idea of the ACM will be discussed for the case of a single conservation law. Let u(x, t) be a solution of the conservation law Ut+ f(u)x = 0 [START_REF] Rusanov | [END_REF] which contains a discontinuity (uL(t), uR(t), S(t)), where uL and uR are the values on the left and right of the jump and S is the speed of the discontinuity. The discontinuity is either a shock or a contact. Assume, without loss of generality that at any given timet the solution u does not take on any values between uL(t) and uR(t). Consider the function g(u, t) with properties

g(u, t) sgn [uR(t) -uL(t)] > 0 g(u, t) ~ 0 for u E (uL(t), uR(t)), for u if= (uL(t), uit)).
This function g will be called an artificial compression flux.

It can be seen that u is also a solution of the conservation law

Ut + (f(u) + g(u, t))x = 0. ( 22
) (23) (24)
By [START_REF] Son | [END_REF] we see that when u is smooth Eq. ( 24) is identical with Eq. ( 21) and the shock speed S(t) remains the same. Finally it is observed (from ( 22)) that if (uL , uR , S) is a shock or contact for Eq. ( 21) then it is a shock for the modified equation [START_REF] Son | A numerical model of unsteady combustion phenomena[END_REF]. The artificial compression method solves the modified equation ( 24) rather than the original equation [START_REF] Rusanov | [END_REF]. For a complete discussion of the implementation of the method see Harten [START_REF] Harten | [END_REF].

Let u~+ 1 represent the approximate solution vector to (5) obtained by using any one of the above first order difference methods. In solving the modified system(analogous to (24)) we use operator splitting. We first define the difference representation g; of the artificial compression flux g, (25) where and [START_REF] Neumann | [END_REF] where k refers to the kth component of the fi, S~+t = fi~~~ -fi~kl. Let Si+i represent the vector whose kth component is sgn (S~H). Then the difference scheme which applies the ACM to the given solution u~+ 1 is

un+l = iin+l-~(g. -g. ) , 2Lix <+1 z-1 (27a) -n+l Lit (Gn Gn ) = ui -2 LI i+ • - i-, , X - - ( 27b
)
where G~H = gin -g~1 -I g~+ 1 -g;n I Si+-t, applied componentwise. See Harten [START_REF] Harten | The Method of Artificial Compression[END_REF].

The method of artificial compression is designed for first-order schemes and cannot be applied directly to higher-order schemes. The idea of ACM is based on the existence of a viscous profile. See Harten [START_REF] Harten | [END_REF]. Higher-order schemes introduce other flux terms so that one obtains different (non physical) speeds of propagation.

Self-Adjusting Hybrid Schemes

The idea of self-adjusting hybrid schemes was introduced by Harten and Zwas [12]. Consider a nonoscillatory first-order scheme L 1 and a kth-order (k ~ 2) scheme Lk , There are many choices for such schemes. The scheme chosen here is discussed in Harten [START_REF] Harten | The artificial compression method for shocks and contact discontinuities[END_REF]. Taking k = 2 we choose MacCormack's scheme and by adding the artificial viscosity term (34) to MacCormack's scheme we obtain the first-order scheme.

The hybridized scheme becomes for system (5)

(35) (36)
The stability condition for the first-order scheme is Lit

3 1 1 2 max(l u I+ c) Llx ~ T,
this being stricter than the stability condition for MacCormack's scheme. So this is the stability condition for the hybrid scheme.

It remains to describe how switch 0 is chosen. There are many possible choices, the one selected is described in Harten [START_REF] Harten | The artificial compression method for shocks and contact discontinuities[END_REF]. Let Lli+t = Pi+ 1 -p;. Define (37) = 0, otherwise.

In this case p = 1 and E > 0 is chosen as a measure of negligible variation in the density p. We define the switch 0 by Since in areas which contain a discontinuity the hybrid scheme is about first order we may apply the artificial compression method discussed above. However, the ACM must not be used in smooth regions. For this purpose the switch is used again, i.e., Eq. (27b) may be replaced with n+1 -n+l

.:1 t (On Gn on Gn )

U; = U; -2 A i+ l i+ > -i-1 i-> • 41X - - • - ( 38
)
10

Antidiffusion Method of Boris and Book

In this section we shall discuss briefly the antidiffusion method developed by Boris and Book [1]. The purpose of this special technique known as ":flux correction" is to achieve high resolution without oscillations.

It can be shown that a first-order difference scheme can be represented by an equation of the form

Ut + j(u)x = Llt [g(u, LltjLlx) Ux]x, ( 39
)
where g(u, LltjLlx) is the coefficient of the diffusion term.

The basis of the antidiffusion method is to use a stable modification of a diffusive difference scheme. Let the original scheme be represented by (39). The modification is represented by The second step is the antidiffusion step, which is unstable by itself since it approximates the backward heat equation. We define It can be seen that if g(u, LltfLlx) -r(u, Llt/Llx) ;;, 0, (43) then the combined scheme AL is stable. However, (43) places more of a restriction on LltjLlx than the stability condition for L.

Ut + f(u)x = Llt [(g(u, LltJLlx) -r(u, LltJLlx)) Ux]x,
We chose for L the two-step Lax-Wenfroff scheme. Following Boris et al.

The parameter YJ is the diffusionjantidiffusion coefficient. The stability condition is max(l u I + c)(L1tjL1x) :(; 1.

Hyman's Predictor-Corrector Method

In [START_REF] Hyman | On robust and accurate methods for the calculation of compressible fluid flows, I[END_REF] Hyman describes a predictor-corrector type scheme. The spatial derivatives are approximated by a second-order difference operator while the time derivative (or time integrator) uses the improved Euler scheme. The improved Euler scheme combines a first-order explicit predictor with a second-order trapazoidal rule corrector.

For stability and to insure proper entropy production an artificial viscosity term is added. The artificial viscosity term used is similar to that used by Rusanov [START_REF] Rusanov | [END_REF].

The scheme is given by where

u~+l = ut -(L1tj2)(DF~H + Pt), q,;H = (l/4L1x)(cx;+ 1 + cxt)(u;+l -ut), cxt = (u + c)? ,
and c is the local sound speed.

(

) (45b) 45a
The stability of the scheme depends on the number of applications of the corrector (45b) and on S. We took as the stability condition max (I u I + c)(L1tjL1x) :(; 1.

In order to maintain stability, the artificial viscosity must not be completely removed in the smooth regions. However, it can be reduced in these regions by using a type of switch. The one chosen was suggested by Hyman [START_REF] Hyman | On robust and accurate methods for the calculation of compressible fluid flows, I[END_REF]. Replace cf>~H in (45a) by f3cpf+l where

R -. l fJ-3' =I, if ex~+ I > cxi n + (Lix/3)
otherwise.

This type of switch greatly reduces the smearing of the contact discontinuity as well as the shock wave. This switch is a type of artificial compression. (region 2) some of the derivatives of the fluid quantities may not be continuous The point x 3 is the position that an element of fluid initially at x 0 has reached by time t. Point x 3 is called a contact discontinuity. It is seen that across a contact discontinuity the pressure and the normal component of velocity are continuous. However, the density and the specific energy are not continuous across a contact discontinuity. Point x 4 is the location of the shock wave (moving to the right). Across a shock all of the quantities (p, m, e, and p) will in general be discontinuous.

THE SHOCK TUBE PROBLEM

In the study of the above numerical methods the following test problem was considered: PI = 1.0, PI = 1.0, ui = 0.0, p 5 = 0.125, p 5 = 0.1, and u 5 = 0. The ratio of specific heats y was chosen to be 1.4. In all of the calculations Llx = 0.01. For the Rusanov scheme the value of w (see Table I) was taken to be 1.0. In the scheme of Boris and Book the parameter 7J was taken to be 0.125. For Hyman's scheme the value of I) was taken to be 0.8. The constant in the artificial viscosity term v was taken to be 1.0 in all but one case. Also the value of a (see Table I) was taken to be 0.9.

In Glimm's original construction a new value of g was chosen for each grid point i and each time level n. The practical effect of such a choice with finite Llx is disastrous since our initial data is not close to constant (which was an assumption made by Glimm). In fact, if g is chosen for each i and n, it is possible that a state will propagate to the left and to the right and thus create a spurious state. An improvement due to Chorin [3] is to choose gn only once per time step (hence the subscript n). The details of the method of selection of the random number are found in Chorin [3] and Sod [START_REF] Son | The Computer Implementation of Glimm's Method[END_REF]. Figure 4 indicates the results using the first-order accurate Godunov scheme. The corners at the endpoints of the rarefaction wave are rounded. The constant state between the contact discontinuity and the shock has not been fully realized. The transition of the contact discontinuity occupies seven to eight zones while the transition of the shock occupies five to six zones.

Figure 5 indicates the results using the Godunov scheme with artificial compression. It should be noted that for this case the constant in the artificial term was taken to be 2.0 to ensure that the solution before application of artificial compression was oscillation free. For the artificial compression cannot be applied in the presence of oscillations. The corners at the endpoints of the rarefaction wave are still rounded, since the artificial compression method is not applied in smooth regions. There is a slight undershoot at the right corner of the rarefaction. Also there are oscillations at the contact discontinuity. The transition of the contact discontinuity occupies three to four zones while the transition of the shock occupies only one to two zones.

Figure 6 shows the results of the two-step Lax-Wendroff scheme. There are very slight overshoots at the contact discontinuity and more noticeable overshoots at the shock. The rarefaction wave is quite accurate. The corners at the endpoints of the rarefaction are only slightly rounded. The transition of the contact discontinuity occupies six to eight zones while the shock wave occupies four to six zones. It is observed that the plots in Fig. 6 are quite similar to those in Fig. 7 obtained by MacCormack's method.

Figure 7 represents the results of the second-order MacCormack scheme. There are slight overshoots at the contact discontinuity and more noticeable overshoots at the shock wave. The rarefaction wave is quite accurate. The corners at the endpoints of the rarefaction are only slightly rounded. The transition of the contact discontinuity occupies seven to eight zones while the transition of the shock occupies five to six zones.

Figure 8 represents the first-order accurate Rusanov scheme. The contact discontinuity is barely visible in the density profile. The corners at the end points of the rarefaction wave are extremely rounded. The constant state between the contact discontinuity and the shock wave is barely existent. The transition of the contact discontinuity occupies 14-16 zones and the transition of the shock occupies 6-8 zones. This scheme is extremely diffusive.

Figure 9 represents the Rusanov scheme with artificial compression. The results with artificial compression are greatly improved. The corners at the endt;>Oints of the rarefaction wave are still rounded since the artificial compression method is not applied in this area. The constant state between the contact discontinuity and the shock is much more visible. The transition of the contact discontinuity occupies two to three zones while that of the shock wave occupies only one to two zones.

Figure 10 represents the upwind difference scheme. It is observed that between the left constant state and the left endpoint of the rarefaction wave is a shock (discontinuity). This is clearly a non physical solution. This is a result of the method used -----+---------+----------- .so . ~5 : .uc

~-~-'• ---+-
"+• -••-----+•--------1---•-------1
..... ;:.cu Fro. 9. Rusanov's first-order method with ACM. to stabilize the scheme, by using centered differences for the pressure term in the momentum equation.

Figure 11 shows the results of the Glimm scheme. The shock wave and the contact discontinuity have been computed with infinite resolution, i.e., the number of zones over which the variation occurs is zero. Due to the randomness of the method the positions of the shock and the contact discontinuity are not exact. However, on the average their positions are exact. The corners at the end points of the rarefaction wave are perfectly sharp. It is observed that the rarefaction is not smooth (due to the randomness), yet it is extremely close to the exact solution. The constant states are perfectly realized.

The Glimm scheme requires between two and three times as much time (see below) as the other finite difference schemes tested. However, the Glimm scheme requires far less spatial grid points for the same resolution. This is displayed in Table II, where nine interior grid points are used. All details are visible.

The Glimm scheme on the average is conservative. One other check on the accuracy is to use the conservation laws (mass, momentum, and energy). For example, the total mass is evaluated by a r+ is the Riemann invariant (c/(y-I))+ (u/2), where c is the local sound speed. .so In Table HI the values of the total mass, momentum, and energy are displayed. The mass and the energy are seen to be conserved on the average, i.e., there are fluctuations but they are contained within a small interval. The momentum is seen to increase linearly on the average (allowing for fluctuations).

Q. = L p(i Ax) Ax.
Figure 12 shows the results of the antidiffusion method of Boris and Book applied to the two-step Lax-Wendroff scheme. There is a slight overshoot at the right corner of the rarefaction. The rarefaction wave is very accurately computed. The corners at the end points of the rarefaction are only slightly rounded. The constant state between the contact discontinuity and the shock wave is only partially realized. The transition of the contact discontinuity occupies five to seven zones and the transition of the shock occupies one to two zones. The resolution is much better than the two-step Lax-Wendroff scheme alone (see Fig. 6).

Figure 13 represents the hybrid scheme ( 35) and (36) of Harten and Zwas. The solution is free of oscillations. The corners at the endpoints of the rarefaction wave are only slightly rounded. The constant state between the contact discontinuity and the shock is only partly realized. The transition of the contact discontinuity occupies eight to nine zones and the transition of the shock occupies five to six zones.

Figure 14 represents the hybrid scheme of Harten and Zwas with the use of artifical compression. Since the artifical compression is not applied in smooth regions the rarefaction is the same as in Fig. 13. The transition of the contact discontinuity occupies three to four zones and the transition ofthe shock wave occupies two to three zones.

Figure 15 represents the results of Hyman's predictor-<:orrector scheme, where the corrector has been applied once. The solution is oscillation free. The corners at the endpoints of the rarefaction are almost perfectly sharp. The constant states between the rarefaction and the contact discontinuity and between the contact discontinuity are extremely well defined. The transition of the contact discontinuity occupies six to eight zones while the transition of the shock occupies three to four zones.

The timing results for all of the methods are listed in Table IV. The times are for 100 spatial grid points. The only substantial difference in timing is between Glimm's scheme and the other finite difference schemes. For Glimm's scheme requires between two and three times as much time. However, Glimm's scheme can give the same resolution with far less points (as seen in Table 11). From the point of view of the least number of grid points per desired resolution, the Glimm scheme can be seen to be much faster.

CONCLUSIONS

Of all the finite difference schemes tested, without the use of corrective procedures, Godunov's and Hyman's methods produced the best results.

It is obvious from the figures that the Glimm scheme gives the best resolution of the shocks and contact discontinuities. Glimm's scheme is at best first-order accurate (see Chorin [3]) so that boundary conditions are easily handled.

It is possible that the rarefaction wave obtained by Glimm's method can be smoothed out by a type of averaging. This is presently being considered.

The hybrid method of Harten and Zwas combines first-and high-order schemes in such a way as to extract the best features of both. The high-order scheme produces better approximations to the smooth parts of the flow.

The corrective procedures of Boris and Book and Harten improve the resolution of a given scheme. The artificial compression method being restricted to first-order schemes except when used in conjunction with the hybrid type schemes produces far better results than the antidiffusion method of Boris and Book. Both methods are easily added to existing programs (as a subroutine). The antidiffusion method requires slightly more storage than the artificial compression method since the former must retain two time levels of information for the computation of intermediate results (Eq. (44c)).

A major disadvantage of the antidiffusion method of Boris and Book, the hybrid scheme of Harten and Zwas, and the artificial compression method of Harten is that there are a number of parameters to be chosen, which depend on the given problem. In the antidiffusion method the coefficient of diffusionjantidiffusion must be chosen. The value of this parameter can greatly affect the results. In the hybrid scheme a tolerance must be chosen for the automatic switch which is taken to be a measure of negligible variation in entropy or density for example. This tolerance depends on the given problem. In the artificial compression method a test must be included to locate the rarefaction (and other smooth regions). Many of the standard tests fail to work well enough for the use of artificial compression. ~ ; With the method described for solving the Riemann problem in the Glimm scheme, it can only be used for the equations of gas dynamics in rectangular coordinates. It is possible to generalize Glimm's method to other coordinate systems and different equations. See Harten and Sod [START_REF] Harten | A generalized version of Glimm's method, to appear[END_REF].

The applicability of Glimm's method to other geometries has only just started to be explored. One successful application is to the equations of gas dynamics for a cylindrically or spherically symmetric flow. See Sod [START_REF] Son | [END_REF].

The usefulness of Glimm's method in the analysis of reacting gas flow is given in Chorin [4] and Sod [START_REF] Son | A numerical model of unsteady combustion phenomena[END_REF]. In Chorin's paper examples are given of deflagration and detonation waves, with infinite and finite rates of reaction.

Since the Godunov iteration procedure provides an exact solution (up to the tolerence prescribed to terminate the iteration) to the Riemann problem, the onedimensional shock tube problem can be viewed as a "sitting duck" for Glimm's or any other numerical method using solutions of Riemann problems. However, it is observed that due to the randomness the results of Glimm's method are not exact.

The true value of a numerical method should be measured by its ability to solve multidimensional problems. For such problems, Glimm's method or others which depend on the solution of one-dimensional Riemann problems may not have any advantage over the other methods considered here. This is presently being studied. Equations ( 46)-( 62) are the key to the programming of Glimm's method. For a summary see the flow chart (Fig. 16).

  + J m dt = 0, Jmdx+ JC(m 2 /p)+p)dt=O, J e dx + J ((m/p)(e + p)) dt = 0.

  Lltf2Llx)(S;+, -s;_,), • where G = (pu, pu 2 , u(e + p))r, S = (0, p, O)T, and s(u) = -1 if ut> 0 if u(' < 0 • a = max(l u 1 + c)Llt(Llx, where c denotes the local sound speed.

  Fig. l, where S 1 and S 2 are either a shock or a centered rarefaction wave. The region S * is a steady state. The lines 1 1 and 1 2 separate the states. The contact discontinuity dx/dt = u* separates the state S* into two parts with possibly different values of p*, but equal values of u* and p*.

FIG. 1 .

 1 FIG. 1. Solution of a Riemann problem.

  u':'+l = fi':'+l + v Lit Ll'[l Ll'fin+I I • Ll'fin+I]

  as not to violate the conservation, hybridize L 1 and Lk through their numerical fluxes. Define the hybrid operator L by (30) where (31) ()i+i is a scalar quantity (called a switch) which satisfies 0 ~ {}i+i ~ I. At discontinuities the automatic switch is such that {} R::i I. Hence at the discontinuities the hybrid scheme is essentially the nonoscillatory first-order scheme. Equation (30) can be written in the form (32) so that if () is o(LfxP) where the solution is smooth, then for p ~ k -1 we have (33)

  r is a positive function. One can introduce the antidiffusion term by operator splitting. The first step consists of solving Ut + f(u)., = 0, (41) with the original difference scheme, say fi7+ 1 = Luin• Then in the second step let A be a difference operator approximating the diffusion equation Ut + Llt[r(u, LltjLlx) u,]., = 0. (42)

  .Ji+t = YJ U;+l -U; ,

Figure 2 1 Reg:Loo 5 FIG. 2 . 5 FIG. 3 .

 215253 Figure 2 represents the initial conditions in a shock tube. A diaphragm at x 0 separates two regions (regions 1 and 5) which have different densities and pressures.The two regions are in a constant state. The initial conditions are p 1 > p 5 , p 1 > p 5 , and u 1 = u 5 = 0; i.e., both fluids are initially at rest. At time t > 0 (see Fig.3) the diaphragm is broken. Consider the case before any wave has reached the left or right boundary. Points x 1 and x 2 represent the location of the head and tail of the rarefaction wave (moving to the left). Although the solution is continuous in this region

FIG. 5 .

 5 FIG. 4. Godunov's method.

FIG. 6 .

 6 FIG.6. Two-step Lax-Wendroffmethod.

  --+--------~

FIG. 7 .FIG. 8 .

 78 FIG. 7. MacCormack's method.

XFIG. 11 .

 11 FIG. 11. Glimm's method.

  FIG. 14. Hybrid method with ACM.

FIG. 15 .

 15 FIG. 15. Hyman's predictor-corrector method.

( 3 )

 3 FIG. 16. Flow chart of Glimm's method.

TABLE li Profiles

 li Obtained by Glimm's Method for Nine Interior Grid Points

	X	p	u	p	e	ra +
	0.1	1.000	0.000	1.000	2.500	2.958
	0.2	1.000	0.000	1.000	2.500	2.958
	0.3	0.869	0.164	0.822	2.363	2.958
	0.4	0.426	0.927	0.303	1.778	2.958
	0.5	0.426	0.927	0.303	1.778	2.958
	0.6	0.426	0.927	0.303	1.778	2.958
	0.7	0.426	0.927	0.303	1.778	2.958
	0.8	0.266	0.927	0.303	2.853	3.624
	0.9	0.125	0.000	0.100	2.000	2.646

TABLE Ill Total

 Ill Mass, Momentum, and Energy for Glimm's Scheme

				g ,;
			FIG. 12. Antidiffusion method.
	tjilt		Qp	Qm	Q.
			0.547	0.018	2.213
	2		0.550	0.019	2.217
	3		0.554	0.032	2.218
	4		0.550	0.039	2.219
	5		0.552	0.047	2.223
	6		0.550	0.059	2.222
	7		0.549	0.070	2.221
	8		0.550	0.079	2.224
	9 10	X	0.545 0.546	0.090 0.097	2.232 2.247
	11		0.548	0.110	2.258
	12		0.545	0.119	2.267
	13		0.549	0.122	2.266
	14		0.552	0.136	2.266
	15		0.549	0.143	2.269
	16		0.553	0.149	2.275
	17		0.550	0.158	2.266
	18		0.546	0.164	2.267
	19		0.550	0.178	2.267
	20		0.543	0.190	2.272

g ' "T FIG. 13. Hybrid method.

TABLE IV Running

 IV Time per Time Step (in sec) on CDC 6600

	Schemes	Without ACM	WithACM
	Godunov	0.226	0.247
	Lax-Wendroff	0.226
	MacCormack	0.224
	Rusanov	0.224	0.240
	Upwind	0.225
	Glimm	0.364
	Anti diffusion	0.242
	Hybrid	0.258	0.269
	Hyman	0.276

a Times include computation of exact solution, calls to printing, and plotting routines, which were the same for all cases.

ACKNOWLEDGMENTS The author would like to thank Professor Alexandre Chorin, Mr. Phillip Colella, Dr. Amiram Harten, and Dr. James M. Hyman for a number of helpful discussions and comments. The author would also like to thank one of the editors, Dr. Raymond Chin, for his assistance. This work was supported in part by the National Science Foundation, Grant MCS76-07039, the Energy Research and Development Administration under Contract No. EY-76-C-02-3077 at the Courant Institute, and the Energy Research and Development Administration under Contract No. W-7405-Eng-48 at Lawrence Livermore Laboratory.

(52) By substitution of (52) into (49) and solving for u we obtain _ 2

By substitution of (53) into (52) c is obtained; by substitution of (52) into the definition of c and solving for p we obtain


Johanna D García-Saldaña
email: johanna@mat.uab.cat

Armengol Gasull
email: gasull@mat.uab.cat

A theoretical basis for the Harmonic Balance Method

Keywords:

series Fixed point theorem

The Harmonic Balance Method provides a heuristic approach for finding truncated Fourier series as an approximation to the periodic solutions of ordinary differential equations. Another natural way for obtaining these types of approximations consists in applying numerical methods. In this paper we recover the pioneering results of Stokes and Urabe that provide a theoretical basis for proving that near these truncated series, whatever is the way they have been obtained, there are actual periodic solutions of the equation. We will restrict our attention to one-dimensional non-autonomous ordinary differential equations, and we apply the obtained results to a concrete example coming from a rigid cubic system.

Introduction and main results

Consider the real non-autonomous differential equation x = X(x, t), [START_REF] Gasull | Limit cycles for rigid cubic systems[END_REF] where the prime denotes the derivative with respect to t, X : Ω × [0, 2π ] → R is a C 2 -function, 2π -periodic in t, and Ω ⊂ R is a given open interval.

There are several methods for finding approximations to the periodic solutions of [START_REF] Gasull | Limit cycles for rigid cubic systems[END_REF]. For instance, the Harmonic Balance Method (HBM), recalled in Section 2.1, or simply the numerical approximations of the solutions of the differential equations. In any case, from all the methods we can get a truncated Fourier series, namely a trigonometric polynomial, that "approximates" an actual periodic solution of the equation. The aim of this work is to recover some old results of Stokes and Urabe that allow the use of these approximations to prove that near them there are actual periodic solutions and also provide explicit bounds, in the infinity norm, of the distance between both functions. To the best of our knowledge these results are rarely used in the papers dealing with HBM.

When the methods are applied to concrete examples one has to manage the coefficients of the truncated Fourier series that are rational numbers which renders the subsequent computations more difficult. See the example of Section 4. At this point we introduce in this setting a classical tool that as far as we know has never been used in this type of problems: we approximate all the coefficients of the truncated Fourier series by suitable convergents of their respective expansions in continuous fractions. This is done in such a way that by using these new coefficients we obtain a new approximate solution that is essentially at the same distance to the actual solution as the starting approximation. With this method we obtain trigonometric polynomials with nice rational coefficients that approximate the periodic solutions.

Before stating our main result, and following [START_REF] Stokes | On the approximation of nonlinear oscillations[END_REF][START_REF] Urabe | Galerkin's procedure for nonlinear periodic systems[END_REF], we introduce some concepts. Let x(t) be a 2π -periodic C 1 -function; we will say that x(t) is noncritical with respect to [START_REF] Gasull | Limit cycles for rigid cubic systems[END_REF]

if 2π 0 ∂ ∂x X x(t), t dt = 0.
(

) 2
Notice that if x(t) is a periodic solution of (1) then the concept of noncritical is equivalent to the one of being hyperbolic; see [START_REF] Lloyd | A note on the number of limit cycles in certain two-dimensional systems[END_REF].

As we will see in Lemma 2.1, if x(t) is noncritical w.r.t. Eq. ( 1), the linear periodic system

y = ∂ ∂x X x(t), t y + b(t)
has a unique periodic solution y b (t) for each smooth 2π -periodic function b(t). Moreover, once X and x are fixed, there exists a constant M such that

y b ∞ M b 2 , ( 3
)
where as usual, for a continuous 2π -periodic function f ,

f 2 = 1 2π 2π 0 f 2 (t) dt, f ∞ = max x∈R f (x) and f 2 f ∞ .
Any constant satisfying (3) will be called a deformation constant associated to x and X. Finally, consider

s(t) := x (t) -X x(t), t . ( 4
)
We will say that x(t) is an approximate solution of (1) with accuracy S = s 2 . For simplicity, if S > S, we also will say that x(t) has accuracy S. Notice that actual periodic solutions of (1) have accuracy 0; in this sense, the function s(t) measures how far is x(t) from being an actual periodic solution of (1).

The next theorem improves some of the results of Stokes [START_REF] Stokes | On the approximation of nonlinear oscillations[END_REF] and Urabe [START_REF] Urabe | Galerkin's procedure for nonlinear periodic systems[END_REF] in the one-dimensional setting. More concretely, in those papers they prove the existence and uniqueness of the periodic orbit when 4M 2 K S < 1. We present a similar proof with the small improvement 2M 2 K S < 1. Moreover our result gives, under an additional condition, the hyperbolicity of the periodic orbit. Theorem 1.1. Let x(t) be a 2π -periodic C 1 -function such that -it is noncritical w.r.t. Eq. (1) and has M as a deformation constant, -it has accuracy S w.r.t. Eq. [START_REF] Gasull | Limit cycles for rigid cubic systems[END_REF].

Given I := [min {t∈R} x(t) -2M S, max {t∈R} x(t) + 2M S] ⊂ Ω, let K < ∞ be a constant such that max (x,t)∈I×[0,2π ] ∂ 2 ∂x 2 X(x, t) K .
Therefore, if , then the periodic orbit x * (t) is hyperbolic, and its stability is given by the sign of this integral.

2M 2 K S <
Once some approximate solution is guessed, for applying Theorem 1.1 we need to compute the three constants appearing in its statement. In general, K and S can be easily obtained. Recall for instance that s 2 , when s is a trigonometric polynomial, can be computed from Parseval's theorem. On the other hand, M is much more difficult to estimate. In Lemma 2.3 we give a result useful for computing it in concrete cases, that is different from the approach used in [START_REF] Stokes | On the approximation of nonlinear oscillations[END_REF][START_REF] Urabe | Galerkin's procedure for nonlinear periodic systems[END_REF][START_REF] Urabe | Numerical computation of nonlinear forced oscillations by Galerkin's procedure[END_REF].

Assuming that a non-autonomous differential equation has a hyperbolic periodic orbit, the results of [START_REF] Urabe | Galerkin's procedure for nonlinear periodic systems[END_REF] also guarantee that, if given a suitable trigonometric polynomial r(t) of a sufficiently high degree, we can apply the first part of Theorem 1.1. Intuitively, while the value of the accuracy S goes to zero when we increase the degree of the trigonometric polynomial, the values M and K remain bounded. Thus, at some moment, it holds that 2M 2 K S < 1.

In Section 4 we apply Theorem 1.1 to study and localize the limit cycles of the rigid cubic system ẋ = -y +

x 10 1 -x -10x 2 , ẏ = x + y 10 1 -x -10x 2 .
In polar coordinates it is written as ṙ = r/10cos(θ )r 2 /10cos 2 (θ)r 3 , θ = 1, or equivalently,

r = dr dt = 1 10 r - 1 10 cos(t)r 2 -cos 2 (t)r 3 , ( 5
)
and it has a unique positive periodic orbit; see also [START_REF] Gasull | Limit cycles for rigid cubic systems[END_REF]. Notice that we have renamed θ as t. We prove the following: As we will see, in this example we will find computational difficulties to obtain the third approximation given by the HBM. Therefore we will get it first by numerically approaching the periodic solution, then by numerically computing the first terms of its Fourier series and finally by using the continuous fractions approach to simplify the values appearing in our computations. We also will see that the same approach works for other concrete rigid systems.

Similar examples for second-order differential equations have also been studied in [START_REF] Urabe | Numerical computation of nonlinear forced oscillations by Galerkin's procedure[END_REF].

Preliminary results

This section contains some technical lemmas that are useful for proving Theorem 1.1 and for obtaining in concrete examples the constants appearing in its statement. We also include a very short overview of the HBM adapted to our interests. See [START_REF] Mickens | Oscillations in Planar Dynamic Systems[END_REF] for a more general point of view on the HBM.

As usual, given A ⊂ R, 1 A : R → R denotes the characteristic function of A: the function takes the value 1 when x ∈ A, and the value is 0 otherwise.

Lemma 2.1. Let a(t) and b(t) be continuous real 2π -periodic functions. Consider the non-autonomous linear ordinary differential equation x = a(t)x + b(t). (6)

If A(2π ) = 0, where A(t) := t 0 a(s) ds, then for each b(t) Eq. ( 6) has a unique 2π -periodic solution x b (t) := 2π 0

H(t, s)b(s) ds, where the kernel H(t, s) is given by the piecewise function

H(t, s) = e A(t) 1 -e A(2π ) e -A(s) 1 [0,t] (s) + e A(2π )-A(s) 1 [t,2π ] (s) . (7) Moreover x b ∞ 2π max t∈[0,2π ] H(t, •) 2 b 2 .
Proof. Since ( 6) is linear, its general solution is

x(t) = e A(t) x 0 + t 0 b(s)e -A(s) ds . (8)
If we impose that the solution is 2π -periodic, i.e. x(0) = x(2π ), we get

x 0 = e A(2π ) 1 -e A(2π ) 2π 0 b(s)e -A(s) ds, ( 9
) then (8) becomes x b (t) = e A(t) 1 -e A(2π ) e A(2π ) 2π 0 b(s)e -A(s) ds + 1 -e A(2π ) t 0 b(s)e -A(s) ds = e A(t)
1e A(2π ) e A(2π

) 2π t b(s)e -A(s) ds + t 0 b(s)e -A(s) ds = 2π 0 H(t, s)b(s) ds.
Therefore, the first assertion follows. On the other hand, by the Cauchy-Schwarz inequality,

x b (t) 2π 0 H 2 (t, s) ds 2π 0 b 2 (s) ds.
Hence,

x b ∞ 2π max t∈[0,2π ] H(t, •) 2 b 2 .
This completes the proof. 2

Corollary 2.2. A deformation constant M associated with x and X is

M := 2π max t∈[0,2π ] H(t, •) 2 ,
where H is given in

(7) with A(t) = t 0 ∂ ∂x X(x(t), t) dt.
Now we prove a technical result that will allow us to compute in practice deformation constants. In fact we will find an upper bound of M that will avoid the integration step needed in the computation of the norm • 2 . First, we introduce some notations.

Given a function A : [0, 2π ] → R, a positive number , and a partition t i = ih with i = 0, 1, . . . , N, of the interval [0, 2π ], where h = 2π /N, we consider the function L : [0, 2π ] → R given by the continuous linear piecewise function joining the points

(t i , A(t i ) -). Notice that L(t) = N-1 i=0 L i (t)1 I i , where I i = [t i , t i+1 ] and L i (t) = A(t i+1 ) -A(t i ) h (t -t i ) + f (t i ) := - 1 2 (α i t + β i ).
We will say that L is an adequate lower bound of A if it holds that L(t) < A(t) for all t ∈ [0, 2π ].

Clearly, smooth functions always have adequate functions that approach them.

For each m = 0, 1, . . . , N we define the function

Ψ m (t) := m-1 i=0 J i + λ 2 N-1 i=m-1 J i + 1 -λ 2 e β m α m e α m t -e α m t m , ( 10
)
where λ = e A(2π ) , and

J i := t i+1 t i e -
N = √ 2π |1 -λ| max t∈[0,2π ] e A(t) N-1 m=0 Ψ m (t)1 I m (t). Proof. Recall that from Lemma 2.1, x b ∞ M b 2 ,
where

M := 2π max t∈[0,2π ] H(t, •) 2 .
Thus, we will find an upper bound of M. Since

H(t, s) = e A(t) 1 -e A(2π ) e -A(s) 1 [0,t] (s) + e A(2π )-A(s) 1 [t,2π ] (s) , it holds that H(t, •) 2 = 1 √ 2π e A(t)
|1 -λ|

G(t)
where

G(t) := t 0 e -2A(s) ds + λ 2 2π t e -2A(s) ds < t 0 e -2L(s) ds + λ 2 2π t e -2L(s) ds, because L(t) < A(t) for all t ∈ [0, 2π ]. Assume that t ∈ I m . Thus, t 0 e -2L(s) ds = m-1 i=0 J i + t t m e -2L m (s) ds, 2π t e -2L(s) ds = N-1 i=m J i + t m+1 t e -2L m (s) ds = N-1 i=m-1 J i - t t m
e -2L m (s) ds.

Therefore, for t ∈ I m ,

G(t) < m-1 i=0 J i + λ 2 N-1 i=m-1 J i + 1 -λ 2 t t m e α m s+β m ds = Ψ m (t).
As a consequence, for t ∈ [0, 2π ],

G(t) < N-1 m=0 Ψ m (t)1 I m (t),
and the result follows. 2

Remark 2.4. Notice that the above lemma provides a way for computing a deformation constant where there is no need of computing integrals. This will be very useful in concrete application, where the primitive of e -2A(t) is not computable, and so Corollary 2.2 is difficult to apply for obtaining M.

In the next result, which introduces the constant K appearing in Theorem 1.1, D • denotes the topological interior of D. Lemma 2.5. Consider X as in [START_REF] Gasull | Limit cycles for rigid cubic systems[END_REF]. Let D be a closed interval, and let x(t) be a 2π -periodic C 1 -function such that {x(t):

t ∈ R} ⊂ D • . Define R(z, t) := X x(t) + z, t -X x(t), t - ∂ ∂x X x(t), t z ( 11
)
for all z such that {x(t) + z: t ∈ R} ⊂ D. Then

(i) |R(z, t)| K 2 |z| 2 , (ii) |R(z, t) -R(z, t)| K max(|z|, |z|)|z -z|, where K := max (x,t)∈D×[0,2π ] ∂ 2 ∂x 2 X(x, t) .
Proof. (i) By using Taylor's formula, for each t it holds that

X x(t) + z, t = X x(t), t + ∂ ∂x X x(t), t z + 1 2 ∂ 2 ∂x 2 X ξ(t), t z 2 for some ξ(t) ∈ x(t), x(t) + z . Therefore R(z, t) = 1 2 ∂ 2 ∂x 2 X ξ(t), t |z| 2 K 2 |z| 2 ,
as we wanted to prove.

(ii) From Rolle's theorem, for each fixed t it follows that there exists η(t

) ∈ z, z such that R(z, t) -R(z, t) ∂ ∂ z R η(t), t |z -z|.
Applying again this theorem, but now to ∂ ∂ z R, and by noticing that

∂ ∂ z R(z, t)| z=0 = 0, we obtain ∂ ∂ z R η(t), t ∂ 2 ∂ z 2 R ω(t), t η(t) = ∂ 2 ∂x 2 X ω(t), t η(t) K η(t) ,
where ω(t) ∈ 0, η(t) . Note also that η(t) max |z|, |z| .

Hence, the result follows combining the three inequalities. 2

The Harmonic Balance Method

In this subsection we recall the HBM adapted to the setting of one-dimensional 2π -periodic nonautonomous differential equations.

We are interested in finding periodic solutions of the 2π -periodic differential equation ( 1), or equivalently, periodic functions which satisfy the following functional equation

F x(t) := x (t) -X x(t), t = 0. ( 12
)
Recall that any smooth 2π -periodic function x(t) can be written as its Fourier series

x(t) =
F y N (t) = A 0 2 + ∞ m=1 A m cos(mt) + B m sin(mt) ,
where A m = A m (r, s) and B m = B m (r, s), m 0, with r = (r 0 , r 1 , . . . , r N ) and s = (s 1 , . . . , s N ). The

HBM consists of finding values r and s such that

A m (r, s) = 0 and B m (r, s) = 0 for 0 m N.

(

) 13
The above set of equations is usually a very difficult nonlinear system of equations, and for this reason in various works including [START_REF] Mickens | Oscillations in Planar Dynamic Systems[END_REF] and the references therein, only small values of N are considered. We also remark that in general the coefficients of y N (t) and y N+1 (t) do not coincide at all.

Notice that Eq. ( 13) is equivalent to The hope of the method is that the trigonometric polynomials found using this approach are "near" actual periodic solutions of the differential equation [START_REF] Gasull | Limit cycles for rigid cubic systems[END_REF]. In any case, as far as we know, the HBM for small N is only a heuristic method that sometimes works quite well.

To end this subsection, we want to comment on a main difference between the non-autonomous case treated here and the autonomous one. In this second situation the periods of the searched periodic orbits, or equivalently their frequencies, are also treated as unknowns. The method then works similarly; see again [START_REF] Mickens | Oscillations in Planar Dynamic Systems[END_REF].

Proof of the main result

Proof of Theorem 1.1. As a first step we prove the following result: consider the nonlinear differential equation

z = X z + x(t), t -X x(t), t -s(t), ( 14
)
where s(t) is given in [START_REF] Stokes | On the approximation of nonlinear oscillations[END_REF]. A 2π -periodic function z(t) is then a solution of ( 14) if and only if z(t)+ x(t) is a 2π -periodic solution of (1). This is a consequence of the following equalities:

z(t) + x(t) = X z(t) + x(t), t -X x(t), t -s(t) + X x(t), t + s(t)
= X z(t) + x(t), t .

By using the function

R(z, t) = X z + x(t), t -X x(t), t - ∂ ∂x X x(t), t z,
introduced in Lemma 2.5, Eq. ( 14) can be written as

z = ∂ ∂x X x(t), t z + R(z, t) -s(t). ( 15
)
Let P be the space of 2π -periodic C 0 -functions. To prove the first part of the theorem it suffices to see that Eq. ( 15) has a unique C 1 , 2π -periodic solution z * (t), which belongs to the set

N = z ∈ P: z ∞ 2M S .
To prove this last assertion, we will construct a contractive map T : N → N . Because N is a complete space with the • ∞ norm, its fixed point will be a continuous function in N that will satisfy an integral equation, equivalent to (15). Finally we will see that this fixed point is in fact a C 1 function, and it satisfies Eq. (15).

Let us define T . If z ∈ N , then T (z) is defined as the unique 2π -periodic solution of the linear

differential equation y = ∂ ∂x X x(t), t y + R z(t), t -s(t).
Notice that this map is well-defined, by Lemma 2.1, because x(t) is noncritical w.r.t. Eq. ( 1). Thus,

z 1 satisfies z 1 = ∂ ∂x X x(t), t z 1 + R z(t), t -s(t).
Let us prove that T maps N into N and that it is a contraction. By Lemmas 2.1 and 2.5 and the hypotheses of the theorem

T (z) ∞ = z 1 ∞ M R z(•), • -s(•) 2 M R z(•), • 2 + S M R z(•), • ∞ + S M K 2 z 2 ∞ + S M 2K M 2 S 2 + S < 2M S,
where we have used in the last inequality that 2M 2 K S < 1.

To show that T is a contraction on N , take z, z ∈ N and denote by

z 1 = T (z), z1 = T (z). Then z 1 = ∂ ∂x X x(t), t z 1 + R z(t), t -s(t), z 1 = ∂ ∂x X x(t), t z1 + R z(t), t -s(t).
Therefore,

(z 1 -z1 ) = ∂ ∂x X x(t), t (z 1 -z1 ) + R z(t), t -R z(t), t .
Again by Lemmas 2.1 and 2.5 and the hypotheses of the theorem,

T (z) -T (z) ∞ = z 1 -z1 ∞ M R z(•), • -R z(•), • ∞ M K max z ∞ , z ∞ z -z ∞ 2M 2 K S z -z ∞ ,
as we wanted to prove, because recall that 2M 2 K S < 1. Thus, the sequence of functions {z n (t)} defined as

z n+1 (t) = ∂ ∂x X x(t), t z n+1 (t) + R z n (t), t -s(t),
with any z 0 (t) ∈ N , and z n+1 (t) chosen to be periodic, converges uniformly to some function x * (t) ∈ N . In fact we also have that

z n+1 (t) = z n+1 (0) + t 0 ∂ ∂x X x(w), w z n+1 (w) + R z n (w), w -s(w) dw. Therefore, x * (t) = x * (0) + t 0 ∂ ∂x X x(w), w x * (w) + R x * (w), w -s(w) dw.
We know that x * (t) is a continuous function, but from the above expression we obtain that it is indeed of class C 1 . Therefore x * (t) is a periodic solution of (15) and is the only one in N , as we wanted to see.

To prove the hyperbolicity of x * (t), it suffices to show that 2π 0 ∂ ∂x X x * (t), t dt = 0, and study its sign; see [START_REF] Lloyd | A note on the number of limit cycles in certain two-dimensional systems[END_REF]. We have that, fixed t,

∂ ∂x X x * (t), t = ∂ ∂x X x(t), t + ∂ 2 ∂x 2 X ξ(t), t x * (t) -x(t) ,
for some ξ(t) ∈ x * (t), x(t) . Therefore, since we have already proved that |x * (t)

-x(t)| < 2M S, ∂ ∂x X x(t), t - ∂ ∂x X x * (t), t 2K M S. Then 2π 0 ∂ ∂x X x(t), t dt - 2π 0 ∂ ∂x X x * (t), t dt 4π K M S < 2π M ,
and the result follows because by hypothesis the first integral is, in absolute value, bigger than 2π /M. 2

Applications

In this section we apply Theorem 1.1 for proving the existence and uniqueness of hyperbolic limit cycles, in a suitable region, of some planar rigid systems, which after some transformations can be converted into differential equation of the form [START_REF] Gasull | Limit cycles for rigid cubic systems[END_REF]. More precisely, we study some concrete cases of the family of rigid cubic systems ẋ = -y -

x a + bx + x 2 , ẏ = x -y a + bx + x 2 , ( 16
)
already considered in [START_REF] Gasull | Limit cycles for rigid cubic systems[END_REF]. In that paper it is proved that (16) has at most one limit cycle, and when it exists, it is hyperbolic. With our point of view we will find an explicit approximation of the limit cycle; see Proposi- We want to find an approximation of the periodic solution of ( 5), which we will use for applying Theorem 1.1.

First attempt: the HBM. Following it, and according to Section 2.1, consider the equation F r(t) = r (t) -

1 10 r + 1 10 cos(t)r 2 + cos 2 (t)r 3 = 0, (17)
which is clearly equivalent to [START_REF] Urabe | Galerkin's procedure for nonlinear periodic systems[END_REF].

Searching for a solution of the form r(t) = r 0 and imposing that the first harmonic of Hence r 0 = √ 5/5 ≈ 0.4472135954 is the first-order solution given by the HBM. We obtain that the positive approximate solution is r = √ 5/5. For applying the second-order HBM we search for an approximation of the form r(t) = r 0 + r 1 cos(t) + s 1 sin(t).

The vanishing of the coefficients of 1, cos(t) and sin(t) in F (r(t)), provides the nonlinear system f (r 0 , r 1 , s 1 ) :=

r 1 = 0, g(r 0 , r 1 , s 1 ) := 9 4 r 2 0 r 1 - 5 8 r 3 1 + 3 8 r 1 s 2 1 + 1 10 r 2 0 + 3 40 r 2 1 + 1 40 s 2 1 - 1 10 r 1 + s 1 = 0, h(r 0 , r 1 , s 1 ) := 3 4 r 2 0 s 1 + 3 8 r 2 1 s 1 + 1 8 s 3 1 + 1 20 r 1 s 1 - 1 10 s 1 -r 1 = 0.
Doing the resultants Res( f , g, r 0 ) and Res( f , h, r 0 ) we obtain respectively 1775r

1 3 + 525r 1 s 1 2 + 240r 1 2 + 20s 1 2 -132r 1 -400s 1 -8 = 0, 105r 1 2 s 1 + 35s 1 3 + 8r 1 s 1 + 80r 1 -4s 1 = 0.
Repeating the resultant between these last two equations with respect to r 1 we have 17 150 000s The approximate real solution of this equation is s1 = -0.0196567414, and then we have the respective approximate solutions r0 = 0.4471066159, r1 = -0.0009814101.

For our purposes we can consider simpler rational approximations of r0 , r1 and s1 , with maintain- ing a similar accuracy. For finding these rational approximations, we seek them by performing the continued fraction expansion of these values. For instance, r0 = [0, Doing all the computations needed to apply Theorem 1.1 we get that we are not under its hypotheses. Therefore we need to continue with the third-order HBM.

Performing the third-order approach we obtain five algebraic polynomial equations that we omit for the sake of simplicity. Unfortunately, neither using the resultant method as in the previous case, nor using the more sophisticated tool of Gröbner basis, our computers are able to obtain an approximate solution to start our theoretical analysis.

A numerical approach. First, we search for a numerical solution of (5) by using the Taylor series method. From this approximation we compute, again numerically, its first Fourier terms obtaining r(t) =

with accuracy 0.00298, quite similar to the one of r(t). Note that (18) is precisely the approximation of the periodic solution of (5) stated in Proposition 1.2.

Proof of Proposition 1.2. We already know that the accuracy of r(t) is S := 0.003.

  cos(mt) + b m sin(mt) , ) cos(mt) dt, and b m = ) sin(mt) dt, for all m 0. Hence, it is natural to try to approximate the periodic solutions of the functional equation (12) by using truncated Fourier series, i.e. trigonometric polynomials.Let us describe the HBM of order N. Consider a trigonometric polynomial y N (t) = cos(mt) + s m sin(mt) with unknowns r m = r m (N), s m = s m (N) for all m N. Compute then the 2π -periodic function F (y N (t)). It has also an associated Fourier series

2π 0 F 0 F

 00 y N (t) cos(mt) dt = 0 and 2π y N (t) sin(mt) dt = 0, for 0 m N.

tion 1 . 2 .

 12 In order to simplify the computations we first consider the case a = -b = 1/10, which in polar coordinates is written as (5): )r 2cos 2 (t)r 3 .

3 k=0r

 3 k cos(kt) + s k sin(kt), where r 0 = 0.4483561517, r 1 = -0.0024133439, s 1 = -0.0193837572,r 2 = -0.0037463296, s 2 = -0.0220176517, r 3 = -0.0012390886, s 3 = 0.0003784656.The accuracy of r(t) is 0.00289. If we take a new simpler approximation, using again some convergents of r k and s k , we obtain

  2L(s) ds = α i t i+1e α i t i .

	t i+1 t i e -2L i (s) ds =	e β i α i

e Lemma 2.3. Let L be an adequate lower bound of A, where A is the function given in Lemma 2.1. Consider the functions Ψ m (t), with m = 0, 1, . . . , N -1. Therefore, also following the notation introduced in that lemma, it holds that x b ∞ N b 2 , where

  2, 4, 4, 2, 2, 2, 4, 2, 1, 1], giving the convergents 1/2, 4/9, 17/38, 38/85, . . . . Similarly r1 gives -1/1018, -1/1019, . . . and s1 gives -1/50, -1/51, . . . . At this point we have the following new candidate to be an approximation

	of the periodic solution
		r(t) =	1 2	-	1 1018	cos(t) -	1 50	sin(t).
	Its accuracy w.r.t. Eq. (5) is
	S = r (t) -	1 10	r(t) +	1 10	cos(t)r(t) 2 + cos 2 (t)r(t) 3	2	≈ 0.046.

  To apply Theorem 1.1 we will compute M and K .

	2 891 685 439 72 733 752 000	-	347 888 350 813 299 559 1 778 094 556 332 494 400	t -	561 179 36 756 720	cos(t) -	685 338 551 8 000 712 720	sin(t)
	-	757 058 717 48 004 276 320	cos(2t) -	40 221 206 418 131 273 447 836 421 760	sin(2t) -	2 923 231 576 974 475	cos(3t)

First we calculate

A(t) = t 0 ∂ ∂r X(r(t), t). A(t) =

Acknowledgments

The authors are partially supported by an MCYT/FEDER grant number MTM2008-03437 and by a CIRIT grant number 2009SGR 410.

Since 2π /M ≈ 0.9, the hyperbolicity of r * (t) follows by applying the second part of the theorem. 2

Notice that the example of the system (16) that we have studied is a = λ and b = -λ with λ = 1/10. With the same techniques we see that the same function r(t) given in the statement of Proposition 1.2 is an approximation of the unique periodic orbit of the system when |λ -1/10| < 1/500, which also satisfies (19).