Dan Stefanoiu

Ioan Doré Landau
email: landau@gipsa-lab.grenoble-inp.fr

Mihai Cornoiu
email: mihai.cornoiu@gmail.com

Identification of an Active Suspension in the Presence of Narrow Band Disturbances

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

The active suspension system is already under study at many great automobile companies. This is the reason some research projects are developed in the academic filed as well. The GIPSA Laboratory from Grenoble (France) is hosting such a system, on both research and didactical purposes. Fig. 1 illustrates a picture of that active system. Fig. 1. Photo of the active suspension system hosted by GIPSA Laboratory from Grenoble (France).

Such a system is quite complex and mainly consists of two sub-systems: a primary one and a secondary one, as Fig. 2 displays. The primary sub-system is actually the classical passive suspension. Its input is the primary force (seen as disturbance), while the output is the first residual force, as returned by the elastomere cone. (Forces are measured by means of accelerometers.) The cone is only able to partially attenuate the primary force. The secondary sub-system is the active extension of passive suspension by means of a controllable inertial actuator. This time, the input is the electrical power applied through the coil winding, in order to provide an attenuation force (the output), which is opposite with respect to the first residual force. By composing both outputs (the first residual force and the attenuation force), the second (final) residual force is obtained. Fig. 2 shows this resulted force only.

This paper focuses on the identification of the secondary path from active suspension, while working under the influence of narrow band disturbances of significant level. Subsequently, optimal controllers have been designed for the actuator, based on several identification models, in order to minimize the intensity of the final residual force. Since the signals between sub-systems are vibrations, the overall system exhibits a characteristic frequency of about 52 Hz and a critical resonance frequency of about 95 Hz. Therefore, the study mainly focuses on the behaviour of active suspension excited by narrow band disturbances, which are centred at different frequencies, including the resonance one. (Further stochastic disturbances (seen as different coloured noises) are under consideration.) This topic has not been addressed in previous contributions (see for example [Landau et.al, 2005]), where identification has been performed in the absence of such disturbances. However, the narrow band disturbances may occur during the commissioning of the suspension system. The article is structured as follows. The next section reviews some identification methods that have been employed in this study. (The ideal model is also emphasized.) Section 3 shortly describes the corresponding controllers. Section 4 is devoted to simulation results. Some concluding remarks and a succinct list of references complete the article.

IDENTIFICATION OF SUSPENSION SYSTEM

Here and hereafter, one assumes that the reader is accustomed with the terminology of System Identification and Digital Control Systems fields [START_REF] Söderström | System Identification[END_REF]], [START_REF] Landau | Identification et commande des systèmes[END_REF], [START_REF] Landau | Digital Control Systems -Design, Identification and Implementations[END_REF], [Stefanoiu et al., 2005]. Before reviewing the appropriate identification methods, one has to mention that the noiseless model of the suspension system has previously been identified by means of long Gaussian pseudo-random inputs, sampled at 800 Hz.

The model is of ARX type, with 14 zeros and 12 poles. In Fig. 3, the model spectrum has been drawn. One can notice the two peaks pointing to the characteristic (primary resonance) frequency (51.66 Hz) and the (secondary) resonance frequency (~95 Hz), respectively. The performance of the noisy identification models will be assessed according to this ideal frequency characteristic. The ARX model actually belongs to ARMAX class, generally expressed by the linear regression equation below:

T [] [] [] yn n en =+ φθ , n * ∀∈ N , (1)
where [] n φ is the regressors vector, θ is the parameters vector (unknown) and e is a (Gaussian) white noise. The two vectors have several compounds that match ARMAX class:

TT T T ,,, yn a un b en c ⎡⎤ ≡-⎣⎦ φφ φ φ , T TTT abc ⎡⎤ ≡ ⎣⎦ θ θθθ , (2)
where

[ ] T , [] [ 1 ] [ ] xm nx n x n m =- - φ A , n *
∀∈ N , for some signal x and structural index m * ∈ N , while θ includes the coefficients of the 3 polynomials A , B and C . Obviously, φ consists of acquired and/or estimated data (since the white noise cannot be measured separately). The ideal ARX model only includes the first two compounds for each vector in (2).

Since the MA model of coloured noise could be non accurate, the ARMAX model has been extended to a more general Filtered Input Filtered Noise (FIFN) model. This time, the noise filter exhibits zeros as well as and poles. Thus the C polynomial of MA model is replaced by the couple of polynomials {C, D} . The FIFN regression equation is however similar to (1), but with a coloured noise v instead of the white noise e . The noise v is obtained by an AR model, with D polynomial playing the filtering role (all poles). Two types of identification methods are bond to ARMAX/FIFN models (1) (of suspension system): open loop and closed loop. When the identification ignores the automatic controller, one can rely on the following methods [START_REF] Dauphin-Tanguy | Modélisation, identification et commande des systèmes[END_REF]: (Recursive) Extended Least Squares (RELS), Open Loop Output Error with Extended Prediction Model (XOLOE), Generalized Least Squares (GLS). As of closed loop methods (when the controller is considered), the following were selected [START_REF] Landau | Recursive Algorithms for Identification in Closed Loop: A Unified Approach and Evaluation[END_REF]: Closed Loop Output Error (CLOE), CLOE with Extended Prediction Model (XCLOE), Generalized CLOE (GCLOE). A very brief description of the 6 above mentioned methods is presented next.

The RELS method actually employs a two stages LS strategy. The regressors vector from ( 2) is not entirely known, because φ , whereas φ becomes φ . Now, the LS procedure is applied once more:

1 T 11 11 ˆˆˆ[ ] [] [] [][] NN nn Nn n n y n NN - == ⎛⎞ ⎛ ⎞ = ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ ⎝⎠ ⎝ ⎠ ∑∑ θφ φ φ (4)
where N * ∈ N is the number of I/O acquired data. The two stages above car recursively be repeated for increasing N .

The XOLOE method has small differences comparing to RELS. Its goal is to increase the convergence speed, by overcoming the bias between the true and estimated parameters faster. This time, C is replaced by CA -. The corresponding coefficients sequence ca δ=is of length max{ , } nn a n c δ =

. In definitions (2), the following replacements are made: The difference ŷ y is actually the output error (OE), which has to be minimized by means of LS method. The OE formally replaces the prediction error (PE) computed like in definition (3) and creates the extended prediction model. The GLS method refers to the FIFN model (eventually with C1 ≡ , i.e. without zeros) and can be implemented by means of the general Minimum Prediction Error (MPE), which is intrinsic recursive. The main problem is to estimate the coloured noise. This can be done by means of an approximate ARX model with long memory, like in (3), unless, now, not the white noise e is estimated, but the coloured noise v . It follows that the noise filter is of ARMA type and can be identified through RELS method, from estimated output data.

The closed loop group of identification methods require the controller configuration. The configuration depicted in Fig. 4 is rather classical, of R-S-T type. Here, ε plays the role of OE between the acquired output y and the simulated output ŷ . The plant and the model include both filters (for input and for noise). However, the model can only work with estimated values of the white noise. The R-S-T controller depends on the identification model. It is designed such that noises are attenuated (or rejected) and the output verifies some performance criteria, according to reference r . Consequently, the , un b φ vector in (2) has to be computed/estimated by means of the controller, starting from the reference signal. Clearly:

( ) ( ) ( ) ( ) 11 11 Tq Rq [] [] [] Sq Sq un rn yn -- -- =- , n * ∀∈ N . ( 5
)
The PAA block stands for Parameter Adaptation Algorithm and performs model updating. (Note that the noise filter can be missing from the scheme, which decreases the controller performance.)

The CLOE method operates with the ARX model ( CD1 ≡≡ ) and is based on minimisation of two types of OE: a priori ( 0 ε ) and a posteriori ( 1 ε ). They are defined by: T

0 [] [] [][] ny n n n ε= -φθ & T 0 [] [] [ 1 ][] ny n n n ε= - - φθ , (6)
for any n * ∈ N . The rightmost terms of definitions (6) are referred to as the a priori output 0 ŷ and the a posteriori output 1 ŷ , respectively. Ignoring the noises leads however to less accurate models, although the convergence is quite fast.

The XCLOE method approaches the ARMAX model ( D1 ≡ ) and operates with the closed loop predictor below:

( ) ( ) ( ) ( ) 1 11 1 Qq ˆ[] 1 Aq [] Bq [] [] Sq yn yn un n - -- - ⎡⎤ = -+ + ε ⎣⎦ , (7)
for any n * ∈ N . In definition (7), the Q polynomial is:

Q1( CA ) SB R ≡ +-⋅ -⋅. ( 8
)
Minimization of extended OE ŷ y -(with ŷ defined in ( 7)) is based on some RELS-like strategy. The GCLOE method deals with the FIFN model, where both polynomials C and D can be non unit. The identification procedure is quite complex and, usually, makes use of MPE algorithm. This is the only method that fully employs the structure of Fig. 4. Several other methods have been tested as well. But their performance was quite low, so they were removed.

ON CONTROLLER DESIGN

Two types of controllers have been adopted. The first one is already described in Fig. 4 and belongs to R-S-T class. This controller covers quite a large panoply of linear systems and can be designed starting from either open loop or closed loop identification models. The design procedure is largely described in various publications such as [START_REF] Dauphin-Tanguy | Modélisation, identification et commande des systèmes[END_REF] or [START_REF] Landau | Digital Control Systems -Design, Identification and Implementations[END_REF]. Before proceeding to the controller design, the user has to set some performance specifications. Usually, such specifications consists of the number and position of poles (real valued and/or complex conjugate pairs), together with their damping factors. For the active suspension system, poles are grouped in 3 categories: (a) two complex conjugate dominant poles that focus on the characteristic frequency (51.66 Hz), with a damping factor of 0.8; (b) four to six complex conjugate pairs of auxiliary poles that cover the frequency band [90,400] Hz with different damping factors (usually smaller than 0.8); (c) eleven to fifteen real valued auxiliary identical poles located on the interval [0.22, 0.4] . The number and location of auxiliary poles depends on the identification model. In case of closed loop approach, updating of controller, on one side, and model estimation/adaptation, on the other side, are alternatively performed, through iterative procedures. The second type of controller belongs to Q-parameterized class, also known as Youla-Kucera class [START_REF] Blondel | Simultaneous Stabilisation of Linear Systems[END_REF]. The corresponding identification and control scheme is drawn in Fig. 5. The controller aims to integrate a model of disturbances by means of Q polynomial. After few manipulations, from the scheme of Fig. 5, one obtains the noise-to-output equation:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 111 11 0 11 1 1 1 00 Aq S q Bq Qq Cq Aq S q Bq R q Dq ye --- -- -- - - - ⎡⎤ - ⎣⎦ ≡ + . (9)
Eq. ( 9) implies that rejection of disturbances poles can be achieved with the following factorization (a Bézout-like equation, in fact):

( ) ( ) ( ) ( ) ( ) The unknown Q and 0 Q polynomials are solutions equation ( 10). Similarly, the disturbances zeros can be removed by selecting 0 S and 0 R polynomials such that ( ) ( ) ( ) ( ) ( ) ( )

11 1 1 1 00 Sq B q Q q Qq D q -- - - - -= . ( 10
11 1 1 1 1 00 0 Aq S q Bq R q P q Cq -- - - - - +=, ( 11
)
where the 0 P polynomial has to be preset (by the controller design specifications). Notice that eq. ( 10) cannot be solved before eq. ( 11), since 0 S needs to be known. Of course, both equations rely on the identification model. So, the controller has to be updated accordingly.

SIMULATION RESULTS AND DISCUSSION

The plant of Fig. 1 has been stimulated with 3 types of narrow band disturbances, at: 75 Hz, 90 Hz and 95 Hz, respectively. The last one falls on the secondary resonance frequency range. The sampling frequency is 800 s F = Hz and the I/O data block counts 10 000 samples. The 6 identification methods provided valid identification models, as described in section 2. The structural indexes na , nb and nc vary in range 10:16, whereas nd is maximum 2. The controllers were designed according to the principles of section 3. From the start, one has to say that both types of controllers led to very similar performance. However, the Q-parameterized controller is easier to implement. To compare the models and especially to emphasize the open loop versus closed loop approaches, frequency characteristics of the identified main filter have been drawn, together with the ideal characteristic of Fig. 3. Also, the distance between each estimated frequency characteristic and the ideal one has been computed by means of error standard deviation (std), denoted by λ . More specifically, assume that H is the ideal frequency characteristic (actually, the spectrum only), while Ĥ is one of its estimates, both expressed in dB. Normally, they are continuous maps. In order to draw their variations, let them uniformly be sampled. The number of samples, K , should be sufficiently large to have a good visual resolution (say at least 500). Then the std λ (as distance between H and Ĥ ) is computed as follows:

( )

2 1 1 ˆ[] [] K k Hk Hk H H K = λ= - - - ∑ , (12)
where Ĥ H is the error between the two characteristics, whereas Ĥ H is the average (the arithmetic mean) of the error. The smaller λ , the closer H and Ĥ each other. samples. The red solid curve is the estimated characteristic ( Ĥ ), while the blue doted curve shows the ideal characteristic ( H ). For each couple of variations, the std λ of error between the two characteristics was computed (like in eq. ( 12)). One can see that the LS methods perform differently than XOLOE, in terms of std. Whilst the LS methods perform the best at 75 Hz disturbance (i.e with minimum std), XOLOE improves as the disturbance approaches the secondary resonance frequency of 95 Hz (i.e. the std decreases). It seems that GLS is the very best open loop method (with the smallest std in all 3 cases). The best structural indexes for open loop models are: 10 na = , 12 nb = and 12 nc = , for ARMAX models identified with RELS or XOLOE methods; 10 na = , 12 nb = , 1 nc = and 12 nd = , for FIFN models identified with GLS method. The corresponded controller was of R-S-T type.

Figures 9, 10 and 11 focus on the results obtained with closed loop models and methods. The same graphical conventions and resolution like in case of open the loop approach are adopted, excepting for the solid curve, which turned in green. Apparently, the XCLOE method leads to the best models in terms of std λ . However, one can see that those models are not as accurate as the FIFN models estimated with GCLOE in 0-200 Hz subband (where disturbances lie). Therefore, one can consider that GCLOE is seemingly the best. This time, the optimal structural indices are: 14 na = , 14 nb = , 2 nc = and 2 nd = . The FIFN identification model also led to an efficient Youla-Kucera controller.

The performance of disturbances rejection is demonstrated in Figures 12 and13. Both figures display the variations ranged on 2 columns. On the left side, the residual forces are drawn. On the right side, power spectral densities (psd) are depicted, with a graphical resolution of 512 samples (like within the previous figures). The dotted (blue) curve corresponds to psd of residual force when no controller is integrated in the active suspension system. The solid variations are actually the psd of residual forces from the left (red for open loop, green for closed loop models). Only the best models were considered, i.e. of FIFN type. The std of error between solid and dotted curves (denoted by σ , this time, but computed just like in eq. ( 12)) was also evaluated. In this context, as overall performance, the higher σ , the better disturbances rejection. However, the quality of disturbances rejection is not entirely revealed by σ . The narrow band of about 50 Hz, centred on each disturbance peak (the dotted blue one in the figures), should be focused as well. Usually, when the disturbance attenuation is quite strong on the narrow band, some other disturbances outside the band could be amplified. This is the reason some identification models are only apparently better than others, when solely considering the overall std, σ . The rejection performance on the narrow band is assessed by nb σ , which is the std computed for the psd error restricted to the narrow band only. In general, nb σ> σ , since the number of psd samples is significantly smaller over the narrow band than over the whole band (the denominator in eq. ( 12) decreases).

Fig. 12 shows the performance of R-S-T controller when using FIFN-GLS (open loop) models. According to std σ , disturbances at 90 Hz are the least rejected, as Fig. 7 already predicted. However, the std nb σ points to the disturbances at 95 Hz as the least rejected, which is a more accurate result.

(Disturbances near the resonance frequency are harder to reject.) Fig. 13 illustrates the superior performance of Q-parameterized controller when using FIFN-GCLOE (closed loop) models, especially for 75 and 90 Hz disturbances. Again, the disturbances at 95 Hz are the least rejected. This time, both σ and nb σ lead to this conclusion.

CONCLUDING REMARKS

The case study developed in this paper concerned the identification and control of an active suspension system in the presence of narrow band disturbances at significant level.

After testing several identification models and two types of controllers, one can say that the FIFN model identified by means of GLS or GCLOE methods allowed us to design the best controller in terms narrow band disturbances rejection. The Q-parameterized controller structure could be preferred to the R-S-T structure only for implementation purposes, because both types of controllers perform about the same. As future research, non linear identification models are under consideration, starting from the analytical system of equations that describe the active suspension dynamics.

Fig. 2 .

 2 Fig. 2. Functional scheme of the active suspension system.

Fig. 3 .

 3 Fig. 3. Ideal spectrum of the active suspension system.

  . Therefore, the first stage is aiming to identify an approximate ARX model, in order to estimate the noise with. Such a model has already been identified: the ideal one, with 12 na = and 14 nb = . Another model (with larger structural indexes) can also be identified from I/O data. Then the white noise can be estimated as the prediction error:

  = φθ , n * ∀∈ N is the noise free simulated output, with φ estimated exactly like in RELS.

Fig. 4 .

 4 Fig. 4. Closed loop configuration with R-S-T controller.

Fig. 5 .

 5 Fig. 5. Closed loop configuration with Youla-Kucera parameterized controller.

Figures 6 ,

 6 Figures 6, 7 and 8 illustrate the frequency characteristics of open loop models, for the three types of disturbances, with a resolution of 512 K =samples. The red solid curve is the estimated characteristic ( Ĥ ), while the blue doted curve shows the ideal characteristic ( H ). For each couple of variations, the std λ of error between the two characteristics was computed (like in eq. (12)). One can see that the LS methods perform differently than XOLOE, in terms of std. Whilst the LS methods perform the best at 75 Hz disturbance (i.e with minimum std), XOLOE improves as the disturbance approaches the secondary resonance frequency of 95 Hz (i.e. the std decreases). It seems that GLS is the very best open loop method (with the smallest std in all 3 cases). The best structural indexes for open loop models are: 10 na = , 12 nb = and 12 nc = , for ARMAX models identified with RELS or XOLOE methods; 10 na = , 12 nb = , 1 nc = and 12 nd = , for FIFN models identified with GLS method. The corresponded controller was of R-S-T type.

Fig. 6 .

 6 Fig. 6. Open loop models for 75 Hz disturbances.

Fig. 7 .

 7 Fig. 7. Open loop models for 90 Hz disturbances.

  Fig. 8. Open loop models for 95 Hz disturbances.

Fig. 9 .

 9 Fig. 9. Closed loop models for 75 Hz disturbances.

Fig. 10 .

 10 Fig. 10. Closed loop models for 90 Hz disturbances.

Fig. 11 .

 11 Fig. 11. Closed loop models for 95 Hz disturbances.

Fig. 12 .

 12 Fig. 12. Residual forces and spectra for FIFN-GLS models.

Fig. 13 .

 13 Fig. 13. Residual forces and spectra for FIFN-GCLOE models.


Vincent Renner
email: vincent.renner@univ-lyon2.fr

Jesús Fernández-Domínguez

Coordinate compounding in English and Spanish 1 Coordinate compounding in English and Spanish

Keywords: word-formation, compounding, coordination, English, Spanish

Coordinate compounding is a process which has been only sporadically considered in wordformation studies on Germanic and Romance languages. In this paper, we compare the situation in English with that in Spanish as far as formal structure and semantics are concerned. To this end, an operational definition of coordinate compounding is first provided, after which a semantic classification of coordinate compounds is developed. For each type, variants are discussed, representative examples are provided, and similarities and dissimilarities between the two languages are highlighted. We finish by raising the issues of recursiveness and productivity in coordinate compounding and we examine the possibilities of expansion and profitability of the different patterns. Our analysis reveals marked similarities between English and Spanish coordinate compounding at various levels, hinting at the possible existence of a cross-linguistic set of common core features in the class of coordinate compounds in both Germanic and Romance.

Introduction

Romance languages are said to resort to the word-formation process of compounding to a much lesser extent than Germanic languages. This is commonly illustrated by the fact that many Germanic compounds correspond to simplex or derived lexemes in Romance (see e. Coordinate compounds represent only a small fraction of the whole class of compounds in both the Romance and Germanic families. [START_REF] Guevara | Searching for universals in compounding[END_REF]'s approximations are respectively 16 and 12 percent, and more conservative estimates, which rely on random samples of dictionary-listed compounds, put the percentage as low as about 2 percent for English (Arnaud 2002: 4;Berg 2009: 134). The subclass of Spanish coordinate compounds being relatively small, it has unsurprisingly received somewhat scant attention from word-formationists, and in this article we endeavour to repair this neglect, to delve into the formal and semantic details of these constructions and to compare the situation in Spanish with that in English, a language for which fine-grained analyses of coordinate compounding are already available (e.g. [START_REF] Renner | On the semantics of English coordinate compounds[END_REF][START_REF] Bauer | Co-compounds in Germanic[END_REF].

The general structure of the article is as follows. We first present a definition of coordinate compounding and associate it with a unitary metalinguistic test sentence (Section 2), and we then attempt a semantic characterization of coordinate compounds based on their lexical input structure, i.e. noun-noun, adjective-adjective, verb-verb and adverb-adverb units (Section 3). For each subtype, several illustrative examples are presented, the English and Spanish variants are contrasted, and the most prominent features are pointed out. In Section 4, we focus on the issues of productivity and recursiveness to evaluate the creative potential of this word-formation device.

Definition of coordinate compounding

There has been a recent surge of cross-linguistic interest in coordinate compounding, which has led to a multiplication of classificatory and terminological discordances in the literature (for a discussion, see Wälchli 2005;[START_REF] Bauer | Dvandva[END_REF][START_REF] Renner | On the semantics of English coordinate compounds[END_REF][START_REF] Arcodia | Coordinating compounds[END_REF], and one should be aware that various authors may use the same terms, e.g. dvandva or co-compound, to refer to different types of compounds. In this study, the stance adopted is to use terms which are as transparent as possible and which closely reflect the semantics of each type and subtype.

The operational definition of coordinate compounds proffered in this research is the following:

those units composed of co-hyponymous lexemes for which the linear order of elements is deemed to be grammatically irrelevant (the factors coming into play are phonological and pragmatic). In some cases, the two variants X+Y and Y+X have been institutionalized (e.g. Eng. lease-lend / lendlease; Sp. oxidación-reducción / reducción-oxidación),1 which leads us to claim that membership in the category of coordinate compounds can be materialized through a unitary metalinguistic test sentence, "X+Y could be named Y+X" for English and "X+Y podría decirse Y+X" for Spanish, as in the following examples:

i. nouns -Fighter-bomber could be named bomber-fighter.

-*Cherry tomato could be named tomato cherry.

-Caza-bombardero2 podría decirse bombardero-caza.

-*Mono araña3 podría decirse araña mono.

ii. adjectives -Blue-green could be named green-blue.

-*Icy-cold could be named cold-icy.

-Obsesivo-compulsivo 4 podría decirse compulsivo-obsesivo.

-*Jurídico-laboral 5 podría decirse laboral-jurídico.

iii. verbs

-Work-study could be named study-work.

-*Blow-dry could be named dry-blow.

-Duermevela6 podría decirse veladuerme.7

Some compounds fail to pass the above test because of their directional semantics. This is the case of asynchronous compounds such as homicide-suicide, drag-drop, servicio-volea ('serve + volley')

and copia-pega ('copy + paste'), which refer to events or actions which are consecutive in time, and also that of units such as English-Russian or español-francés ('Spanish + French') -dubbed "translative" by Bauer (2001: 700) -when modifying nouns such as dictionary, glossary, manual / diccionario, glosario, manual. The two elements being on an equal footing semantically (neither can be said to be the sole head of the compound), the non-subordinate nature of such compounds is, however, unambiguous.

4 'obsessive + compulsive'.

5 'legal + labor (adj.) = pertaining to labor law'.

Semantic typology

In this section, we propose a semantic typology of coordinate compounds making use of test sentences. Most -but not all -types are attested in the two languages; Spanish most notably contrasts with English by exhibiting a somewhat richer semantic typology of verb-verb compounding.

noun-noun compounds

Three semantic classes of noun-noun compounds are identified:

i. multifunctional compounds such as owner-operator, speaker-hearer, gymnasiumauditorium / compositor-arreglista ('composer + arranger'), lavadora-secadora ('washer + dryer'), sofá-cama ('sofa + bed'), which can be felicitously inserted in the test sentence "An X+Y is an X who/which is also a Y" / "Un(a) X+Y es un(a) X que también es un(a) Y";

ii. hybrid compounds such as troutperch, mashie-niblick, nickel-iron, folk punk / aguanieve ('rain + snow = sleet'), gomorresina ('gum + resin'), cerapez ('wax + pitch'), pasitrote 8 ('walk + trot = ambling movement'), which can be felicitously inserted in the test sentence "X+Y is a blend of X and Y" / "X+Y es una mezcla de X y de Y";

iii. additional compounds such as robbery-homicide, space-time, tractor-trailer, Luke-Acts9 / compraventa ('buying + selling'), aguaviento ('rain + wind'), casatienda ('house + store'), 8 The vowels i and o are occasionally used compound-internally in Spanish (e.g. capisayo < capa + sayo tarde-noche ('evening + night'), which can be felicitously inserted in the test sentence "An X+Y is an X plus a Y"/ "Un(a) X+Y es la suma de un(a) X y de un(a) Y".

Multifunctionals and hybrids have a strong affinity -they share the property of being constituted of elements whose denotata are perceptually merged and are, in this respect, unlike additionals, whose denotata are simply juxtaposed -, but we have chosen to distinguish the two classes due to the clear difference in acceptation between the test sentences. In Spanish, there is a formal dissonance between multifunctionals, which never seem to admit a medial compound marker, and hybrids, which include constructions with a marker (e.g. pasitrote, ceromiel, gomorresina), which also goes against the idea of merging the two categories.

In the two languages, exocentric compounds are attested only in the class of additionals: in English, a penny-farthing, also known as an ordinary, is a bicycle with a very large ("penny-like") front wheel and a small ("farthing-like") rear wheel; in Spanish, capigorra ('cloak + cap') denotes an idle vagrant through two conceptually salient sartorial details and ajoqueso ('garlic + cheese') refers to a dish prepared with garlic and cheese, but also onion and pepper.

adjective-adjective compounds

Two semantic classes of adjective-adjective compounds are identified:

i. additional compounds such as deaf-blind, manic-depressive, pale-dry, solar-lunar, sweetsour / agripicante ('sour + spicy'), anchicorto ('wide + short'), fisicoquímico ('physical + chemical'), tontiloco ('dumb + crazy'), which can be felicitously inserted in the test sentence "To be X+Y is to be X and to be Y" / "Ser X+Y es ser X y ser Y";

ii. hybrid compounds such as blue-red (= purplish), medium-rare, dressy-casual / verdiazul ('green + blue'), verdiseco ('green + dry = half-dry'), which can be felicitously inserted in the test sentence "To be X+Y is to be about midway between X and Y" / "Ser X+Y es estar a medio camino entre X e Y".

Exocentric adjective-adjective compounding is attested only in the case of Spanish, in claroscuro ('light + dark = chiaroscuro').

verb-verb compounds

Five The semantics of some verb-verb compounds are not beyond debate among morphologists and lexicographers. This is for example the case of arrancasiega in Spanish, Bustos Gisbert (1986: 318) and Alvar Ezquerra (1993: 23) considering that it has a hybrid reading, and the editors of the Diccionario de la lengua española a disjunctive one, or that of freeze-dry in English, Renner (2006: 59-62) and the editors of the American Heritage dictionary of the English considering that the compound is a troponym of preserve, and therefore a coordinate compound, and Bauer (2010: 204) and the editors of the Oxford English Dictionary (= OED) considering it is a troponym of dry, and hence a subordinate compound. 11 Verb-verb compounds also markedly differ in the two languages in so far as Spanish accepts only exocentric nominal outputs whereas English accepts verbal, adjectival and nominal ones (respective examples are scrunch-dry, cook-chill and lendlease).

adverb-adverb compounds

Adverb-adverb compounding is highly marginal in English, and it is unattested in Spanish. The 11 Troponymy is an asymmetrical relation of sense between two verbs. This relation is explained in terms of manner elaboration, as in the formula "To run is to move in some manner", where run is said to be a troponym of move (for a discussion, see [START_REF] Fellbaum | On the semantics of troponymy[END_REF].

Coordinate compounding in English and Spanish 10 only unequivocal English example is softly-softly, an exocentric reduplicative compound with an adjectival output. Willy-nilly is an adverb and it could be argued that there is a coordinate relation between the two elements -the OED gives nilly-willy as an established variant -, but the construction is not a compound in the narrow sense of the term, as nilly is not an English lexeme. It is rather a pseudo-compound, which has arisen through univerbation (< "will I, nill I"). Similarly, following Cook (1969: 135), Jackson and Zé Amvela (2000: 83) list in-to and through-out as instances of coordinate adverbial compounding, but this classification is, in our view, strongly objectionable as the first example is not known to be an adverb, and the second one does not fit the definition given in Section 2 -the OED gives out-through as an obsolete, chiefly Scottish variant of throughout, but linear reversibility is synchronically impossible.

Productivity and recursiveness

English and Spanish compounds are said to differ widely in their general proclivity to both productivity, which can be broadly defined as the property whereby a word-formation process has the potential for the rule-governed creation of an indefinite number of words, and recursiveness, which refers to the possibility of lengthening a compound by using it as the base in a new compounding operation. This is true in the general case, i.e. that of subordinate compounds, but the situation is more complex and strikingly similar in the two languages as far as coordinate compounding is concerned. The relative productivity of the types listed in the previous section is as follows:

a high degree of productivity for multifunctional nouns, and also for some classes of additional adjectives -mainly neoclassical compounds, e.g. oronasal, Franco-British, hispano-inglés ('Hispano-English'), but also colour adjectives such as blanquinegro ('white + black'), azuliverde ('blue + green') and grisirojo ('grey + red') in Spanish; 12 -a low degree of productivity for hybrid and additional nouns;

a near-zero degree of productivity for verb-verb compounds and hybrid adjectives.

We would like to posit as a general principle that the degree of recursiveness of a given coordinate compounding pattern is correlated to the degree of productivity of that same pattern, which would explain why the distribution of recursive types runs parallel to that of productive types:

relatively frequent recursiveness for multifunctional nouns -e. absence of recursiveness for verb-verb compounds and hybrid adjectives.

12 The attestation of the last two adjectives goes against Val Álvaro (1999: 4813)'s contention that compounds with a medial marker cannot contain a final-stressed and consonant-final left element such as azul-('blue') or gris-('grey').

13 Several triple units with a double -i-compound marker (e.g. rojinegriblanco ['red + black + white'],

blanquirojiverde ['white + red + green']) are attested in Spanish, which disproves [START_REF] Rainer | Compounding in Spanish[END_REF] 139-140)'s claim that such formations are absent in the language.

14 North-northwest and sur-sureste are exceptional examples of hierarchical (right-branching) coordinate compounding.

In the case of multifunctionals and additional adjectives, longer units, like the following quadruplets, are also attested: singer-songwriter-producer-arranger, mother/father/nurturer /provider, Anglo-Franco-Turco-Russian / actor-director-escritor-estrella ('actor + director + writer + star'), marxista-leninista-estalinista-maoísta ('Marxist + Leninist + Stalinist + Maoist'), socioeconómico-político-cultural ('socio-economico-politico-cultural').

One substantial difference between subordinate and coordinate compounds lies in the relative ease of their semantic interpretation. Recursiveness in subordinate compounds gives rise to long units which may prove difficult to interpret, e.g. Berg (2006: 198)'s six-noun nonce compound child language acquisition research group member. Our claim is that this is not the case to the same extent with coordinate compounds because their internal structure and semantics remain relatively transparent as long as all the members of the compound are obvious co-hyponyms, as in the case of, for example, singer/composer/pianist/organist or actor-bailarín-mimo.

Conclusion

This study has brought to the fore the fact that English and Spanish coordinate compounds display striking similarities: a vast majority of the different lexical and semantic types are attested in both English and Spanish, and the relative proclivity to productivity and recursiveness of these types is remarkably analogous in the two languages. Noteworthy dissimilarities centre around verb-verb compounding: in Spanish, a somewhat richer semantic typology can be established and the existence of non-nominal outputs seems to be categorically excluded.

The remarkable number of shared traits in the two languages hints at the possible existence of a cross-linguistic set of common core features in the class of coordinate compounds in both Germanic and Romance. Investigating this thoroughly seems to us a worthy avenue of inquiry, but it is left for future research.

  g. Paillard 2000: 55-56 for French vs. English); examples of this contrast in Spanish and English are gallina ~ gallinero / hen ~ henhouse and plátano ~ platanal / banana ~ banana plantation.

[

  'cloak + smock = bishop's mantelletta']; ceromiel < cera + miel ['wax + honey = mixture of wax and honey']). We indicate the presence of such compound markers (for a broad discussion of the concept, see Ralli 2008) through underlining.

  semantic classes of verb-verb compounds are identified, the last two being restricted to

	boil = busybody') and picapica ('itch + itch = itching powder');
	v. hybrid compounds, which fit the test sentence "X+Y refers to a state/process which is about
	midway between Xing and Ying" / "X+Y se refiere a un estado/proceso que está a medio
	camino entre X e Y"; to the best of our knowledge, English has no such thing as hybrid verb-
	verb compounding; Spanish has duermevela ('sleep + be awake = state intermediate between
	sleeping and being awake').
	Spanish:
	iv. reduplicative-intensive compounds, in which X is formally equivalent to Y, and whose
	semantics are "To X a lot" / "X mucho"; to the best of our knowledge, English has no such
	thing as reduplicative-intensive verb-verb compounding; 10 Spanish has bullebulle ('boil +

i. synchronous compounds, which fit the test sentence "X+Y refers to the concept of 'Xing and Ying at the same time'" / "X+Y refiere al concepto de 'X e Y al mismo tiempo'", as in work-study / ganapierde ('win + lose = game of loser-takes-all'); ii. asynchronous compounds, which fit the test sentence "X+Y refers to the concept of 'Xing and then Ying'" / "X+Y se refiere al concepto de 'X y luego Y'", as in copy-paste, stop-start, fly-drive / corta-pega ('cut + paste = cut-paste operation'), quitapón ('remove + put = headstall ornament for mules');

iii. disjunctive compounds, which fit the test sentence "X+Y refers to the concept of 'Xing or Ying'" / "X+Y se refiere al concepto de 'X o Y'", as in pass-fail / arrancasiega ('uproot + cut = process of uprooting or cutting');

'oxidation + reduction'.

'fighter + bomber'.

'monkey + spider = spider monkey'.

'sleep + be awake = state intermediate between sleeping and being awake'.

To the best of our knowledge, there is no such thing as subordinate verb-verb compounding in Spanish.

In Biblical studies, Luke-Acts refers to the hypothetical composite work of the Gospel of Luke and the Acts of the Apostles in the New Testament.

The adjectival compound win-win arises from the juxtaposition of two identical verb forms, but its semantics are not reduplicative-intensive; it belongs among synchronous compounds.

Wälchli, B. 2005. Co-compounds and natural coordination. Oxford: Oxford University Press.


Claire David

Pierre Sagaut

Spurious caustics of Dispersion Relation Preserving schemes

Keywords: Spurious caustics of Dispersion Relation Preserving schemes Dispersion, numerical schemes, spurious caustics, numerical group velocity, spurious error pile-up Dispersion, numerical schemes, spurious caustics Dispersion, numerical schemes, spurious caustics

published or not. The documents may come

Introduction

The analysis and the control of numerical error in discretized propagation-type equations is of major importance for both theoretical analysis and practical applications. A huge amount of works has been devoted to the analysis of the numerical errors, its dynamics and its influence on the computed solution (the reader is referred to classical books, among which [6,14,9,10]). The emergence of Dispersion-Relation-Preserving (DRP) schemes [4]), which have the same dispersion relation as the original partial difference equations, enables one to have very accurate high order finite difference schemes. The two sources of numerical error are the dispersive and dissipative properties of the numerical scheme, which are very often investigated in unbounded or periodic domains thanks to a spectral analysis. It appears that existing works are mostly devoted to linear, one-dimensional numerical models, such as the linear advection equation

∂u ∂t + c ∂u ∂x = 0 ( 1
)
where c is a constant uniform advection velocity. The two sources of numerical error are the dispersive and dissipative properties of the numerical scheme, which are very often investigated in unbounded or periodic domains thanks to a spectral analysis. Following this approach, a monochromatic wave is used to measure the accuracy of the scheme. Such a tool is very efficient and provides the user with a deep insight into the discretization errors. But some results coming from practical numerical experiments still remain unexplained, despite the linear character of the discrete numerical model.

As an example, let us note the sudden growth of the numerical error for long range propagation reported by Zingg [16] for a large set of numerical schemes, including optimized numerical schemes.

The usual modal analysis is almost always applied to monochromatic reference solutions, with the purpose of analyzing the error committed on both their amplitude and their phase, leading to classical plots of the relative error as the function of the Courant number and/or the number of grid points per wavelength. Therefore, dispersive phenomena associated to polychromatic solutions are usually not taken into account.

The present paper deals with the analysis of linear dispersive mechanism which results in local error focusing, i.e. to a sudden local error burst in the L ∞ norm for polychromatic solutions. This phenomena is reminiscent of the physical one referred to as the caustic phenomenon in linear dispersive physical models [15], and will be referred to as the spurious caustic phenomenon hereafter. It extends our previous work [2] to DRP schemes. The present analysis is restricted to interior stencil, and the influence of boundary conditions will not be considered. The paper is organized as follows. Main elements of caustic theory of interest for the present analysis are recalled in section 2. DRP schemes are presented in section 3. Their caustical analysis is exposed in section 4. A numerical example is presented in section 5.

Caustics

The solution of Eq. ( 1) is taken under the form:

u(x, t, k) = e j (k x-ω t) (2)
where j denotes the complex square root of -1, ω = ξ ω + i η ω is the complex phase, and k the real wave number: where λ 0 denotes the wavelength, i.e. For dispersive waves, it is recalled that the group velocity V g (k) is defined as

k = 2 π λ 0 (3)
V g (k) ≡ ∂ξ ω ∂k ( 4
)
A caustic is defined as a focusing of different rays in a single location. The equivalent condition is that the group velocity exhibits an extremum, i.e. there exists at least one wave number k c such that

∂V g ∂k (k c ) = 0 (5)
The corresponding physical interpretation is that wave packets with characteristic wave numbers close to k c will pile-up after a finite time and will remain superimposed for a long time, resulting in the existence a region of high energy followed by a region with very low fluctuation level. The linear continuous model Eq. ( 1) is not dispersive if the convection velocity c is uniform, and therefore the exact solution does not exhibits caustics since the group velocity does not depend on k. The discrete solution associated with a given numerical scheme will admit spurious caustics, and therefore spurious local energy pile-up and local sudden growth of the error, if the discrete dispersion relation is such that the condition ( 5) is satisfied. For a uniform scale-dependent convection velocity, such spurious caustics can exist in polychromatic solutions only, since they are associated to the superposition of wave packets with different characteristic wave numbers. Set:

k = ϕ σ c τ ( 6
)
where σ is the cfl number, defined as σ = cτ /h, h and τ denoting respectively the mesh size and time step.

The non-dimensional wave number is defined as ϕ = kh, where k is the wave number of the signal under consideration. The general dispersion relation associated with the discrete scheme enables us to obtain the corresponding group velocity, given by:

V g = h ∂ξ ω ∂ϕ (7)
The numerical solution will therefore admit spurious caustics if

∂V g ∂k = ∂V g ∂ϕ ∂ϕ ∂k = 0 ⇐⇒ ∂V g ∂ϕ = 0 (8)
Since we deal with a constant speed transport equation:

∂V g ∂k = ∂V g ∂ϕ ∂ϕ ∂k = 0 (9)
∂V g ∂ϕ = 0 (10)
The numerical solution will therefore admit spurious caustics if

∂V g ∂ϕ = 0 (11)
The corresponding values of ϕ and k will be respectively denoted ϕ c and k c . Spurious caustics are associated with characteristic lines given by

x t = U c ( 12
)
where

U c = V g (ϕ c ) (13)

DRP schemes

The Burgers equation:

u t + c u u x -µ u xx = 0, (14)
c, µ being real constants, plays a crucial role in the history of wave equations.

It was named after its use by Burgers [1] for studying turbulence in 1939.

For µ = 0, (14) reduces to the linear advection equation

u t + c u x = 0, (15)
which will be the object of our study.

For given natural integers i ∈ [0, n x ], n ∈ [0, n t ], a linear finite difference scheme for this equation can be written under the form:

F u m n li = 0 (16)
where:

u mn l i = u (l i h, m n τ ) (17)
and where F is a linear function of the u mn li ,

l i ∈ [1, n x ], m n ∈ [0, n t ]. Common values for l i are i -1, i, i + 1, common values for m n are n -1, n, n + 1.
F depends on the mesh size h, and the time step τ .

A numerical scheme is specified by selecting an appropriate expression of F. Then, depending on them, one can obtain optimum schemes, for which the error will be minimal. m being a strictly positive integer, the first derivative ∂u ∂x is approximated at the l th node of the spatial mesh by:

( ∂u ∂x ) l m k=-m γ k u n i+k ( 18
)
Following the method exposed by C. Tam and J. Webb in [4], the coefficients γ k are determined requiring the Fourier Transform of the finite difference scheme (18) to be a close approximation of the partial derivative ( ∂u ∂x ) l . ( 18) is a special case of:

( ∂u ∂x ) l m k=-m γ k u(x + k h) ( 19
)
where x is a continuous variable, and can be recovered setting x = l h. Denote by ω the phase. Applying the Fourier transform, referred to by , to both sides of (19), yields:

j ω u m k=-m γ k e j k ω h u (20)
Comparing the two sides of (20) enables us to identify the wavenumber λ of the finite difference scheme (18), i.e. the number of wavelengths per unit distance, and the quantity

1 j m k=-m
γ k e j k ω h , one obtains:

λ = -j m k=-m γ k e j k ω h (21)
To ensure that the Fourier transform of the finite difference scheme is a good approximation of the partial derivative ( ∂u ∂x ) l over the range of waves with wavelength longer than 4 h, in order to have enough grid points per wavelength (one must bear in mind that the error is a decreasing function of the number of grid points per wavelength); the a priori unknowns coefficients γ k must be choosen so as to minimize the integrated error:

E = π 2 -π 2 |λ h -λ h| 2 d(λ h) = π 2 -π 2 |λ h + j m k=-m γ k e j k ω h h| 2 d(λ h) = π 2 -π 2 |ζ + j m k=-m γ k {cos( k ζ) + j sin( k ζ)} | 2 dζ = π 2 -π 2      ζ - m k=-m γ k sin( k ζ)   2 +   m k=-m γ k cos( k ζ)   2    dζ = 2 π 2 0      ζ - m k=-m γ k sin( k ζ)   2 +   m k=-m γ k cos( k ζ)   2   F o r P e e r R e v i e w O n l y
The conditions that E is a minimum are:

∂E ∂γ i = 0 , i = -m, . . . , m (23)
i. e.:

π 2 0 -ζ sin( i ζ) + m k=-m γ k cos ( (k -i) ζ) dζ = 0 ( 24
)
Changing i into -i, and k into -k in the summation yields:

π 2 0 ζ sin( i ζ) + m k=-m γ -k cos ( (-k + i) ζ) dζ = 0 (25)
i. e.:

π 2 0 ζ sin( i ζ) + m k=-m γ -k cos ( (k -i) ζ) dζ = 0 (26)
Thus:

π 2 0 m k=-m {γ -k + γ k } cos ( (k -i) ζ) dζ = 0 (27)
which yields:

π 2 {γ -i + γ i } + m k =i, k=-m γ -k + γ k k -i sin (k -i) π 2 = 0 ( 28
)
which can be considered as a linear system of 2 m + 1 equations, the unknowns of which are the γ -i +γ i , i = -m, . . . , m. The determinant of this system is not equal to zero, while it is the case of its second member: the Cramer formulae give then, for i = -m, . . . , m:

γ -i + γ i = 0 (29)
or:

γ -i = -γ i ( 30
)
For i = 0, one of course obtains:

γ 0 = 0 (31)
All this ensures: m being a strictly positive integer, a 2m + 1-points DRP scheme ( [4]) is thus given by:

u n+1 i -u n i + c τ h m k=-m γ k u n i+k = 0 (33)
where the γ k , k ∈ {-m, m} are the coefficients of the considered scheme, and satisfy the relations (30).

At this point, we need to recall that the optimization procedure described in the above to obtain DRP schemes does not always lead to consistent schemes. The analyzis of the related consistency error, by means of the Lie group theory, can be found in [3]. As expectable, DRP schemes require as many points as possible, which explains why it still raises lots of interests from scientists.

General study of DRP schemes

The dispersion relation related to a general DRP -scheme (33) is given by:

τ h m k=-m γ k e j (k ϕ+ξω τ )-B τ +e j ξ ω τ -B τ -1 = 0 (34)
from which it comes that

j ξ ω τ = B τ -ln        1 + τ m k=-m γ k e j k ϕ h        (35)
The group velocity, i.e. the velocity with which the overall shape of the wave's amplitudes propagates through space, can be expressed as

Vg = i m k=-m i k γ k e j k ϕ ω       σ m k=-m γ k e j k ϕ c + 1       (36)
from which it comes that

∂Vg ∂ϕ = i c 2 τ     c + σ m k=-m γ k e j k ϕ   m k=-m i k 2 γ k e j k ϕ -σ   m k=-m γ k e j kϕ i k   2   σ   c + σ m k=-m γ k e j k ϕ   2 (37)
Through identification of the real and imaginary part of (37), we obtain:

γ k γ i+l k 2 sin [ (k + l)ϕ ] -k l cos [ (k + l)ϕ ] = -c m k=-m k 2 γ k sin(k ϕ) (38)
and

σ m k,l=-m γ k γ l {-cos [ (k + l) ϕ] -k l sin [ (k + l) ϕ ]} = c m k=-m k 2 γ k cos(k ϕ) (39)
Due to (30), ( 40) and ( 41) respectively become:

σ m k,l=-m γ k γ l k 2 sin [ (k + l)ϕ ] -k l cos [ (k + l)ϕ ] = -2 c m k=1 k 2 γ k sin(k ϕ) (40)
and

σ m k,l=-m γ k γ l {-cos [ (k + l) ϕ] -k l sin [ (k + l) ϕ ]} = 0 (41)
Denote by T p , p ∈ IN * the Chebyshev polynomial of the first kind, and by U p , p ∈ IN * the Chebyshev polynomial of the second kind:

cos(p x) = T p (cos(x)) = p 2 p 2 k=0 (-1) k (p -k -1)! k! (p -2 k)! (2 cos(x)) p-2k (42) sin(p x) = sin(x) U p (cos(x)) (43)
where:

U p (cos(x)) = p 2 k=0 (-1) k (p -k)! k! (p -2 k)! (2 cos(x)) p-2k ( 44
)
p 2 denotes the integer part of p 2 . Equations ( 40), (41) can thus be written as:

σ m k,l=-m γ k γ l k 2 sin(ϕ) U k+l (cos(ϕ)) -k l T k+l (cos(ϕ)) = -c m k=-m k 2 γ k sin(ϕ) U k (cos(ϕ)) (45) and σ m k,l=-m γ k γ l {T k+l (cos(ϕ)) + k l sin(ϕ) U k+l (cos(ϕ))} = 0 (46)
Using the relation:

sin(ϕ) = 1 -cos 2 (ϕ) (47)
equations ( 45), (46) can be written as: where, for all θ ∈ IR:

f 1 (θ) = σ m k,l=-m γ k γ l k 2 1 -θ 2 U k+l (θ) -k l T k+l (θ) + c m k=-m k 2 γ k 1 -θ 2 U k (θ) (50)
i.e.:

f 1 (θ) = σ m k,l=-m γ k γ l k 2 1 -θ 2 U k+l (θ) -k l T k+l (θ) + 2 c m k=1 k 2 γ k 1 -θ 2 U k (θ) (51)
and

f 2 (θ) = m k,l=-m γ k γ l T k+l (θ) + k l 1 -θ 2 U k+l (θ) (52)
Due to:

T j (1) = 1 ∀ j ∈ IN * (53)
it is worth noting that:

f 1 (1) = -σ m k,l=-m γ k γ l k l (54)
and

f 2 (1) = m k,l=-m γ k γ l ( 55
)
The knowledge of the scheme coefficients γ k , k ∈ {-m, m}, enables one to study their variations, and to determine wether the equations ( 48), (49) admit a solution. One can thus know wether ∂V g ∂ϕ = 0 admits real roots, i. e. wether the schema has spurious caustics.

Numerical application: the 3-points DRP scheme

The 3-points DRP scheme is given by:

γ 1 = 0.63662 (56)
We thus have:

f 1 (1) = -2 σ γ 2 1 -γ 2 1 = 0 (57) 9 F o r P e e r R e v i e w O n l y and f 2 (1) = 2 γ 2 1 -γ 2 1 = 0 (58)
For the 3-points DRP scheme, the dispersion relation is:

e i ϕ -e -η ω τ + e i τ ξ ω + e iτ ξ ω -0.63662 + 0.63662 e 2 i ϕ σ = 0 (59)
which leads to:

e i τ ξω = e i ϕ e -ηω τ e i ϕ + 0.63662 σ (e 2 i ϕ -1) (60)
It yields:

ξ ω = 1 τ Arctan -(1 + 0.63662 σ ) sin(ϕ) 1 + (0.63662 σ -1) cos(ϕ) (61)
The derivative ∂V g ∂ϕ of the group velocity V g vanishes for ϕ = 0, ϕ = ±π 2 , and ϕ = 0.950935. The 3-points DRP scheme thus admits spurious caustics. We now illustrate the caustic phenomenon considering the two following sinusoidal wave packets:

u 1 = e -α (x-x 1 0 -c t) 2 cos [ k 1 (x -x 1 0 -c t) ] (62)
u 2 = e -α (x-x 2 0 -c t) 2 cos [ k 2 (x -x 2 0 -c t) ] (63)
where α > 0. The two wave packets are initially centered at x 1 0 and x 2 0 , respectively. The group velocity of the two wave packets are V 1 = V g (k 1 ) and V 2 = V g (k 2 ), respectively, where the function V g (x) is associated to the numerical scheme used to solve Eq. ( 1). If the solution obeys the linear advection law given by Eq. ( 1), the initial field is uniformly advected at speed c, while, if the advection speed is scale-dependent (as in numerical solutions), the two packets will travel at different speeds, leading to the rise of discrepancies with the constant-speed solution. Another dispersive error is the shape-deformation phenomenon: due to numerical errors, the exact shape of the wave packets will not be exactly preserved. This secondary effect will not be considered below, since it is not related to the existence of spurious caustics. It is emphasized here that the occurence of spurious caustics originates in the differential error in the group velocity, not in the fact that shapes of the envelope of the wave packets are not preserved. The issue of deriving shapepreserving schemes for passive scalar advection has been adressed by several authors (e.g. [7,8]). The spurious caustics will appear if the two wave packets happen to get superimposed. During the cross-over, the L ∞ norm of the error (defined as the difference between the constant-speed solution and the dispersive one) will exhibit a maximum. The characteristic life time of the caustic, t * , depends directly on the difference between the advection speeds of the two wave packets and the wave packet widths. Denoting l 1 and l 2 the characteristic length of the two wave packets, the time during which they will be (at least partially) superimposed can be estimated as

t * = l 1 + l 2 |V 1 -V 2 | (64)
It is seen that, since caustics are defined as solutions for which ∂V g /∂k = 0, t * will be large if

|k 1 -k 2 | 1. Noting k 1 = k c + δk and k 2 = k c -δk, one obtains t * l 1 + l 2 2(δk) 2 ∂ 2 Vg ∂k (k c ) (65) leading to t * ∝ (δk) -2 .
Neglecting shape-deformation effects and assuming that the numerical scheme is non-dissipative, the numerical error E is given by:

E = | e -α(x-x 1 0 -c t) 2 cos[k 1 (x -x 1 0 -ct) ] -e -α(x-x 1 0 -t V 1 ) 2 cos[k 1 (x -x 1 0 -tV 1 )] + e -α(x-x 2 0 -c t) 2 cos[k 2 (x -x 2 0 -ct)] -e -α(x-x 2 0 -t V2) 2 cos [ k 2 (x -x 2 0 -t V 2 ) ]| (66) A simple analysis show that lim t→+∞ L ∞ (E(t)) = L ∞ (u 1 (t = 0)), max t L ∞ (E(t)) = 2L ∞ (u 1 (t = 0)) (67)
The time history of the L ∞ norm of E for the 3-points DRP scheme scheme, is displayed in Fig. 1, showing the occurence of the caustic and the sudden growth of the L ∞ error norm. As another illustration of spurious caustics phenomena, Figures 3,4 display snapshots of the computed solution at different times: one easily notices that, from t = 20 τ , the amplitude of the numerical solution begins to grow in a non-admissible way, in conjunction with the appearance of oscillations.

Concluding remarks

In the above, we have set a general method that enables one to determine wether a DRP scheme admits or not spurious caustics. The existence of spurious numerical caustics in linear advection DRP schemes has been proved. This linear dispersive phenomenon gives rise to a sudden growth of the L ∞ norm of the error, which corresponds to a local focusing of the numerical error in both space and time. In the present analysis, spurious caustics have been shown to occur in polychromatic solutions. The energy of the caustic phenomenon depends on the number of spectral modes that will get superimposed at the same time. As a consequence, the spurious error pile-up will be more pronounced in simulations with very small wave-number increments. It has been shown that a popular existing scheme, as the 3-points DRP -scheme, allows the existence of spurious caustics. Please find the revised version of our paper "Spurious caustics of Dispersion Relation Preserving schemes". We thank you for the time you must have spent on this manuscript, and for your comments, which were of great benefit, and helped us correct imprecisions. We took them into consideration, and did our best to use them. Our detailed answer is given in the following. 2. In page 3, we have specified, just after eq. ( 8) that ∂Vg ∂k reduces ∂Vg ∂ϕ since we deal with a constant speed transport equation.

Respectfully yours,

Claire

3. In page 4, we have clarified the notation for a linear finite difference scheme.

4. In page 4, we have explained that for µ = 0, the Burgers equation reduces to the linear transport advection equation, which is the object of our study.

5. In page 5, we have defined the term wavenumber, and, in page 7, the term group velocity.

6. In page 5, we have explained that the condition "over the range of waves with wavelength longer than 4h enables one to have enough grid points per wavelength, the error being a decreasing function of the number of grid points per wavelength.

7. In page 7, we have corrected the signs, and our relation effectively depends on the physical velocity c:

u n+1 i -u n i + c τ h m k=-m γ k u n i+k = 0
Also, we have recalled that the optimization procedure described in the above to obtain DRP schemes does not always lead to consistent schemes. The analyzis of the related consistency error, by means of the Lie group theory, can be found in our paper 1 .

8. As another illustration of spurious caustics phenomena, we have plotted snapshots of the computed solution at different times in order to show how the spatial profile is affected by caustics.

Minor changes

1. Page 2 (and in all the text), we have used the notation j for the complex square root of -1.

2. Page 2 (line 18): "powerful" has been replaced with "efficient".

3. Page 2: the term "monochromatic references solutions" has been explained.

4. Page 2: "gammak" has been replaced with γ k .

5. Page 3: "admits" has been replaced with "admit".

6. Page 5: "unknows" has been replaced with "unknown". 8. Page 10: "DRP scheme admits thus" has been replaced by "DRP scheme thus admits". 9. Page 9 (and elsewhere): "Cos" has been been changed into cos.

10. Page 10: "passively" has been replaced by "uniformly".

11. Page 10: "occurance" has been replaced by "occurrence".

Introduction

The analysis and the control of numerical error in discretized propagation-type equations is of major importance for both theoretical analysis and practical applications. A huge amount of works has been devoted to the analysis of the numerical errors, its dynamics and its influence on the computed solution (the reader is referred to classical books, among which [6,14,9,10]). The emergence of Dispersion-Relation-Preserving (DRP) schemes [4]), which have the same dispersion relation as the original partial difference equations, enables one to have very accurate high order finite difference schemes.

The two sources of numerical error are the dispersive and dissipative properties of the numerical scheme, which are very often investigated in unbounded or periodic domains thanks to a spectral analysis. It appears that existing works are mostly devoted to linear, one-dimensional numerical models, such as the linear advection equation

∂u ∂t + c ∂u ∂x = 0 (1)
where c is a constant uniform advection velocity. The two sources of numerical error are the dispersive and dissipative properties of the numerical scheme, which are very often investigated in unbounded or periodic domains thanks to a spectral analysis. Following this approach, a monochromatic wave is used to measure the accuracy of the scheme. Such a tool is very efficient and provides the user with a deep insight into the discretization errors. But some results coming from practical numerical experiments still remain unexplained, despite the linear character of the discrete numerical model.

As an example, let us note the sudden growth of the numerical error for long range propagation reported by Zingg [16] for a large set of numerical schemes, including optimized numerical schemes.

The usual modal analysis is almost always applied to monochromatic reference solutions, with the purpose of analyzing the error committed on both their amplitude and their phase, leading to classical plots of the relative error as the function of the Courant number and/or the number of grid points per wavelength. Therefore, dispersive phenomena associated to polychromatic solutions are usually not taken into account.

The present paper deals with the analysis of linear dispersive mechanism which results in local error focusing, i.e. to a sudden local error burst in the L ∞ norm for polychromatic solutions. This phenomena is reminiscent of the physical one referred to as the caustic phenomenon in linear dispersive physical models [15], and will be referred to as the spurious caustic phenomenon hereafter. It extends our previous work [2] to DRP schemes. The present analysis is restricted to interior stencil, and the influence of boundary conditions will not be considered. The paper is organized as follows. Main elements of caustic theory of interest for the present analysis are recalled in section 2. DRP schemes are presented in section 3. Their caustical analysis is exposed in section 4. A numerical example is presented in section 5.

Caustics

The solution of Eq. ( 1) is taken under the form:

u(x, t, k) = e j (k x-ω t) (2)
where j denotes the complex square root of -1, ω = ξ ω + i η ω is the complex phase, and k the real wave number: where λ 0 denotes the wavelength, i.e. For dispersive waves, it is recalled that the group velocity V g (k) is defined as

k = 2 π λ 0 (3
V g (k) ≡ ∂ξ ω ∂k (4)
A caustic is defined as a focusing of different rays in a single location. The equivalent condition is that the group velocity exhibits an extremum, i.e. there exists at least one wave number k c such that

∂V g ∂k (k c ) = 0 (5)
The corresponding physical interpretation is that wave packets with characteristic wave numbers close to k c will pile-up after a finite time and will remain superimposed for a long time, resulting in the existence a region of high energy followed by a region with very low fluctuation level. The linear continuous model Eq. ( 1) is not dispersive if the convection velocity c is uniform, and therefore the exact solution does not exhibits caustics since the group velocity does not depend on k. The discrete solution associated with a given numerical scheme will admit spurious caustics, and therefore spurious local energy pile-up and local sudden growth of the error, if the discrete dispersion relation is such that the condition ( 5) is satisfied. For a uniform scale-dependent convection velocity, such spurious caustics can exist in polychromatic solutions only, since they are associated to the superposition of wave packets with different characteristic wave numbers. Set:

k = ϕ σ c τ ( 6
)
where σ is the cfl number, defined as σ = cτ /h, h and τ denoting respectively the mesh size and time step.

The non-dimensional wave number is defined as ϕ = kh, where k is the wave number of the signal under consideration. The general dispersion relation associated with the discrete scheme enables us to obtain the corresponding group velocity, given by:

V g = h ∂ξ ω ∂ϕ (7)
The numerical solution will therefore admit spurious caustics if

∂V g ∂k = ∂V g ∂ϕ ∂ϕ ∂k = 0 ⇐⇒ ∂V g ∂ϕ = 0 (8)
Since we deal with a constant speed transport equation: The numerical solution will therefore admit spurious caustics if

∂V g ∂k = ∂V g ∂ϕ ∂ϕ ∂k = 0 (
∂V g ∂ϕ = 0 (11)
The corresponding values of ϕ and k will be respectively denoted ϕ c and k c . Spurious caustics are associated with characteristic lines given by

x t = U c ( 12
)
where

U c = V g (ϕ c ) ( 13
)

DRP schemes

The Burgers equation:

u t + c u u x -µ u xx = 0, (14)
c, µ being real constants, plays a crucial role in the history of wave equations. It was named after its use by Burgers [1] for studying turbulence in 1939.

For µ = 0, (14) reduces to the linear advection equation

u t + c u x = 0, (15)
which will be the object of our study.

For given natural integers i ∈ [0, n x ], n ∈ [0, n t ], a linear finite difference scheme for this equation can be written under the form:

F u mn li = 0 (16)
where:

u mn li = u (l i h, m n τ ) (17)
and where F is a linear function of the u mn li ,

l i ∈ [1, n x ], m n ∈ [0, n t ]. Common values for l i are i -1, i, i + 1, common values for m n are n -1, n, n + 1.
F depends on the mesh size h, and the time step τ .

A numerical scheme is specified by selecting an appropriate expression of F. Then, depending on them, one can obtain optimum schemes, for which the error will be minimal. m being a strictly positive integer, the first derivative ∂u ∂x is approximated at the l th node of the spatial mesh by:

( ∂u ∂x ) l m k=-m γ k u n i+k (18)
Following the method exposed by C. Tam and J. Webb in [4], the coefficients γ k are determined requiring the Fourier Transform of the finite difference scheme (18) to be a close approximation of the partial derivative ( ∂u ∂x ) l . ( 18) is a special case of:

( ∂u ∂x ) l m k=-m γ k u(x + k h) (19)
where x is a continuous variable, and can be recovered setting x = l h. Denote by ω the phase. Applying the Fourier transform, referred to by , to both sides of (19), yields:

j ω u m k=-m γ k e j k ω h u (20)
Comparing the two sides of (20) enables us to identify the wavenumber λ of the finite difference scheme (18), i.e. the number of wavelengths per unit distance, and the quantity

1 j m k=-m
γ k e j k ω h , one obtains:

λ = -j m k=-m γ k e j k ω h (21)
To ensure that the Fourier transform of the finite difference scheme is a good approximation of the partial derivative ( ∂u ∂x ) l over the range of waves with wavelength longer than 4 h, in order to have enough grid points per wavelength (one must bear in mind that the error is a decreasing function of the number of grid points per wavelength); the a priori unknowns coefficients γ k must be choosen so as to minimize the integrated error: The conditions that E is a minimum are:

E = R π 2 -π 2 |λ h -λ h| 2 d(λ h) = R π 2 -π 2 |λ h + j m X k=-m γ k e j k ω h h| 2 d(λ h) = R π 2 -π 2 |ζ + j m X k=-m γ k {cos( k ζ) + j sin( k ζ)} | 2 dζ = R π 2 -π 2 8 < : 2 4 ζ - m X k=-m γ k sin( k ζ) 3 5 2 + 2 4 m X k=-m γ k cos( k ζ) 3 5 2 9 = ; dζ = 2 R π 2 0 8 < : 2 4 ζ - m X k=-m γ k sin( k ζ) 3 5 2 + 2 4 m X k=-m γ k cos( k ζ)
∂E ∂γ i = 0 , i = -m, . . . , m (23)
i. e.:

π 2 0 -ζ sin( i ζ) + m k=-m γ k cos ( (k -i) ζ) dζ = 0 (24)
Changing i into -i, and k into -k in the summation yields:

π 2 0 ζ sin( i ζ) + m k=-m γ -k cos ( (-k + i) ζ) dζ = 0 (25)
i. e.:

π 2 0 ζ sin( i ζ) + m k=-m γ -k cos ( (k -i) ζ) dζ = 0 (26)
Thus:

π 2 0 m k=-m {γ -k + γ k } cos ( (k -i) ζ) dζ = 0 (27)
which yields:

π 2 {γ -i + γ i } + m k =i, k=-m γ -k + γ k k -i sin (k -i) π 2 = 0 (28)
which can be considered as a linear system of 2 m + 1 equations, the unknowns of which are the γ -i +γ i , i = -m, . . . , m. The determinant of this system is not equal to zero, while it is the case of its second member: the Cramer formulae give then, for i = -m, . . . , m:

γ -i + γ i = 0 (29)
or:

γ -i = -γ i (30)
For i = 0, one of course obtains:

γ 0 = 0 (31)
All this ensures: m being a strictly positive integer, a 2m + 1-points DRP scheme ( [4]) is thus given by:

u n+1 i -u n i + c τ h m X k=-m γ k u n i+k = 0 (33)
where the γ k , k ∈ {-m, m} are the coefficients of the considered scheme, and satisfy the relations (30).

At this point, we need to recall that the optimization procedure described in the above to obtain DRP schemes does not always lead to consistent schemes. The analyzis of the related consistency error, by means of the Lie group theory, can be found in [3]. As expectable, DRP schemes require as many points as possible, which explains why it still raises lots of interests from scientists.

General study of DRP schemes

The dispersion relation related to a general DRP -scheme (33) is given by:

τ h m X k=-m
γ k e j (k ϕ+ξω τ )-B τ +e j ξω τ -B τ -1 = 0 (34)

from which it comes that

j ξω τ = B τ -ln 0 B B B B B @ 1 + τ m X k=-m γ k e j k ϕ h 1 C C C C C A (35)
The group velocity, i.e. the velocity with which the overall shape of the wave's amplitudes propagates through space, can be expressed as

Vg = i m X k=-m i k γ k e j k ϕ ω 0 B B B B @ σ m X k=-m γ k e j k ϕ c + 1 1 C C C C A (36)
from which it comes that

∂Vg ∂ϕ = i c 2 τ 0 @ 0 @ c + σ m X k=-m γ k e j k ϕ 1 A m X k=-m i k 2 γ k e j k ϕ -σ 0 @ m X k=-m γ k e j kϕ i k 1 A 2 1 A σ 0 @ c + σ m X k=-m γ k e j k ϕ 1 A 2 (37)
Through identification of the real and imaginary part of (37), we obtain:

γ k γ i+l ˘k2 sin [ (k + l)ϕ ] -k l cos [ (k + l)ϕ ] ¯= -c m X k=-m k 2 γ k sin(k ϕ) (38)
and

σ m X k,l=-m γ k γ l {-cos [ (k + l) ϕ] -k l sin [ (k + l) ϕ ]} = c m X k=-m k 2 γ k cos(k ϕ) (39)
Due to (30), ( 40) and ( 41) respectively become:

σ m X k,l=-m γ k γ l ˘k2 sin [ (k + l)ϕ ] -k l cos [ (k + l)ϕ ] ¯= -2 c m X k=1 k 2 γ k sin(k ϕ) (40)
cos(p x) = T p (cos(x)) = p 2 p 2 k=0 (-1) k (p -k -1)! k! (p -2 k)! (2 cos(x)) p-2k (42)
sin(p x) = sin(x) U p (cos(x))

where:

U p (cos(x)) = p 2 k=0 (-1) k (p -k)! k! (p -2 k)! (2 cos(x)) p-2k (44)
p 2 denotes the integer part of p 2 . Equations ( 40), (41) can thus be written as:

σ m X k,l=-m γ k γ l ˘k2 sin(ϕ) U k+l (cos(ϕ)) -k l T k+l (cos(ϕ)) ¯= -c m X k=-m k 2 γ k sin(ϕ) U k (cos(ϕ)) (45)
and

σ m X k,l=-m γ k γ l {T k+l (cos(ϕ)) + k l sin(ϕ) U k+l (cos(ϕ))} = 0 (46)
Using the relation:

sin(ϕ) = 1 -cos 2 (ϕ) (47)
equations ( 45), (46) can be written as: where, for all θ ∈ IR:

f 1 (θ) = σ m X k,l=-m γ k γ l n k 2 p 1 -θ 2 U k+l (θ) -k l T k+l (θ) o + c m X k=-m k 2 γ k p 1 -θ 2 U k (θ) (50)
i.e.:

f 1 (θ) = σ m X k,l=-m γ k γ l n k 2 p 1 -θ 2 U k+l (θ) -k l T k+l (θ) o + 2 c m X k=1 k 2 γ k p 1 -θ 2 U k (θ) (51)
and

f 2 (θ) = m X k,l=-m γ k γ l n T k+l (θ) + k l p 1 -θ 2 U k+l (θ) o (52)
Due to:

T j (1) = 1 ∀ j ∈ IN * (53)
it is worth noting that:

f 1 (1) = -σ m k,l=-m γ k γ l k l (54)
and

f 2 (1) = m k,l=-m γ k γ l (55)
The knowledge of the scheme coefficients γ k , k ∈ {-m, m}, enables one to study their variations, and to determine wether the equations (48), (49) admit a solution. One can thus know wether ∂Vg ∂ϕ = 0 admits real roots, i. e. wether the schema has spurious caustics.

Numerical application: the 3-points DRP scheme

The 3-points DRP scheme is given by:

γ 1 = 0.63662 (56)
We thus have: For the 3-points DRP scheme, the dispersion relation is:

f 1 (1) = -2 σ γ 2 1 -γ 2 1 = 0 ( 57
e i ϕ
" -e -ηω τ + e i τ ξω " + e iτ ξω `-0.63662 + 0.63662 e 2 i ϕ ´σ = 0 (59) which leads to:

e i τ ξω = e i ϕ e -ηω τ e i ϕ + 0.63662 σ (e 2 i ϕ -1) (60)
It yields:

ξ ω = 1 τ Arctan -(1 + 0.63662 σ ) sin(ϕ) 1 + (0.63662 σ -1) cos(ϕ) (61)
The derivative ∂Vg ∂ϕ of the group velocity V g vanishes for ϕ = 0, ϕ = ±π 2 , and ϕ = 0.950935. The 3-points DRP scheme thus admits spurious caustics. We now illustrate the caustic phenomenon considering the two following sinusoidal wave packets:

u 1 = e -α (x-x 1 0 -c t) 2 cos [ k 1 (x -x 1 0 -c t) ] (62)
u 2 = e -α (x-x 2 0 -c t) 2 cos [ k 2 (x -x 2 0 -c t) ] (63)
where α > 0. The two wave packets are initially centered at x 1 0 and x 2 0 , respectively. The group velocity of the two wave packets are V 1 = V g (k 1 ) and V 2 = V g (k 2 ), respectively, where the function V g (x) is associated to the numerical scheme used to solve Eq. ( 1). If the solution obeys the linear advection law given by Eq. ( 1), the initial field is uniformly advected at speed c, while, if the advection speed is scale-dependent (as in numerical solutions), the two packets will travel at different speeds, leading to the rise of discrepancies with the constant-speed solution. Another dispersive error is the shape-deformation phenomenon: due to numerical errors, the exact shape of the wave packets will not be exactly preserved. This secondary effect will not be considered below, since it is not related to the existence of spurious caustics. It is emphasized here that the occurence of spurious caustics originates in the differential error in the group velocity, not in the fact that shapes of the envelope of the wave packets are not preserved. The issue of deriving shapepreserving schemes for passive scalar advection has been adressed by several authors (e.g. [7,8]). The spurious caustics will appear if the two wave packets happen to get superimposed. During the cross-over, the L ∞ norm of the error (defined as the difference between the constant-speed solution and the dispersive one) will exhibit a maximum. The characteristic life time of the caustic, t * , depends directly on the difference between the advection speeds of the two wave packets and the wave packet widths. Denoting l 1 and l 2 the characteristic length of the two wave packets, the time during which they will be (at least partially) superimposed can be estimated as

t * = l 1 + l 2 |V 1 -V 2 | (64)
It is seen that, since caustics are defined as solutions for which ∂V g /∂k = 0, t * will be large if |k 1k 2 | 1. Noting k 1 = k c + δk and k 2 = k cδk, one obtains

t * l 1 + l 2 2(δk) 2 ∂ 2 Vg ∂k (k c ) (65)
leading to t * ∝ (δk) -2 . Neglecting shape-deformation effects and assuming that the numerical scheme is non-dissipative, the numerical error E is given by: E = | e -α(x-x 1 0 -c t) 2 cos[k 1 (xx 1 0ct) ]e -α(x-x 1 0 -t V1) 2 cos[k 1 (xx 1 0 -tV 1 )] + e -α(x-x 2 0 -c t) 2 cos[k 2 (xx 2 0ct)]e -α(x-x 2 0 -t V2) 2 cos [ k 2 (xx 2 0t V 2 ) ]| (66) A simple analysis show that

lim t→+∞ L ∞ (E(t)) = L ∞ (u 1 (t = 0)), max t L ∞ (E(t)) = 2L ∞ (u 1 (t = 0)) (67)
The time history of the L ∞ norm of E for the 3-points DRP scheme scheme, is displayed in Fig. 1, showing the occurence of the caustic and the sudden growth of the L ∞ error norm. As another illustration of spurious caustics phenomena, Figures 3,4 display snapshots of the computed solution at different times: one easily notices that, from t = 20 τ , the amplitude of the numerical solution begins to grow in a non-admissible way, in conjunction with the appearance of oscillations.

Concluding remarks

In the above, we have set a general method that enables one to determine wether a DRP scheme admits or not spurious caustics. The existence of spurious numerical caustics in linear advection DRP schemes has been proved. This linear dispersive phenomenon gives rise to a sudden growth of the L ∞ norm of the error, which corresponds to a local focusing of the numerical error in both space and time. In the present analysis, spurious caustics have been shown to occur in polychromatic solutions. The energy of the caustic phenomenon depends on the number of spectral modes that will get superimposed at the same time. As a consequence, the spurious error pile-up will be more pronounced in simulations with very small wave-number increments. It has been shown that a popular existing scheme, as the 3-points DRP -scheme, allows the existence of spurious caustics.

7 Page

 7

10 Page 10

 1010 of 31 URL: http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.

Figure 1 :Figure 2 Figure 2 :

 122 Figure1: Time history (L ∞ norm) of the numerical error for the two-wave packet problem (shape deformation and dissipative errors are neglected to emphasize the linear focusing phenomenon). Numerical parameters are α = 0.0005, h = 0.01, V 1 = -2.68381, V 2 = -2.51381, corresponding to the properties of the 3-points DRP scheme, for σ = 0.9.

Figure 3 :

 3 Figure 3:The computed solution at t = 20 τ and t = 50 τ for the 3-points DRP scheme, for cf l = 0.9.

12 Page

 12

Figure 4 :

 4 Figure 4:The computed solution at t = 100 τ and t = 150 τ for the 3-points DRP scheme, for cf l = 0.9.

  Claire David, Pierre Sagaut E-mail: Claire.David@upmc.fr January 10, 2010Professor Referee for "Spurious caustics of Dispersion Relation Preserving schemes" Dear Professor Referee,

1 .

 1 David and Pierre Sagaut Page 15 of 31 URL: http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.In page 3, we have introduced the cfl number σ, and the time step, τ , instead of dt.

2 Page 7 .

 27 Page 6: it was a mystiping, we have kept just one equation.

1 Page

 1

4 Page

 4

7 Page

 7

  -m γ k γ l {cos [ (k + l) ϕ]k l sin [ (k + l) ϕ ]} = 0 (41)Denote by T p , p ∈ IN * the Chebyshev polynomial of the first kind, and by U p , p ∈ IN * the Chebyshev polynomial of the second kind:

Figure 1 :Figure 2 Figure 2 :

 122 Figure1: Time history (L ∞ norm) of the numerical error for the two-wave packet problem (shape deformation and dissipative errors are neglected to emphasize the linear focusing phenomenon). Numerical parameters are α = 0.0005, h = 0.01, V 1 = -2.68381, V 2 = -2.51381, corresponding to the properties of the 3-points DRP scheme, for σ = 0.9.

Figure 3 :

 3 Figure 3:The computed solution at t = 20 τ and t = 50 τ for the 3-points DRP scheme, for cf l = 0.9.

12 Page

 12

Figure 4 :

 4 Figure 4:The computed solution at t = 100 τ and t = 150 τ for the 3-points DRP scheme, for cf l = 0.9.