Fabien Delaleux
email: fabien.delaleux@u-pec.fr

Olivier Riou

Vincent Guiheneuf

Pierre-Olivier Logerais

Jean-Félix Durastanti

Outil de contrôle non destructif pour la quantification du vieillissement des concentrateurs solaires

Dans cet article, nous présentons une méthode de quantification du vieillissement de miroirs verre des héliostats par leur émissivité apparente. La méthode est raccordée aux propriétés spectrales du verre. Elle présente l'avantage d'être sans contact et permet de distinguer des altérations équivalentes à 3 ans de durée de vie standard. Cet outil permet de cartographier la surface étudiée en émissivité apparente puis de convertir ces valeurs en variation de propriétés optiques dans le visible (absorbance, réflectance ou transmittance) suivant la fonction attendue.

Nomenclature

t température d'émission, °C L thermosignal, OS LWIR long wavelength IR MWIR mid wavelength IR SMIR small wavelength IR Symboles grecs ελ émissivité spectrale ε Δλ émissivité apparente Indices et exposants app apparent env environnement

Introduction

L'objectif de ce travail est de concevoir un outil de contrôle non destructif et non intrusif pour l'évaluation du vieillissement des systèmes solaires à concentration. Celui-ci impose la définition préalable des indicateurs associés puis la constitution de moyens de mesurage permettant leur quantification in situ. Cette approche permet de décider du maintien ou du remplacement des systèmes selon des critères technico-économiques. La diminution du rendement de concentration au cours du temps est étroitement liée au vieillissement des matériaux constitutifs du concentrateur solaire (verre, couches minces etc..). Ce dernier est généralement évalué en laboratoire en appliquant des protocoles de vieillissement accéléré normalisés simulant des périodes d'exploitation de l'ordre de 20 à 30 ans. Cette approche ne permet pas une évaluation du vieillissement en exploitation par des moyens non destructifs puisque nécessitant un échantillonnage. Le vieillissement naturel des matériaux entraînent des altérations notamment en surface, modifiant ainsi leur émissivité spectrale. Etant donné que l'émissivité apparente est reliée à l'émissivité spectrale, cet indicateur doit être capable de détecter les altérations de surface. La quantification de l'émissivité apparente met en oeuvre une caméra thermique commerciale. D'un usage simple et rapide, cet outil permet de cartographier la surface étudiée en émissivité apparente puis de convertir ces valeurs en variation de propriétés optiques dans le visible (absorbance, réflectance ou transmittance) suivant la fonction attendue.

Emissivité apparente et émissivité spectrale

L'émissivité apparente intervient dans l'équation radiométrique qu'utilisent implicitement tous les imageurs thermiques : ce paramètre mesure le thermosignal d'émission de l'objet rapporté à celui d'un corps noir à la même température. Nous rappelons l'expression de l'équation radiométrique, mise en oeuvre dans les systèmes infrarouges :

) ( ) ) ( 1 ( ) ( ) ( ) ( . 0 0 . env app t L t t L t t L                (1)
et sont respectivement les thermosignaux de corps noirs associés aux températures d'émission et d'environnement

. Le thermosignal fourni par la caméra est associé à la température apparente Ces thermosignaux sont rapportés à la bande spectrale de détection de la caméra. Une fois la validité du modèle de l'émetteur admise, l'expression de l'émissivité apparente devient :

) ( ) ( ) ( ) ( ) ( . 0 0 . 0 . env env app t L t L t L t L t               (2)
Son mesurage est alors possible à partir des trois températures (apparente, d'environnement, émission) et de la courbe d'étalonnage de la caméra IR. La caractérisation des systèmes est détaillée dans la référence [START_REF] Riou | A self-method for resolving the problem of apparent LWIR emissivity for quantitative thermography up to 130[END_REF]. Cet indicateur est spécifique à la technique de détection mise en oeuvre dans la caméra, de la bande spectrale de détection et de l'émissivité spectrale de la cible.

L'émissivité apparente est liée aux propriétés spectrales de l'émetteur. Chrzanowski [START_REF] Chrzanowski | Problem of determination of effective emissivity of some materials in MIR range[END_REF] propose une formulation générique permettant de raccorder les deux émissivités :

    2 1 2 1 ) ( ) , ( ) ( ) , ( ) , ( ) ( 0 0               d R T L d R T L T T (3)
ε(λ,T) désigne l'émissivité spectrale de la cible, L 0 (λ,T) est la luminance du corps noir et R(λ) est la réponse spectrale normalisée du système. Le produit de ces quantités est enfin intégré sur la bande spectrale [λ 1 -λ 2 ] du système.

Une quantification indirecte de l'émissivité apparente est donc possible moyennant la connaissance de chaque terme de la relation [START_REF] Riou | Accurate methods for single-band apparent emissivity measurement of opaque materials[END_REF]. La connexion de l'émissivité apparente et de l'émissivité spectrale obtenue par réflectométrie IR a été vérifiée expérimentalement pour différents matériaux (acier, verre, polycarbonates) dans la référence [START_REF] Riou | Accurate methods for single-band apparent emissivity measurement of opaque materials[END_REF]. De nombreux mécanismes de vieillissement naturel induisent une modification de l'émissivité spectrale suite à des mécanismes de corrosion, de recristallisation ou d'altération chimique. Les effets du vieillissement naturel peuvent être quantifiés par émissivité apparente si les dégradations produisent des variations d'émissivité spectrale dans la bande spectrale du système infrarouge utilisé. Pour illustrer les capacités de détection de telles variations, nous traitons un corps gris (émissivité spectrale de 0,5) présentant des variations d'émissivité spectrale supposées résulter de modifications structurelles de la surface produites pendant la durée de vie du système. Une forme gaussienne est utilisée pour reproduire de telles variations. Elle est infiniment différentiable et elle permet de simuler le passage d'une variation de type bande Le barycentre de l'anomalie est calculé par: [START_REF] Corwin | Temperature error in radiation thermometry caused by emissivity and reflectance measurement error[END_REF] La position de par rapport au milieu de la fenêtre de détection ( ) détermine la dépendance en température de l'émissivité apparente [START_REF] Ibos | Infrared emissivity measurement device : principle and applications[END_REF]: a. aussi longtemps que λ c < λ m , l'émissivité apparente croît avec la température; b. il est stable avec la température pour λ c =λ m et c. elle diminue avec la température lorsque λ c > λ m . Comme le barycentre de l'anomalie change au cours du processus de vieillissement accéléré, la dépendance en température de l'émissivité apparente peut varier d'un échantillon à l'autre.

               2 1 d ) ( c . rect

Caractérisation des miroirs en émissivité apparente

Nous utilisons une caméra ThermaCAM E300 pour quantifier l'émissivité apparente des miroirs. Elle met en oeuvre une matrice VOx comportant 320 × 240 détecteurs. Sa bande spectrale est [7,9 -12,7μm] [START_REF] Riou | A self-method for resolving the problem of apparent LWIR emissivity for quantitative thermography up to 130[END_REF]. L'incertitude de mesure de température apparente est de 0,1°C. L'écart d'émissivité apparente rend compte du vieillissement accéléré. Les trois miroirs présentent des valeurs quasi-constantes pour des températures d'émission comprise entre 40°C et 110°C, ce qui signifie que la dépendance en température de l'émissivité apparente est conservée pendant le processus de vieillissement accéléré. D'autre part, le processus de vieillissement accéléré est de nature différente selon l'âge initial des miroirs. Il est d'autant plus prononcé que le miroir est initialement jeune. La résolution globale δε Δλ de mesure de l'émissivité apparente est la somme quadratique des résolutions associées à chaque incertitude de température. La répétabilité dépend à la fois de la qualité de la correction de non uniformité (typiquement 0,1°C à 35°C) et de la répétabilité de la température d'émission. L'incertitude combinée dépend de la qualité de la compensation de la dérive thermique. Sans précaution particulière, l'incertitude de la compensation chez la caméra ThermaCAM E300 est typiquement de 0,3°C.

Discussion

Est-il possible d'obtenir les mêmes performances en conditions extérieures? Dans une approche non-intrusive, nous ne pouvons pas envisager une régulation de la température d'émission. Sa valeur est imposée par les conditions ambiantes. Les trois paramètres clés pour atteindre des performances similaires sont a. le contraste t app -t env ; b. l'incertitude de mesure de la température apparente δt app et c. l'uniformité de la température d'émission δt. L'utilisation en extérieur d'un système infrarouge standard avec δt app = 0,1°C et une nonuniformité δt = 0,5°C de la température d'émission (20°C) permettrait d'atteindre une résolution de 1,7% en exposant la cible à un environnement type ciel clair (t env comprise entre -30 et -20°C). La recherche de meilleures conditions de contrôle (δt = 0,25°C) améliorerait la résolution (0,9%). On aurait également avantage à thermographier des miroirs « chaud » i.e. t = 40°C ce qui permettrait d'atteindre une résolution de 0,6% soit 6 ans d'exposition standard. Une étude poussée des performances théoriques de la méthode a été réalisée [START_REF] Riou | Accurate methods for single-band apparent emissivity measurement of opaque materials[END_REF].

Conclusion et perspectives

L'émissivité apparente LWIR est un indicateur pertinent pour quantifier le vieillissement du verre. Il présente des variations de l'ordre de 1,8% au cours du vieillissement accéléré. En comparant l'écart d'émissivité apparente avec l'incertitude combinée de mesure, une approche de laboratoire est capable de discriminer sans contact un vieillissement absolu avec une résolution de 8 années standard d'exploitation. Cette approche peut également distinguer des modifications équivalentes à 3 ans standards sur le même échantillon. L'approche peut être utilisée en situation extérieure. L'utilisation d'un système infrarouge standard (δt app = 0,1°C) associée à non-uniformité δt = 0,25°C de la température d'émission (20°C), permettrait d'atteindre une résolution de 0,9% (10 ans d'exploitation standard) en exposant la cible à un environnement type ciel clair (t env comprise entre -30 et -20°C).

La mesure d'émissivité apparente est sans contact et bénéficie des techniques d'imagerie développées pour la thermographie infrarouge. La cartographie de l'émissivité apparente est possible par l'élaboration d'un code de post-traitement spécifique traitant les thermosignaux. Enfin, il est inutile de développer un système IR spécifique dans la mesure où la caméra commerciale répond aux exigences de répétabilité (qualité de de correction de nonuniformité) et d'incertitude de combinée (qualité de la compensation de dérive thermique).

Un point important est la compréhension du vieillissement accéléré. Il implique des modifications chimiques dues à l'hydratation forcée du verre. Les effets sur la réflectivité des miroirs sont faibles. Ils sont significatifs sur de grandes périodes de temps généralement supérieures à 20 ans d'exploitation. D'autres processus de vieillissement peuvent être quantifiés : c'est le cas de l'érosion. L'érosion produit une fissuration et un écaillage de la surface, rendant les surfaces poreuses au voisinage de l'impact des aérosols solides. Les dommages vus sur le verre après un test de sandstorm [START_REF] Kolli | HF etching effect on sandblasted soda-lime glass properties[END_REF] sont similaires à ceux observés sur la céramique d'alumine de porosité croissante [START_REF] Rozenbaum | Texture and Porosity Effects on the Thermal Radiative Behavior of Alumina Ceramics[END_REF]. La porosité est connue pour modifier les propriétés émissives de la matière, en particulier dans la bande LWIR [START_REF] Rozenbaum | Texture and Porosity Effects on the Thermal Radiative Behavior of Alumina Ceramics[END_REF]. Par exemple, l'émissivité spectrale d'une céramique d'alumine de porosité croissante montre des variations de l'ordre de 80 µm% dans la bande LWIR et devrait produire une variation d'émissivité apparente de l'ordre de 18%. Si un seuil de détection de 1% est considéré, elle permettrait de quantifier les effets précurseurs dus à l'érosion équivalent à 5% de la surface du réflecteur en verre correspondant à 3 années d'érosion naturelle dans les sites océaniques / désert au Maroc [START_REF] Corwin | Temperature error in radiation thermometry caused by emissivity and reflectance measurement error[END_REF].

  étroite à une anomalie de type large bande en faisant varier son écart-type . La largeur de l'anomalie est définie par tandis que son amplitude est quantifiée par (figures 1a, b et c). Le produit définit la surface typique de l'anomalie qui diffère de la surface gaussienne dans un rapport . L'émissivité apparente est calculée en utilisant la relation (3), dans laquelle la réponse spectrale est connue sur la bande normalisée LWIR [8 µm -13μm]. Le résultat est montré sur la figure 1 pour des températures d'émetteur comprises entre -20°C et 180°C. Les effets sont significatifs si nous sommes capables de quantifier des variations d'émissivité apparente de l'ordre du pourcent.

Figure 1 :

 1 Figure 1 : Connexion de l'émissivité apparente avec l'émissivité spectrale des surfaces

Figure 2 :

 2 Figure 2 : Emissivité spectrale hémisphérique initiale des miroirs verre Odeillo

Figure 3 :

 3 Figure 3 : Emissivité spectrale hémisphérique d'Odeillo #1 1 à l'état initial et après 1000h de vieillissement en chaleur humide

Figure 4 :

 4 Figure 4 : Ecart d'émissivité spectrale hémisphérique des miroirs verre Odeillo. a. Odeillo #1; b. Odeillo #2 et c. Odeillo #3. En trait plein: écart [%]; ligne horizontale pointillée: valeur moyenne % calculée sur la bande spectrale [7,9 -12,7 µm] de la caméra ThermaCAM E300; ligne pointillée verticale: barycentre λ C de l'anomalie.

Figure 5 :

 5 Figure 5 : Emissivité apparente LWIR des miroirs Odeillo dans leur état initial.

Figure 6 :

 6 Figure 6 : Emissivité apparente d'Odeillo #1 avant et après 1000h de vieillissement accéléré La caractérisation de l'émissivité apparente est effectuée dans la plage [20°C -120°C] au moyen d'un banc de caractérisation détaillé dans la référence [3]. La figure 5 montre l'émissivité apparente des trois échantillons de verre dans leur état initial. A partir d'une température d'émetteur supérieure à 70°C, la discrimination fondée sur l'âge des échantillons est envisageable. Les valeurs d'émissivité apparentes augmentent de façon monotone avec l'état de vieillissement initial des miroirs. Odeillo #3 est clairement discriminé: l'écart

Figure 7 :

 7 Figure 7 : Ecart relatif % en émissivité apparente LWIR des miroirs verre Odeillo: a. Odeillo #1; b. Odeillo #2 et c. Odeillo #3

LTableau 1 :

 1 'émissivité apparente LWIR est un indicateur pertinent pour quantifier le vieillissement des miroirs verre des héliostats. Les limitations de l'approche résident dans ses capacités de mesurage ainsi que la connaissance des températures d'émission et d'environnement. Les performances typiques de mesures obtenues en situation intérieure sont détaillées dans le tableau 1. Incertitude de type A et B des mesures de températures pour la quantification de l'émissivité apparente LWIR en situation intérieure.


domain: scco scco

text:
Quality of life (QoL) is a predictor of life expectancy and institutionalization among the elderly (Petit & al., 2014). Social relationships play an important role (Gabriel & Bowling, 2014), but they also convey prejudices and stereotypes (Adam & al., 2013) that influence the perceptions attributed to one's own experience of ageing and attitudes towards ageing in its 3 dimensions (physical change; psychological growth; psychosocial losses) (Laidlow & al., 2007). These ones contribute to the development of identity in relation to body image in its two dimensions of global self-perception and body satisfaction (Mucchielli, 2011), for which seniors attach more importance to the functional dimension compared to appearance (Reboussin & al., 2000). This raises the question of the effects of a mind-body practice such as yoga on the quality of life of senior citizens through attitudes to ageing and body image.

Yoga & ageing : a study on attitudes to ageing, body image and quality of life among seniors over 65.

Nathalie Guillaume*, Carole Lefèvre** * Université Paris 8, nathalie.guillaume@icloud.com; ** Université Paris 8-IED, ER IPC Denfert-Rochereau, carole.lefevre@iedparis8.net

Conclusion

Yoga could stand out as an activity that improves quality of life, but it would seem that it is involvement in a group activity that is beneficial. These results echo cognitivist theories which suggest that whatever the activity, a reorganisation of cognitive resources contributes to a better QoL and body image. This study suggests that a mind-body approach could improve QoL by giving older people back an empowerment & helping mitigate the negative effects of self-objectification.

  Adam, S., Joubert, S., & Missotten, P. (2013). L'âgisme et le jeunisme : conséquences trop méconnues par les cliniciens et chercheurs ! Revue de neuropsychologie, 5(1), 4-8. Evers, L., & Verbanck, P. (2010). Création d'un questionnaire de satisfaction corporelle et de perception globale de soi: le QSCPGS. Recherche de normes et début de validation. L'Encéphale, 36(1), 21-27. Gabriel, Z., & Bowling, A. (2004). Quality of life from the perspectives of older people. Ageing & Society, 24(5), 675-691. doi:10.1017/S0144686X03001582 Laidlaw, K., Power, MJ et Schmidt, S. (2007). The attitudes to ageing questionnaire (AAQ): development and psychometric properties. International journal of geriatric psychiatry, 22 (4), 367-379. Marquet, M. et al (2016). A Validation of the French Version of the Attitudes to Aging Questionnaire (AAQ): Factor Structure, Reliability and Validity. Psychologica Belgica, 56 (2), pp. 80-100, DOI: http://dx.doi.org/10.5334/pb.301 Mucchielli, A. (2011). L'identité. Presses Universitaires de France. https://doi.org/10.3917/puf.mucch.2011.01 Petit, S., Bergua, V., Peres, K., Bouisson, J., & Koleck, M. (2014). Élaboration et validation d'une échelle de qualité de vie adaptée aux personnes âgées (EQVPA). Gériatrie et Psychologie Neuropsychiatrie du Vieillissement, 12(4), 379-386. Reboussin, B. A., Rejeski, W. J., Martin, K. A., Callahan, K., Dunn, A. L., King, A. C., & Sallis, J. F. (2000). Correlates of satisfaction with body function and body appearance in middle-and older aged adults: The Activity Counseling Trial (ACT). Psychology and Health, 15(2), 239-254.

  : No mediation of the attitudes to ageing and their dimensions in the link between activity and QoL (p > 0.05) H3 not validated : No mediation of the body image and its dimensions in the link between activity and QoL not validated : no significance in favour of the yoga or recreational group for a better QoL (p > 0.05) H1-b : no significance in favour a QoL T1 > T0 for recreational group (p > 0.05) H1-b validated : Trend towards significance for QoL T1 > T0 for the yoga group. Journées d'Étude du Vieillissement Cognitif TOURS JUIN 2023


text:
Delphine Machot

Carole Lefèvre

Life satisfaction and emotional regulation in aging : temporal perspective and effect of gender

Introduction: Paradoxically, older people report higher levels of life satisfaction than adults under 65 (Urry & Gross, 2010). For some authors, this finding can be explained by the emotional regulation strategies adopted by seniors. Gross and John (2003) postulate that it is the greater use of cognitive reappraisal (though more cognitively costly) and the lesser recourse to expressive suppression that enable access to greater life satisfaction. At the same time, motivational theory postulates that an increasingly narrow temporal perspective would induce the use of more efficient strategies (Carstensen et al., 1999). However, the roles of gender and perception of future time remain to be clarified. Openness: Understanding the specifics of emotional regulation remains a difficult undertaking. However, effective management of negative moods in old age remains a priority, and would seem to benefit from the use of innovative technologies (Virtual Reality, Therapeutic Trainer, Magic Table).

Issues in this study

Limitations : This study was restricted to psychosocial models and selfassessments to examine regulation strategies.

Positive correlation Negative correlation

Negative correlation

No correlation

  : -Does gender influence life satisfaction and emotional regulation strategies in the elderly? -Do gender and perception of future time influence the emotional regulation of older people, and thus their life satisfaction? With Life» Diener et al.(1985). * «Future Time Perspective» (FTP) Castersen & Lang, (1996) Randomized order, individual interviews. Bibliography: • Carstensen, L. L., Isaacowitz, D. M., & Charles, S. T. (1999). Taking time seriously: a theory of socioemotional selectivity. American psychologist, 54(3), 165. • Gross, J. J., & John, O. P. (2003). Individual differences in two emotion regulation processes: implications for affect, relationships, and well-being. Journal of personality and social psychology, 85(2Gender and temporal perspective are not inter-individual differences to be considered in order to better understand the role of emotional regulation in the age paradox. It is the use of compensatory strategies that enables us to continue to be satisfied with our lives despite advancing age.


Jean-Paul Chehab
email: paul.chehab@u-picardie.fr

The Symbol Approximation Method: a numerical approach to the approximation of the symbol of self-adjoint operators June 24, 2023

Keywords: Operator Symbol, Trace of a matrix estimate, inverse problems, Preconditioning, Hydrodynamics damping model AMS Classification[2023]: 65C05, 65F15, 65M55, 65N06, 65N30, 65N55

We propose a simple numerical procedure to approach the symbol of a self-adjoint linear operator A by using trace estimates with numerical data. The Symbol Approximation Method (SAM) is based on an adaptation of the matrix trace estimator to successive distinct numerical spectral bands in order to build a piece-wise constant function as an approximation of the symbol σ of A. The decomposition of the spectral interval into band of frequencies is proposed with several approaches, from the formal spectral to the multi-grid one. We apply the new method to different operators when discretized in finite differences or in finite elements. The SAM is also proposed as a tool for the modeling of waves equations, and is presented a means to capture an additional linear damping term in hydrodynamics models such as Korteweig-de Vries or Benjamin Ono equations.

Introduction

The knowledge of the the symbol of an operator is an important and practical question which arises in different applications, e.g., when finding a missing (linear) term in a PDE to take into account a phenomena and then to enrich the model ; it can be also related to preconditioning techniques for differential operators (find a simple operator to approach a given one in order to build a simple preconditioner).

Indeed, in a number of situations, a physical model has to be corrected (or enriched) by taking into account phenomena that have been ignored at first. For instance, in hydrodynamics, the damping of waves is observed experimentally, and a still challenging question is its mathematical representation, see, e.g. [START_REF] Chehab | Numerical Study of a family of dissipative KdV equations[END_REF][START_REF] Chehab | On Damping Rates of dissipative KdV equations[END_REF][START_REF] Dumont | Numerical investigation of asymptotical properties of solutions to models for water waves with non local viscosity[END_REF][START_REF] Dumont | Numerical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach[END_REF][START_REF] Dutykh | Viscous-potential free-surface flows and long wave modeling[END_REF][START_REF] Ott | Nonlinear theory of ion acoustic wave with Landau damping[END_REF][START_REF] Ott | Damping of solitary waves[END_REF]. The most common approach is to try to express the damping force with a linear monotone self-adjoint operator A; in particular cases it can be obtained by a formal derivation, but it can give rise to tricky expressions difficult to be handled and implemented numercially [START_REF] Dutykh | Derivation of a Viscous Serre-Green-Naghdi Equation : An Impasse ?[END_REF][START_REF] Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF]. Another approach on which we concentrate here consists in approaching A by using experimental or numerical samples.

However, approximating the eigenvalues of A is not sufficient to approach the symbol σ of A since σ relies the frequency number to the eigenvalues. Also σ is intrinsic to the operator. Practically, we propose to extract information on the symbol of the operator through sequences of linear systems

AW (k) = R (k) , k = 1, , N,
where A is a n × n discretization matrix of A; W (k) and R (k) are given; (λ k A ) n k=1 and (σ A (k) k k=1 ) denote respectively the eigenvalues and the symbol of A. Computing the coefficients of A is of course too costly, it is also to strongly related to the discretization scheme, and notably its dimension. We propose to approach σ by a piece-wise constant function on successive band of frequencies. This necessitates to produce an approximation of a mean value of the symbol in a given band of frequencies.To this end we apply trace estimates techniques based on sample vectors and scalar product (see also [START_REF] Ubaru | Applications of trace estimation techniques[END_REF]):

T race(A) = N k=1 λ k A ≈ n N k=1 < AW (k) , W (k) > N k=1
< W (k) , W (k) > , giving directly the mean value µ of the eigenvalues as

µ = N k=1 λ k A n = N k=1 σ A (k) n ≈ N k=1 < AW (k) , W (k) > N k=1 < W (k) , W (k) > .
The approximation of the symbol on band of frequencies is proposed considering the decomposition of the numerical frequencies interval I = [1, n] of A into m sub-intervals as

I = m-1 j=0 I j = m-1 j=0 [k j + 1, k j+1 ].
If we have at our disposal linear filters Π j which allow to extract the eigen-components of a vector associated to eigenvalues belonging in I j , then we can consider the local mean value estimates

µ j ≈ N k=1 < AΠ j W (k) , Π j W (k) > N k=1 < Π j W (k) , Π j W (k) > , j = 1, • • • , m.
The number µ j is then proposed as a estimate of the mean value of eigenvalues of A in I j , and we can build the piecewise constant function

σ m (x) = m j=1 µ j χ I j (x),
as and approximation of k → σ(k), where χ I j (x) denotes the characteristic function of the interval I j , j = 1, • • • m. The separation of band of frequencies (exact or approached) can be done with different techniques: spectral methods, multilevel methods (including multi-grid and hierarchical methods).

The article is organized as follows. In Section 2 we first recall basic and useful results on the estimate of the trace of a matrix using samples that we relay to a mean value estimate of the eigenvalues. Then, we can the define Symbol Approximation Method by first considering the spectral decomposition, applying the mean value estimates of to eigenvalues to successive band of frequencies. In Section 3 we adapt SAM to finite differences or finite elements discretization by using sequences of discretization (such as multi-grid techniques) to produce numerical filter band of frequency. We give illustrations of the multi-grid SAM (MSAM) when applied to different self-adjoint operators, their symbol is captured in a satisfactory way. In Section 4 we consider applications of multilevel SAM to time dependent problems. We first present the way the multilevel SAM can be applied to approach the symbol of a damping differential operator. We apply the method on damped dispersive models as Korteweig-de Vries (KdV) and Benjamin Ono (BO) equations. The numerical computations have been done using Matlab R and Freefem++ [25].

2 Numerical Approximation of the symbol: the spectral case 2.1 Approximation of a matrix and of its trace using samples

Approximation of matrix coefficients

Let A be a matrix in M n (R). We would like to compute the coefficients of A using a sequence of linear systems

AW (k) = R (k) , k = 1, • • • , N, (1)
where (W (k) ) N k=1 and (R (k) ) N k=1 are two given sequences of vectors of R n , related by [START_REF] Abdelouhab | Nonlocal models for nonlinear dispersive waves[END_REF]. This problem is related to the the least square problem min

A∈Mn(R) N k=1 AW (k) -R (k) 2 .
(

The functional Φ(A) = N k=1 AW (k) -R (k) 2 is convex and coercive, the minimum values are reached at the critical points satisfying

N k=1 (AW (k) -R (k) )W T k = 0, or A N k=1 (W (k) (W (k) ) T ) = N k=1 (R (k) (W (k) ) T ), which is equivalent to (2) iff W = [W (1) W (2) • • • W (N )
] is full rank.

We begin with the following simple but useful result:

Proposition 2.1 Let (W (k)
) N k=1 and (R (k) ) N k=1 be two given sequence of vectors of R n There exists a unique matrix A ∈ M n (R) such that

AW (k) = R (k) , k = 1, • • • , N if and only if the matrix [W (1) W (2) • • • W (N )
]. is full rank. In such a case, A is given by A = T S -1 where

S = N k=1 W (k) (W (k) ) T and T = N k=1 R (k) (W (k) ) T Proof. The proof is classical.
Of course, as examples of sequence W k we can consider the case S = W W T = γId, where γ ∈ R * + .This occurs trivially when W K = e k , the k-th element of the canonical basis of R n , it is obviously too costly in practice to consider the whole collection of vector e k . To reduce the computational cost of the problem, it is natural to consider the approximations (A j ) j≥1 of A using the truncated sum

A j S j = T j , j = 1, • • • , N, where S j = [W (1) W (2) • • • W (j) ][W (1) W (2) • • • W (j) ] T and T j = [R (1) R (2) • • • R (j) ][W (1) W (2) • • • W (j) ] T .
Note that W (k) can be a rectangular matrix n × s matrix.

Another strategy consists in choosing (W (k) ) N k=1 as a Bernoulli test vector (or matrix) or a Hadamard matrix (for which W (k) (W (k) ) T = nId). We recall the following two choices of samples of interest:

• Bernouilli test vector or matrices (also as i.i.d Rademacher random variables). They are defined by (P r((W

(k) i = ±1) = 1/2).
• Hadamard Matrices. Following the Sylvester construction, we start from the partitioned matrix

H H H -H .
Starting from H 1 = [1] it generates the sequences of 2n × 2n matrices In Figure [START_REF] Abdelouhab | Nonlocal models for nonlinear dispersive waves[END_REF], we give an illustration of the approximation of the coefficient of the laplacian matrix when using both Hadamard and Bernouilli samples; we observe that the convergence is slow in both the cases.

H 2 k = H 2 k-1 H 2 k-1 H 2 k-1 -H 2 k-1 = H 2 ⊗ H 2 k-1 , k ≥ 1.

Approximation of the trace of a matrix using samples

Let A ∈ M n (R). Consider the problem min D∈Mn(R),D Diagonal A -D 2 F , (3)
which solution is D = diag(A). Here . F denotes the Frobenius norm. We consider the simplest case when D = αId so we look to

min α∈R A -αId 2 F . (4)
This last problem solves easily: we set φ(α) = A -αId 2 F , we find that the minimum is reached for

α = < A, Id > F < Id, Id > F = trace(A) n = n k=1 λ k n .
If A is only known trough a sequence of matrix-vector products (1), α will be approached numerically such as minimizing

N k=1 AW (k) -αW (k) 2 .
We find

α opt = N k=1 < AW (k) , W (k) > N k=1 < W (k) , W (k) >
. We then obtain the following approximation of the trace

T race(A) n N k=1 < AW (k) , W (k) > N k=1 < W (k) , W (k) > .
From the previous relation, taking (W (k) ) i = ±1 with equiprobability, 1 2 , we find

T race(A) 1 N N k=1 < AW (k) , W (k) >,
and we recover the method given proposed in [START_REF] Hutchinson | A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines[END_REF], then studied and applied, e.g. in [START_REF] Avron | Ramdorized algorithms for Estimating the Trace of an implicit symmetric Positive Semi-Definite Matrix[END_REF][START_REF] Bai | Some large-scale matrix computation problems[END_REF][START_REF] Ubaru | Applications of trace estimation techniques[END_REF]; more generally, let (W (k) ) N k=1 can be a sequence matrices of size n × s and we have

T race(A) ≈ 1 N s N k=1 < (AW (k) ), W (k) > F . (5)
We now recall the result from [START_REF] Avron | Ramdorized algorithms for Estimating the Trace of an implicit symmetric Positive Semi-Definite Matrix[END_REF][START_REF] Hutchinson | A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines[END_REF] which gives a justification of the above formula:

Lemma 2.2 Let A be an n × n symmetric matrix with trace(A) = 0. Let z be a random vector whose entries are i.i.d Rademacher random variables (P r(zi = ±1) = 1/2). z T Az is an unbiased estimator of trace(A) i.e., E(zT Az) = trace(A), and

V ar(z T Az) = 2 A 2 F - n i=1 A 2 ii .
As an illustration, we display in Figure [START_REF] Avron | Ramdorized algorithms for Estimating the Trace of an implicit symmetric Positive Semi-Definite Matrix[END_REF], the approximation of the trace of the laplacian matrix when using both Hadamard and Bernouilli samples ; we observe that the the convergence is fast in both the cases. Relativ error We make now a simple but important remark regarding the applications we look to: when the matrix A is known only through matrix-vector products AW k = R k the trace estimation (5) writes as

Computation of trace(A) using Hadamard and Bernouilli samples

T race(A) N N k=1 < R (k) , W (k) > m k=1 < W (k) , W (k) > . ( 6
)
and particularly when (W (k) ) i = ±1 with equiprobability, 1 2 , we obtain

T race(A) 1 N N k=1 < R (k) , W (k) >.

Estimates of spectral mean values on band of frequencies

Let A be a symmetric matrix of M n (R). We introduce the following notation:

• σ A (k) is the value of the numerical symbol of the matrix A for the frequency number A

• σ is the symbol of the operator A.

Before presenting the method, we make the following remark: the eigenvalues λ k A of A coincide with the values of σ A , however, the numbers σ A (k), k = 1, • • • , n are numbering in the increasing order of the frequency number k. Hence

n k=1 λ k A = n k=1 σ A (k).
We first consider the so-called spectral approach. Let A be a symmetric matrix. A being diagonalizable. We decompose the numerical frequency interval I = [1, n] as

I = ∪ m-1 j=0 I j ,
where

I j = [k j + 1, k j+1 ], k 0 = 0, k m = n, k j < k j+1 .
Let (w j ) n j=1 be the eigenvectors of A defined by the relations Aw j = σ A (j)w j , j = 1, • • • , n. We insist on the fact that

{λ k A , k = 1, • • • , n} = {σ A (k) k = 1, • • • , n}
, however the eigenvalues of A are not necessarily numbered in the increasing order of the numerical frequencies. We denote by

L j = span(w k j +1 ..w k j+1 )
and by P j the (k j+1 -k j ) × n-matrix P j = [w k j +1 ..w k j+1 ] T . The orthogonal projector onto L j is then

Π j = P T j P j . For any u ∈ R n , we can write u = n i=1 ûi w i , so Π j (u) = k j+1 i=k j +1
ûi w i and, since A and Π j commute, and Π 2 j = Π j , j = 1, • • • , m, we have

< AΠ j u, Π j u > < Π j u, Π j u > = < AP j u, P j u > < P j u, P j u > = k j+1 i=k j +1 û2 i σ A (i) k j+1 i=k j +1 û2 i . Now, for a sequence of vectors u ( ) ∈ R n , = 1, • • • , N , we define µ j = N =1
< AP j u ( ) , P j u ( ) > N =1 < P j u ( ) , P j u ( ) > .

We have then

µ j = k j+1 i=k j +1 N =1 ( û
( ) i ) 2 σ A (i) k j+1 i=k j +1 N =1 ( û
( ) i ) 2
, so

µ j = k j+1 i=k j +1 γ i σ A (i), with γ i = N =1 ( û
( ) i ) 2 k j+1 i=k j +1 N =1 ( û
( ) i ) 2 ∈ [0, 1], k j+1 i=k j +1 γ i = 1.
R j is then a convex combination of λ i and is proposed as an estimator of

k j+1 =k j +1 σ A k j+1 -k j
, the mean value of the partial symbol in the numerical frequency number interval [k j + 1, k j+1 ]. Indeed, we have

µ j = N =1 < (P T j AP j )u ( ) , u ( ) > N =1 < P T j P j u ( ) , u ( ) > .
When u is a random vector whose entries are i.i.d Rademacher random variables, we have

µ j = N =1 < (P T j AP j )u ( ) , u ( ) > N (k j+1 -k j ) .
At this point, we observe that, when u is a random vector whose entries are i.i.d Rademacher random variables, we obtain, applying Lemma 2.2,

E(u T A j u) = trace(P T j AP j ) = k j+1 =k j +1 σ A ( ), j = 1, • • • , m,
where

A j = P T j AP j , j = 1, • • • , m. Therefore µ j = k j+1 =k j +1 σ A ( ) k j+1 -k j .
Remark 2.3 This approach applies to the approximation of matrices expressed as B -1 A, which is related to generalized eigenvalues problems Au = λBu. This situation arises when considering finite elements discretization, as we will see later on. The ration µ j is then replaced by

µ j = N =1 < (P T j AP j )u ( ) , u ( ) > N =1 < P j Bu ( ) , P T j u ( ) >
As pointed out above, when the matrix A is known only through matrix-vector products A(P j u ( ) ) = v ( ) j the trace estimation writes as

µ j = m =1 < v ( ) , (P j u ( ) ) > m =1 < (P j u ( ) ), B(P j u ( ) ) > . ( 7
)
We resume as follows the spectral Symbol Approximation Method (SAM) to approach the spectrum of the generalized eigenvalue problem Algorithm 1 : Spectral SAM 1: Given the projectors

P j , j = 1, • • • , m and a sequence of samples (u ( ) ) N =1 and (v ( ) ) N =1 ∈ R n satisfying the relations Bu ( ) = v ( ) , 2: for = 1, • • • , N do 3: for j = 1, • • • m, do 4: Compute v ( ) j = A(P j u ( ) ) 5: Set N j = N j + < v ( ) , (P j u ( ) ) > 6:
Set D j = D j + < (P j u ( ) ), (P j u ( ) ) >

7:

end for 8: end for 9:

for j = 1, • • • m, do 10: Set µ j = N j D j 11: end for 12: Set σ m (x) = m j=1 µ j χ I j (x)
At this point we give the following result:

Proposition 2.4 Let A ∈ M n (R) be a symmetric matrix with symbol σ A . Decompose the numerical frequency interval I = [1, n] as I = ∪ m-1 j=0 [k j + 1, k j+1 ] where I j = [k j + 1, k j+1 ], k 0 = 0, k m = n, k j < k j+1 .
Then, the numbers µ j , j = 1, • • • , m defined by Algorithm (1) satisfy

σ A j ≤ µ j ≤ σAj , j = 0, • • • , m.
where σ A j = min p∈I j σ A (p) and σAj = max p∈I j σ A (p)

Proof. Let (w k ) n k=1 be the orthonormal eigenvectors of A, with the numbering

Aw k = σ(k)w k , k = 1, • • • , n. Let now j ∈ 0, 1, • • • , m be fixed. We have Π j u ( ) = k j+1 p=k j +1 û( ) p w p and (AΠ j u ( ) ) = k j+1 p=k j +1 λ p û( ) p w p . Therefore σ A j N =1 k j+1 p=k j +1 (û ( ) p ) 2 ≤ N =1 < A(Π j u ( ) ), Π j u ( ) > = N =1 k j+1 p=k j +1 σ(p)(û ( ) p ) 2 ≤ σAj N =1 k j+1 p=k j +1 (û ( ) p ) 2 .
We have then

σ A j ≤ µ j = N =1 < A(P j u ( ) ), P j u ( ) > N =1 < (P j u ( ) ), (P j u ( ) ) > = k j+1 p=k j +1 σ(p)(û ( ) j,p ) 2 ≤ σAj 2.

Illustrations

Finite differences discretization

We first consider A as the discretization matrix of the laplacian on [0, 1] with periodic boundary conditions. We take n = 256 and a bandwidth m = 4 (the length of each band-with is then 32). In The coefficient of x 2 is 9.48 π 2 , which is rather good and fits very well with the laplacian's symbol leading term . The large values of the other coefficients are attribuable to the size of the bandwidth. The sample vectors are of Bernoulli type.

In space dimension 2, the symbol defines a surface that we fit with a polynomial of the two variables x and y. The computed symbol of the operator is σm (x, y) = 1.0e + 04 (6.702518213669119 -4.960659057862875x -5.612707457168954y) +1.0e + 04 0.672528694583729(

x 2 + y 2 ) ,
which is a correct result allowing to propose the laplacian as the operator A, see Figure (4).

Finite Elements discretization

We represented in Figure [START_REF] Bona | Singularity formation in the generalized Benjamin-Ono equations[END_REF] the approximation of the symbol of A = -∆ in a square domain with holes.

The SAM method allows to capture in a very satisfactory way σ. 3 The Multigrid Symbol Approximation Method (MSAM)

Derivation of the method

In the practical cases, the matrix A will be the discretization of a (local or nonlocal) differential operator which symbol can be assumed as a function defined on the spectral domain of -∆. Hence, the eigenvalues of A are related to the natural frequencies of the domain, says to the eigenvalues λ i of -∆. Computing these numbers can be costly, and a way to select a band of frequencies is to use numerical filters that will replace the projectors P j . A classical procedure to realize numerically a separation of band of frequencies is to use grids or meshes of different characteristic mesh-size: this consists in handling different level of discretization (and associated spaces of approximation), from the coarsest to the finest, each one with a limited capability to capture a set of frequencies (by Nyquist's frequency sampling theorem). The coarse spaces allow to represent only signals supported by low frequencies while the fine ones cover a large band. Consider the solution of the linear system

Ax = b,
Here A corresponds to the discretization matrix of an operator A but its coefficients are not known, A is only known through matrix-vector products; this e.g. is the case when A is a Schur complement. In such a situation the use of Krylov methods is recommended and the numerical resolution is of course more efficient with the use of a pre-conditioner. The problem here is then to build a pre-conditioner of A without using the coefficient of A. A way to overcome this difficulty is then to approach the symbol σ of A by a simple polynomial (or rational function) σm . The pre-conditioner of A is then built as a discretization matrix of the operator of symbol σ m .

To this end, we can of course apply the SAM method as presented above.

We first introduce here a simple variant of SAM when a sequence of matrices A is known, say when we have at our disposal discretization of A with different d.o.f. Let A 1 and A 2 be two discretization matrices of A of size n 1 × n 1 and n 2 × n 2 , respectively, with n 2 > n 1 . On the one hand the trace approximation coming from the consistency of the discretization gives:

trace(A 1 ) n 1 i=1 σ(k), trace(A 2 ) n 1 i=1 σ(k) + n 2 i=n 1 +1 σ(k),
where k → σ(k) is the symbol of A.

On the other hand, the trace estimation of the SAM gives

σ1 = trace(A 1 ) n 1 N =1 < A 1 W ( ) 1 , W ( ) 1 > N i=1 < W ( ) 1 , W ( ) 1 > , σ2 = trace(A 2 ) n 2 N =1 < A 2 W ( ) 2 , W ( ) 2 > N i=1 < W ( ) 2 , W ( ) 2 >
.

Here W ( )

1 ∈ R n 1 and W ( ) 2 ∈ R n 2 , = 1, • • • , N , are two sequences of sample vectors.
Now, applying the trace estimation, we deduce directly from above

µ 1 = σ1 , (8)
µ 2 = 1 n 2 -n 1 (n 2 σ2 -n 1 σ1 ) (9)
More generally let n j , j = 1, • • • , m be integers labelled in the increasing order. We consider the m matrices A j of sizes n j × n j . We denote by λ k , k ≥ 1 the eigenvalues of A and λ

(j) k , k = 1, • • • , n j the ones of A j .
The trace approximation results gives:

Algorithm 2 : MSAM method 1: Compute σj = N =1 < A j W ( ) j , W ( ) j > N i=1 < W ( ) j , W ( ) j > 1 n j trace(A n j ) , j = 1, • • • , m 2: Initialization µ 1 = σ1 3: for j = 2, • • • m, do 4: Set µ j = 1 n j -n j-1
(n j σj -n j-1 σj-1 )

5: end for

6: Set σ m (x) = m j=1 µ j χ [n j-1 +1,n j ] (x)
Here W ( )

j ∈ R n j , = 1, • • • , N
, are the sequences of sample vectors.

Remark 3.1 The expression giving µ j from σj and σj-1 in Algorithm ( 4) is nothing else but a (simple) extrapolation formula; more general extrapolations techniques could be considered to improve the predictive formula, we refer to [START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF][START_REF] Brezinski | Extrapolation Methods. Theory and Practice[END_REF], e.g., for general techniques and surveys.

Proposition 3.2 Let A be a self-adjoint positive definite operator and A n j a sequence of n j × n j discretization SPD matrices of A. We note by σ(k) > 0 the symbol values of A and σ A j (k) those of A n j . We introduce the error e (k,j) = σ(k) -σ A j (k).

We have the relations:

i.

µ j - 1 n j -n j-1 n j k=n j +1 σ(k) = 1 n j -n j-1 n j k=1 e (k,j) - n j-1 k=1 e (k,j-1)
ii.

1 - 1 n p σAp ≤ np j=1 n j -n j-1 n p µ j ≤ σAp , p = 2, • • • , m.
Proof. The proof is established by direct computations.

Applications in spectral Fourier discretization

As a first illustration, we look to the performance of the MSAM method for approaching the symbol of a band-limited filter which symbol is θ(k) = 1 if 50 < k < 100, 0 else, say, a characteristic function in Fourier space. We display the comparison between the initial symbol and its approximation in Figure [START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF].

The result is quite satisfactory, the symbol being difficult to approach; FFT are used to build the matrix vector product with θ in the physical space.

Applications in finite differences

As a first illustration, we compare in Figure [START_REF] Brezinski | Extrapolation Methods. Theory and Practice[END_REF] the capability of the MSAM method to approximate the symbol of the laplacian matrix, when using spatial discretization schemes of order 2, 4 and 6, the two last ones being built with compact schemes. In all the cases, we remark that the approached symbol is very closed to the exact symbol for the low frequencies: this is due to the very good approximation of the first eigenvalues by the discretization scheme, the accordance is even better when considering the high accurate compact schemes (order and 6), as expected by formula i. of Proposition 3.2.

We now consider the operator

A = 100(-∂ 2 x )(1000 -∂ 2 x ) -1 which symbol is σ(k) = 100π 2 k 2 1000 + π 2 k 2
, and a sequence of finite differences discretization matrices. We plot in Figure [START_REF] Cabral | Chaos for a damped and forced KdV equation[END_REF] the approximated spectrum, the fitted curve, and the exact symbol. The nonlinear least-square fitting gave

σ m (k) = 8.2614 + 2.0832.π 2 k 2 2.5260 + 0.0208π 2 k 2
The first part of the symbol, which corresponds to the low frequencies, is not as well represented as in the laplacian case (see preceding remarks on Proposition 3.2). However, the result can be considered as relatively satisfactory since only matrix vector product have been used.

x ) -1 , with 4th order discretization schemes. Nmax=275, N=50, m=6

As a second type of illustration we give application of the MSAM method to approximation of the spectrum of the Schur complement (or Uzawa's discrete operator) of the generalized Stokes problem:

αu - 1 Re ∆u + ∇p = f, in Ω =]0, 1[ 2 (10) divu = 0 in Ω, (11)
u = 0 ∂Ω, (12)
here u = (u, v) is the velocity field and the scaler p is the pression. The Uzawa operator is formally

U = div αu -1 Re ∆ -1
∇. When α = 0, U it is closed to the Re Id but for large values of α, U can be approached by a laplacian with Neumann Boundary conditions (see [START_REF] Peyret | Computational Methods for Fluid Flow[END_REF]), and an the Schur complement of the Uzawa-type method needs to be pre-conditionned in this case. The (second order) discretization of the Stokes system on a MAC mesh with finite differences reads as

αU + 1 Re A u .U + B x P = 0, ( 13
)
αV + 1 Re A v .V + B y P = 0, (14)
B T x U + B T y V = 0, (15)
with natural notations. The Schur Complement is

S = B T x (αId + 1 Re A U ) -1 B x + B T y (αId + 1 Re A V ) -1 B y .
The results of the MSAM method are displayed in Figure ( 9) for α = 0. We observe that expect the lowest eigenvalues close to 0, 83 Re most of the other ones are concentrated around Re = 1000; this agrees with [START_REF] Crouzeix | Approximation et méthodes itératives de résolution d'inéquations variationelles et de problèmes non linéaires, I.R.I.A[END_REF][START_REF] Crouzeix | On an operator related to the convergence of Uzawa?s algorithm for the Stokes equation[END_REF] . Now we take α = 1000, so U -∆. In Figure [START_REF] Calgaro | Séparation des échelles et schémas multiniveaux pour les équations d'ondes non-linéaires[END_REF] we display the computed MSAM approximation of the symbol of S, and we compare it with the one of the discrete Laplacian, the eigenvalues being ranged in the increasing order. Then, we use the well-known pre-conditioner given by Cahouet and Chabard [START_REF] Cahouet | Some fast 3d finite element solvers for the generalized Stokes problem[END_REF] as a symmetric pre-conditioner of S build of (incomplete) Cholesky factors of the negative Laplacian matrix (associated to Homogeneous Diirichlet Boundary conditions), and we observe that it is concentrated around the point (1, 0) in the complex plane.

Applications in Finite Elements

First of all, we must adapt the trace estimation formula. We consider for the sake of simplicity the homogeneous Dirichlet problem on the regular open and bounded domain Ω ⊂ R 2,3 ; the problem is discretized in finite elements on a (regular) triangulation T h on which the finite elements space V h is built. We obtain the linear system

A h u h = M h f h ,
where A h is the stiffness matrix and M h the mass matrix; h is the characteristic mesh-size of the triangulation. The trace estimation formula writes here as

.eps T race(

A h ) n N k=1 < A h W k , W k > N k=1 < M h W k , W k > .
where W k ∈ R n , k = 1, • • • , N are sample vectors representing sample functions of V h . We have set n = dim(V h ). Now, to display the SAM method in a multi-grid like framework, we can consider an embedded sequence of finite element spaces, from the coarsest V h 1 to the finest V hm , with

V h 1 ⊂ V h 2 ⊂ • • • ⊂ V hm ,
but also, in a more general way, increasing dimensional finite elements spaces,

dim(V h 1 ) < dim(V h 2 ) < • • • < dim(V hm ),
not necessarily embedded, and it will be the case for the numerical illustrations we present hereafter. We denote by A h j the associated stiffness and by M h j , the mass matrices j = 1, • • • , m. We set n j = dim(V h j ).

Algorithm 3 : MSAM method in Finite Elements

1: Compute λj N k=1 < A h j W (j) k , W (j) k > N k=1 < M h j W (j) k , W (j) k >
.

2: Initialisation µ 1 = λ1 3: for j = 2, • • • m, do 4: Set µ j = 1 n j -n j-1 n j λj -n j-1 λj-1 5: end for 6: Set σ m (x) = m j=1 µ j χ [n j-1 +1,n j ] (x)
Remark 3.3 It is not necessary in practice to have embedded finite elements spaces to apply the MSAM method, say T h j ⊂ T h j+1 , j = 1, • • • , m -1, as illustrated below.

The reader will find in [START_REF] Kuttler | Eigenvalues of the laplacian in two dimensions[END_REF] a extensive description of the eigenvalues and the eigenfunctions of the laplacian in space dimension two, for a number of domain shapes. Identitykit Spectrum Exact Operator Spectrum

Hierarchical methods

Another way to decompose a signal into a set of components, each one attached to a band of frequencies is the hierarchical approach. The separation of the components (or filtering) is realized considering several levels of discretization, as in the multi-grid methods in finite differences or finite elements. The leading idea is based on a numerical filtering of the solution; this principle can however be applied to other discretization such as wavelets [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], or spectral methods [START_REF] Costa | Time marching techniques for the nonlinear Galerkin method[END_REF]. We recall here briefly hereafter the hierarchy process in finite differences. We first consider for simplicity 2 levels of grid in the 1D case. The use of different levels of discretization allows to decompose a signal (or values of a function at given grid points) into mean part and oscillatory (fluctuent) part. As in multi-grid frameworks we consider two levels of discretization: given h = 1 2n , the coarse grid

G 2h = {x i = 2hi, i = 0, • • • n -1} and the fine grid G h = {x i = ih, i = 0, • • • 2n -1}.
The filtering consists in computing proper local average of the function u at the grid points of G h ⊂ G 2h and to leave unchanged the values of u at the grid points of G 2h ; the resulting vector ū carries low frequencies components of u the high ones are smoothed by the average procedure; the resulting signal z = u -ū is a high modes correction of ū to u. The local average is based on a interpolation scheme and that the separation is more efficient when a high order compact scheme is applied, see [START_REF] Chehab | Stabilization and Numerical Filtering for the Modeling and the Simulation of time dependent PDEs[END_REF][START_REF] Chehab | Incremental Unknowns Method and Compact Schemes[END_REF] and the references therein.

We can resume the process with the two following steps:

• Step1 (Hierarchization): rewrite the vector of unknown as u = (y,

u f ) t ; y ∈ G 2h , u f ∈ G h \ G 2h • Step 2 (linear change of variable) replace u f by z = u f -R(Y ), where R : G 2h → G h \ G 2h
is a p-th order linear interpolation scheme; we have then z = O(h p ) for regular signals, the simplest one being the mid-point interpolation:

z 2i+1 = u 2i+1 -1 2 (u 2i + u 2i+2 ) = O(h 2
). Higher interpolation schemes can be implemented using compact schemes.

As an illustration we represent hereafter by a crux × the grid points belonging to G 2h and those of G h \ G 2h by a circle o , and by . the points at the boundaries:

× o × o × o × o . . . . . . . . . . o o o o o o o . . o × o × o × o . . o o o o o o o . . o × o × o × o . . o o o o o o o . . o × o × o × o . . o o o o o o o . . . . . . . . . .
G 2 d h ), G 2 d h ⊂ • • • G 2h ⊂ G h ,
we can rely the original (nodal) components to the new ones as

     y u f 1 . . . u f d      = S      y z 1 . . . z d     
, with obvious notations. Matrix S is called the transfer matrix. A a general construction in space dimension 2 and 3 is presented in [START_REF] Chehab | Incremental Unknowns Method and Compact Schemes[END_REF].

We give hereafter as an illustration in Figures [START_REF] Chehab | On Damping Rates of dissipative KdV equations[END_REF] and ( 15) the decomposition of approximation of u(x) = sin(100x(1 -x)) at the fine grid points of G h into several levels of oscillatory components, when second and fourth order interpolation formula are used respectively for the change of variable S. We observe that the fine grid corrections are smaller when using a higher order interpolation scheme. To apply the SAM method, we must handle vectors of the same dimensions. We write

u =      y u f 1 . . . u f d      = S           y 0 . . . 0      +      0 z f 1 . . . 0      + • • • +      0 0 . . . z f d           = ũ1 + • • • ũm (16)
We can now present the method as follows:

Algorithm 4 : Hierarchical SAM (HSAM)

1: Initialization: S 1,j = S 2,j = 0, j = 1, • • • , m 2: for k = 1, • • • N do 3: Decompose W k as W k = W (1) k + • • • + W (m)
k following formula [START_REF] Chen | Decay of solutions to a water wave model with nonlocal viscous dispersive term[END_REF] 4:

for j = 1, • • • m do 5: Compute S 1,j = S 1,j + < A W (j) k , W (j) k > 6: Set S 2,j = S 2,j + < W (j) k , W (j) k > 7:
end for 8: end for 9:

for j = 1, • • • m do 10: Compute µ j = S1 j S2 j 11: end for 12: Set σ m (x) = m j=1 µ j χ [n j-1 +1,n j ] (x)
We give hereafter in Figures [START_REF] Chen | Decay of solutions to a water wave model with nonlocal viscous dispersive term[END_REF] and ( 17) the approximation of the spectrum of the negative laplacian matrix in space dimensions 1 and 2 respectively.

Sample mean Spectrum mean data1

Figure 17: Finite differences discretization. HSAM for the 2D Laplacian n = 127 2 , A is n × n-2d order discretization scheme (column 1), 4th order discretization scheme (column 2). Second order Filter (2d order interpolation scheme).

On the influence of the accuracy of the discretization schemes

We here give a description of the error of the SAM method. Let A n ∈ M n (R) be a symmetric matrix, obtained by a consistent discretization of the operator A. We note by σ An the numerical symbol of A n and by σ A the one of A. We consider the decomposition of the numerical frequency interval [1, n] as

[1, n] = ∪ m-1 j=0 [k j + 1, k j+1 ].
For a given sequence of sample vectors (W (k) ) N k=1 ∈ R n , the error of the SAM method satisfies the estimate:

|σ A (k) - N =1 < AW ( ) j , W ( ) j > N i=1 < W ( ) j , W ( ) j > | ≤ |σ An (k) -σ A (k)| + |σ An (k) - N =1 < AW ( ) j , W ( ) j > N i=1 < W ( ) j , W ( ) j > |, k j + 1 ≤ k ≤ k j+1 , ∀j = 1, • • • , m.
This is a simple consequence of the triangular inequality. It allows to decompose the error

E = E 1 + E 2 into E 1 = σ An (k) -σ A (k)
which depends only on the discretization scheme (and its accuracy), and

E 2 = N =1 < AW ( ) j , W ( ) j > N i=1 < W ( ) j , W ( ) j >
-σ An (k), which depends on the "quality" of the sample vectors and of the capability of the numerical filtering to separate group of frequencies. We now can give the following result:

Proposition 3.4 Let A n ∈ M n (
R) be a symmetric matrix, obtain by a consistent discretization of the operator A. We note by σ An the numerical symbol of A n and by σ A the one of A. We consider the decomposition of the numerical frequency interval [1, n]

as [1, n] = ∪ m-1 j=0 [k j +1, k j+1 ].
For a given sequence of sample vectors (W (k) ) N k=1 ∈ R n , the error of the SAM method satisfies the estimate: If we assume in addition that for any fixed M ∈ N * , we have:

lim n→+∞ max 1≤k≤M |σ An (k) -σ A (k)| = 0, then, if n is large enough, |σ A (k) - N =1 < AW ( ) j , W ( ) j > N i=1 < W ( ) j , W ( ) j > | ≤ 2|σ An (k) - N =1 < AW ( ) j , W ( ) j > N i=1 < W ( ) j , W ( ) j > |, ∀k ≤ M.
Proof. The inequality shows that the more accurate is the discretization scheme, the smaller is the error term E 1 so, for n large enough and for a fixed M , we obtain the result.

We now give an illustration of this result by considering the approximation of the symbol of A = -∆ on ]0, 1[ associated to Homogeneous Dirichlet Boundary conditions. We consider the discretization of A by 2d, 4th and 6th order discretization schemes, the two last ones being obtain by using compact scheme, [START_REF] Lele | Compact Difference Schemes with Spectral Like resolution[END_REF], we used N = 50 sample Rademacher vectors. Figure [START_REF] Crouzeix | Approximation et méthodes itératives de résolution d'inéquations variationelles et de problèmes non linéaires, I.R.I.A[END_REF] gives an illustration of Proposition 3.4: we observe that passing to the 4th and 6th order discretization schemes of -∆, MSAM gives a very satisfactory approximation σ. Also, in this context, the multi-grid filtering is effective. Symbol of -

SAM for time dependent problems

In a number of situations, a physical model has to be corrected (or enriched) by taking into account phenomena that have been ignored at first. For instance, in hydrodynamics, the damping of waves is observed experimentally, and a still challenging question is its mathematical representation, see, e.g. [START_REF] Chehab | Numerical Study of a family of dissipative KdV equations[END_REF][START_REF] Chehab | On Damping Rates of dissipative KdV equations[END_REF][START_REF] Dumont | Numerical investigation of asymptotical properties of solutions to models for water waves with non local viscosity[END_REF][START_REF] Dumont | Numerical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach[END_REF][START_REF] Dutykh | Viscous-potential free-surface flows and long wave modeling[END_REF][START_REF] Ott | Nonlinear theory of ion acoustic wave with Landau damping[END_REF][START_REF] Ott | Damping of solitary waves[END_REF]. The most common approach is to try to express the damping force with a linear operator B; in particular cases it can be obtained by a formal derivation, but it can give rise to tricky expressions difficult to be handled and implemented numercially [START_REF] Dutykh | Derivation of a Viscous Serre-Green-Naghdi Equation : An Impasse ?[END_REF][START_REF] Meur | Derivation of a viscous Boussinesq system for surface water waves[END_REF]. Another approach on which we concentrate here consists in approaching B by using experimental or numerical samples.

Derivation of the method

We adapt here the numerical approximation to the symbol developed above. Consider the evolution system

∂u ∂t + Au + F (u) = 0, t ∈)0, T ), (17)
which is valid in absence of damping phenomena. We would like to identify a damping operator by fitting a set of experimental data with the enhanced model

∂u ∂t + Au + F (u) + Bu = 0, (18)
where B represents the (unknown) damping operator; B is assumed to be time independent, auto-adjoint and positive definite. When discretizing in space the last system we get

du dt + Au + F (u) + Bu = 0, t ∈ [0, T ]. ( 19
)
Assume that we have at our disposal measured physical data W (t) at discrete times t k ∈ [0, T ]; It is important to precise that the dimension of the differential system is supposed to be fixed, e.g. limited in practice by the size of exprimental mesures W (t). We would like to compute B in such a way to fit with W . To this end we consider a stable time marching scheme, e.g.

W (k+1) -W (k) ∆t + AW (k+1) + F (W (k+1) ) + BW (k+1) = 0,
or, equivalently,

BW (k+1) = - W (k+1) -W (k) ∆t -AW (k+1) -F (W (k+1) ) = R (k+1) , k = 0, • • • , N -1, (20)
with N ∆t = T . We propose to apply the multilevel SAM as follows:

• The sequence W (k) is supposed to be known, in practice it can be given by experimental measures.

• Compute R (k+1) , k = 0, • • • , N -1 by formula (20),

• Decompose each W (k+1) and R (k+1) using frequency filters Π j , j = 1, • • • , m (spectral, multi-grid or so) into

W (k) j = Π j W (k) and R (k) j = Π j R (k) , j = 1, • • • , m, • Compute S (k) j =< R (k+1) j , W (k) j > and T (k) j =< W (k+1) j , W (k) j >, j = 1, • • • , m,
• Set

µ j = N -1 k=0 S (k) j N -1 k=0 T (k) j , j = 1, • • • , m.

Application to the Korteweg-de Vries equation

We consider the Korteweg-de-Vries model (KdV), which is obtained from Euler's equations by selecting a particular physical regime: small amplitude elevation, large wavelength, unidirectional propagation, see [START_REF] Miranville | Mathematical Modeling in Continuum Mechanics[END_REF]; it addresses then to low frequency regimes. The long time behavior of dissipative asymptotic models is still an important issue: the capture of damping rates in several norms, the measure of regularization effects, the evidence of complex asymptotic dynamics, just to name but a few, are important questions to be considered when trying understand natural phenomena. Mostly, several of these questions are still open and the numerical simulation is a way to capture some properties, to select pertinent models and to develop theoretical as well as practical strategies.

We here focus on KdV equations on the torus T = T(0, L); other dispersive models such as BBM (Bona-Benjamin-Mahony) or BO (Benjamin-Ono) equations could also studied following a similar approach. Damped Korteweg-de Vries equations appear in different physical situations, they can be expressed in a large generality as

u t + B(u) + u xxx + uu x = 0, x ∈ T, t > 0, (21)
where B is a linear operator, defined on a Hilbert space V , subspace of L 2 with values in L 2 , and satisfying

L 0 B(v)vdx ≥ 0, (22)
for all function v ∈ V , regular enough, in such a way the L 2 -norm of the solution is decreasing in time as 1 2

d|u| 2 L 2 dt + L 0 B(u)udx = 0. ( 23
)
We find in the literature different choices for L, depending on the physical situations, see [START_REF] Cabral | Chaos for a damped and forced KdV equation[END_REF][START_REF] Chehab | Numerical Study of a family of dissipative KdV equations[END_REF][START_REF] Chehab | On Damping Rates of dissipative KdV equations[END_REF][START_REF] Chehab | Long-time behavior of solutions of a BBM equation with generalized damping[END_REF][START_REF] Chen | Decay of solutions to a water wave model with nonlocal viscous dispersive term[END_REF][START_REF] Dias | Viscous potentiel free-surface flows in a fluid layer of finite depth[END_REF][START_REF] Dumont | Numerical investigation of asymptotical properties of solutions to models for water waves with non local viscosity[END_REF][START_REF] Dumont | Numerical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach[END_REF][START_REF] Dutykh | Viscous-potential free-surface flows and long wave modeling[END_REF][START_REF] Ott | Nonlinear theory of ion acoustic wave with Landau damping[END_REF][START_REF] Ott | Damping of solitary waves[END_REF] and [START_REF] Jäger | Experiments on KdV Solitons[END_REF] for physical experiments.

A sampling we consider the solution generated by the damped KdV model, starting from a soliton initial datum. The initial datum is

u 0 (x, 0) = α sec 2 (β(x -x0) -C.t), with = 1, L = 100, x 0 = L/2, α = 0.8, C = 2α/6, β = 1 2 √ C .
In Figure [START_REF] Crouzeix | On an operator related to the convergence of Uzawa?s algorithm for the Stokes equation[END_REF], we represent the energy spectrum of u 0 ; we observe that the bandwidth is relatively important. We first display the results when approaching the symbol of the damping operator by local spectral means. The results show a very good agreement and the symbol of the damping operator is very correctly captured. This is also observed for operator like Id, √ -∆.

The results of SAM when using spectral Fourier space discretization are reported in Figure [START_REF] Dias | Viscous potentiel free-surface flows in a fluid layer of finite depth[END_REF]: the

symbols of B = νId, B = -ν ∂ 2 ∂x 2 , B = ν -∂ 2 ∂x 2 and B = ν(Id -∂ 2 ∂x 2
) -1 are very correctly approached by the method. BO equation is rewritten into its conservative form:

u t -∂ x H(u x ) + 1 p + 1 u p+1 = 0 x ∈ R, t > 0, u(x, 0) = u 0 (x) x ∈ R. (25)
When p = 1, BO possess an infinity of invariants, see [START_REF] Abdelouhab | Nonlocal models for nonlinear dispersive waves[END_REF], from which

I -1 (u) = udx, I 0 (u) = u 2 dx, I 1 (u) = 1 3 u 3 dx -uHu x dx, I 2 (u) = 1 4 u 4 dx - 3 2 u 2 Hu x , dx + 2 u 2 x dx.
When p ≥ 2, they are only 3 invariants:

I -1 (u) = udx, I 0 (u) = u 2 dx, I 1 (u) = 2 (p + 1)(p + 2)
u p+2 dx -uHu x dx, see [START_REF] Bona | Singularity formation in the generalized Benjamin-Ono equations[END_REF] .

We displayed a Sanz-Serna Scheme to simulate the equation and computed approximation to the symbol of the damping operator. We use both the spectral approach and the multi-grid one in finite differences. In Figure [START_REF] Dumont | Numerical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach[END_REF],we represent the estimate of the damping operator B = √ -∆ when spectral Fourier discretization and spectral SAM methods are applied. We observe that the symbol is correctly approached only for low modes: this is due to the spectral profile of the soliton (build only low frequencies). In Figure [START_REF] Dutykh | Viscous-potential free-surface flows and long wave modeling[END_REF], we have represented HSAM applied to finite differences discretization and hierarchical filtering. We observe as above, that the symbol is only well approached for very low frequencies. This is not surprising since the soliton is to poor in frequencies, the approximation of the symbol in the Fourier spectral case is also very limited to the low frequencies, as showed in Figure [START_REF] Dumont | Numerical analysis of a water wave model with a nonlocal viscous dispersive term using the diffusive approach[END_REF]. The numerical procedures proposed in the present work (SAM, MSAM and HSAM) allow to produce a simple, fast and satisfactory approximation of the symbol of dissipative and damping operators, using different technique of discretization (Spectral Fourier, Finite Differences and Finite Elements). Of course natural extensions to orthogonal polynomials of wavelet discretization are possible since they allow to produce a separation of set of frequencies for a given data.

It appears (Proposition 3.4 and illustrations) that the three important ingredients from which depends the efficiency of the approach are:

• The discretization scheme

• The choice of the filtering procedure (Multi-grid or so)

• The choice of the sample (Rademacher, experimental data) Working with samples data that cover an important band of frequencies seems essential to capture the part of the symbol associated to high frequency numbers. This was pointed out when considering KdV then BO equations.

As further developments, we can consider more models in hydrodynamics (such as from Boussinesq's family), and also the application of extrapolation formula to improve the SAM method, as pointed in Remark 3.1.

We can also consider the case in which the damping operator B is time dependent and derive corresponding SAM approximation:

du dt + Au + F (u) + B(t)u = 0, t ∈ [0, T ]. (26)
The object here to be approached is the "spectrogram" of B. We decompose the time interval [0, T ] as [0, T ] = ∪ q-1 p=0 [t p , t p+1 ] with t 0 < t 1 < • • • < t P = T, Let ∆t > 0 be given, we consider for simplicity the case t p+1 -t p = cst = r∆t for a given integer r. The sequence w (k) satisfies

B k W (k+1) = - W (k+1) -W (k) ∆t -AW (k+1) -F (W (k+1) ) = R (k+1) , k = 0, • • • , N = qr.
We propose to apply the identikit spectral method as follows: let p ∈ 1, • • • , q be fixed • Compute the sequence R (k+1) = -W k+1 -W (k) ∆t -AW (k+1) -F (W (k+1) ), k = (p -1)r + 1...pr

• We decompose W (k+1) ) and R (k+1) using frequency filters Π j , j = 1, • • • , m (spectral, multigrid or so)

• We compute S (k,p) j =< Π j R (k+1) , Π j W k > and T k j =< Π j W (k+1) , Π j W (k) >, k = (p -1)r + 1...pr

Figure 1 :

 1 Figure 1: Matrix estimation with Hadamard and Bernoulli samples n × s, random matrix n × n, n = 4096, s = 16

Figure 2 :

 2 Figure 2: Trace estimation with Hadamard and Bernoulli samples n × s, random matrix n × n, n = 4096, s = 16

Figure ( 3 Figure 3 :

 33 Figure 3: Symbol Identikit and its polynomial fitting : the 1D case n = 127, A is n × n

Figure 4 :Figure 5 :

 45 Figure 4: Spectral Identikit and its polynomial fitting : the 2D case n = 64, A is n 2 × n 2

Figure 6 :

 6 Figure 6: Approximation of the symbol of band-limited filter Nmax=290, N=30, m=15

Figure 7 :Figure 8 :

 78 Figure 7: Approximation of the Spectrum of the 1D Finite differences laplacian matrix, with 2d and 4th order discretization schemes. Nmax=256, N=20, m=6

Figure 9 :

 9 Figure 9: Stokes Problem with α = 0. Spectrum Identikit of the Schur Complement. Re = 1000, m = 12, dof=3025, N=50. The approach symbol is σ m (x, y) = 932.79382.1243 x + 2.1243 y -0.0278(x 2 + y 2 )

Figure 10 :

 10 Figure 10: Stokes Problem with α = 1000. Line 1: Spectrum Identikit of the Schur Complement in 2D (left) and 1D (right); Line 2 SAM estimate of the symbol of the preconditioned Uzawa Matrix in the complex plane. Re = 1000, m = 12 , n d.o.f.=3025, N=50. The approach symbol is σ m (x, y) = 8.1173 -0.2429 x -0.2429 y + 0.0111(x 2 + y 2 )

Figure 11 :

 11 Figure 11: Top : the successive meshes. Bottom: Homogeneous Dirichlet Problem on Ω =]0, 1[ 2 . Comparison of Identikit Spectrum of the stiffness matrix and of the exact symbol of the negative laplacian. P 1 elements, n.d.o.f = 3721

Figure 12 :

 12 Figure 12: Top : the successive meshes. Bottom: Homogeneous Dirichlet Problem on the unit Disk. Comparison of SAM of the stiffness matrix and of the exact symbol of the negative laplacian. P 1 elements, n.d.o.f = 348

Figure 13 :

 13 Figure 13: Two levels hierarchical decompositions on the interval with periodic B.C. (left) and on the unit square with Dirichlet BC (right)

Figure 14 :Figure 15 :

 1415 Figure 14: Multilevel Decomposition of the signal u(x) = sin(100x(1 -x)), N = 127. Second order interpolation

Figure 16 :

 16 Figure 16: Finite differences discretization. HSAM for the 1D Laplacian, n = 511, A is n × n-4th order discretization scheme

Figure 18 :

 18 Figure 18: Comparisons of the SMA method when using 2d, 4th and 6th order discretization schemes. N = 50, n = 127

Figure 19 :

 19 Figure 19: KdV equation-Spectral Fourier discretization in space. Energy spectrum of the initial soliton datum. ν = 1, ∆t = 0.001, L = 100, n = 128 ; Line 1 γId (left) and ∆ (right), Line 2 √ -∆ (left), (1 -∆) -1 (right)

Figure 20 :

 20 Figure 20: KdV equation. Approximation of the symbol of the damping operators by spectral local means. ν = 1, ∆t = 0.001, L = 100, N = 128 ; Line 1 γId (left) and ∆ (right), Line 2 √ -∆ (left), (1 -∆) -1 (right)

Figure 22 :

 22 Figure 22: Benjamin Ono equation. Approximation of the symbol of the damping operator √ -∆ by spectral local means (left), Spectrum energy of the initial datum (right) n = 64.

Figure 23 :

 23 Figure 23: Benjamin Ono equation. Approximation of the symbol of the damping operator √ -∆ by spectral local means. HSMA and 4th order finite differences discretization. n = 64

  p) j , j = 1, • • • , m.

We now consider the multi-grid case in 1D, base on a 4th order interpolation. The discretization is realized with 4th order compact finite differences schemes, we proceed as in see [START_REF] Calgaro | Séparation des échelles et schémas multiniveaux pour les équations d'ondes non-linéaires[END_REF].

, the symbol is only well approach for the lower and the higher frequencies, however it appears to be a low-pass filter.

Application to the Benjamin-Ono equation

The Benjamin-Ono equation (BO) describes one-dimensional internal waves in deep water, see [START_REF] Benjamin | Internal waves of permanent form in fluids of great depth[END_REF][START_REF] Ono | Algebraic solitary waves in stratified fluids[END_REF]. It reads as

with p ≥ 1; here H is the Hilbert transform:

The spectrum of B in the interval [t p , t p+1 ] will be approached by the piece-wise constant function based in the sequence λ (p) j , j = 1, • • • , m.


Michel Paul Newts

come