Luana Russo

Laurent Beauguitte

S Bourdin

M Bussi

B Cautrès

C Colange

S Freire-Diaz

A Jadot

J Rivière

Aggregation Level Matters: Evidences from French Electoral Data

Keywords: Abstention, Ecological Inference, Linear Regression Model, Multicollinearity, Parisian Agglomeration, Political Geography

This working paper, issued from the Cartelec project, raises a methodological issue regarding ecological inference. Studying abstention on two recent French elections in the Parisian agglomeration, a linear regression is launched on three different levels of aggregation from the finest one (polling station) until cities and electoral districts (circonscriptions). Tests on multicollinearity clearly indicate that employing the finest level of aggregation allows more variables being significant, enhances the explanatory power of the model, and improves collinearity amongst dependant variables.

Introduction

When dealing with ecological data the level of aggregation is a key point (Oumlil and Balloun 1998[10]). It is well known, at least since Robinson (1950[13]), that using different levels of aggregation leads to different results. In this paper, we will show that the finest the level of aggregation, the better the quality of the estimates and the less the multicollinearity in the model. Indeed, it was already demonstrated that, when using individual data, increasing the sample size is a viable remedy for collinearity [START_REF] Leahy | Multicollinearity: when the solution is the problem[END_REF]). In this case we do not use individual data nor a sample, but all the polling station available for a given area of the Parisian agglomeration. Hence, the number of cases is only involving the aggregation level employed, and the higher aggregation levels are only a sum of the lowest aggregation level available. This is a crucial difference with respect to the individual data, because in this case having more observation does not mean adding additional observations per se, but simply descending at a lowest aggregation level, and then indirectly having more information about the demographic, economical and social context.

In order to verify whether using different aggregation levels when dealing with ecological data leads to a difference in the final estimates, we basically use a linear regression model in which the dependent variable is the abstention and the independent ones are key variables according to both theoretical assumptions and previous empirical works (all variables are precisely described and justified below).

The paper proceeds as follows: we start with a very inclusive linear regression model that we apply only at the polling station level. Then, we test the multicollinearity and, by using a stepwise regression, we try to eliminate the multicollinearity problem. Once obtained a final, parsimonious and multicollinearity-free model, we apply the same model at the same data again, but this time we aggregate the data at two different levels: cities and electoral districts (see map 1 on the following page). We do this same procedure for two French elections: the first round of the 2007 Presidential Elections and the first round of the 2010 Regional Elections.

Since the goal of this paper is mainly methodological, we will not focus on the interpretation of the models in substantial terms. However, we try to build plausible and reliable models regarding abstention.

Method and Data

In order to verify if the level of aggregation was affecting the final estimates, we run the same linear regression model on three different aggregation levels for two elections: the 2007 French Presidential elections and the 2010 French Regional elections (first round in both cases). We choose in purpose a low (Régionales) and a high intensity election (Présidentielles) to check if the variables produce the same effect on abstention in these both different political configurations1 . As we consider the first round in both cases, the political offer can be compared without any bias as all political parties, from extreme left to extreme right, are present. The data are referring to the area around Paris called la petite couronne, which includes four departments: the 75 (Paris), the 92 (Hauts-de-Seine), the 93 (Seine-Saint-Denis), and the 94 (Val-de-Marne). These four departments are part of the Parisian agglomeration and, even if it presents strong contrasts (the 92 being the wealthiest department in France while the 93 is one of the poorest), this urban system forms in itself a relevant unit of analysis, especially for electoral and socio-economical works (Rivière, 2012[12]).

We will employ a weighted regression were the weight will be the voters who have the right to vote per each polling station (inscrits1 ). The weighted least squares regression has the advantage of taking into account the behaviour of the random errors in the model. In fact, by optimizing the weighted fitting criterion, it is possible to find the parameter estimates. Hence, the parameter estimates allow the weights to determine the contribution of each observation to the final parameter estimates.

With regard to the R-squared (R 2 ), we will report only the adjusted R 2 . In fact, while the R 2 increases when a new variable is added, the adjusted R 2 increases only if the new term improves the model more than would be expected by chance. This feature will be crucial when determining which regression model to employ.

The electoral data are provided at the polling station level by the Minister of the Interior, and in this specific case they consist in the absolute numbers of the abstention per polling station. These numbers were converted in percentages (obtained on the base of the number of members of the polling station, the inscrits) in order to make the analysis clearer, more readable and effective.

The independent variables 2 firstly regard the age of the population divided in three classes (18-24, 25-65 and more than 65 -thereafter named Young, Working age and Old). Several studies, starting with Lancelot in 1968 [START_REF]L'abstentionnisme électoral en France[END_REF], showed that old adults were more reluctant than young ones to abstain. So we expect to get a higher level of abstention where the rate of 18-24 years is higher (which of course does not mean that a relation on an aggregate level is true at the individual one).

We tested the importance of variables regarding nationality, even if only French citizens are able to vote. We considered three variables: percentage of French persons, percentage of foreigners (people living in la petite couronne but not owning French nationality), and finally percentage of immigrants (people born abroad but which now get French nationality).

Socio-economic variables used are from three different types: percentage of unemployed persons, level of education (no diploma, low level of diploma, high level of diploma1 ) and finally housing status (owners vs renters -the category free housing was introduced as a control variable). The literature regarding abstention generally highlights the relation between a high level of abstention and a low level of diploma, a high level of unemployment and the rate of renters [START_REF] Braconnier | La démocratie de l'abstention. Aux origines de la démobilisation électorale en milieu populaire[END_REF] [START_REF] Subileau | Les chemins de l'abstention : une comparaison franco-américaine[END_REF].

In all cases, the absolute number where converted in percentages obtained on the base of the total amount of the population in each polling station.

The following list sums up the dependant variables used for the different models

• Age:

young: 18-24 years old;

working age: 25-64 years old; old: 65 onward years old.

• National origin:

French; strangers;

immigrants.

• Occupational status:

unemployed.

• Education:

any diploma;

from CAP to BEP (secondary education); all types of BAC (end of the secondary education /university).

• Housing: owners; renters; free housing.

Starting from the polling station level (which is the smallest level of aggregation for electoral ecological data), we aggregated the data at other two levels: cities and electoral districts (circonscriptions législatives in French). Table 1 summarizes the number of observations for each level. As you can notice, the number of observations is different with regard the polling station and the city levels for the two elections considered. Actually, for both the elections, the real number of the polling stations was 3326, but it was necessary to exclude some cities from the analysis because some of them had a re-design of the polling stations and it was not possible to match the electoral and the sociodemographic data. As you can see, the problem was not severe, as for the Presidential Election we were able to maintain the 92.87% of the polling stations, and for the Regional Elections elections the 98.25%.

With regard to the city level, since Paris was significantly bigger than the other cities in the dataset1 , following the already administrative division, we maintained the twenty smaller areas called arrondissements. Then, in the dataset we kept each arrondissement as a separate observation.

Since we are dealing with ecological data, it is useful to underline how to interpret the dataset. Each unit of observation is a polling station, and for each variable we have a percentage. This percentage is the total amount of a given variable in a given polling station.

Findings

In order to build a reliable and effective model, we first checked the correlations (see appendix). We noticed that some of the independent variables were medium/highly correlated among them. Since this might cause some problem with the multicollinearity (Lewis-Beck 1980 [START_REF] Lewis-Beck | Applied Regression. An introduction[END_REF]; Upton and Cook 2008 [START_REF] Upton | Oxford Dictionary of Statistics[END_REF]), we will use a specific diagnostic to check the multicollinearity and try to keep it under control.

We started with a very inclusive regression model, with all the thirteen variables selected. The model employed is a weighted linear regression (in which the weights are the inscrits). Tables 2 and3 show the first model applied to both the elections considered. As you can see, this model is more effective in terms of explained variance when applied to the 2007 Presidential Election 2007. Indeed, in the 2010 Regional Election the abstention rate hits 53.64 percent1 in other words, it seems that for the 2010 election, the abstention has different dynamics that this model does not capture.

In order to check the multicollinearity for this model, we employ the "estat vif" procedure in STATA (Table 4).

As [START_REF] Hamilton | Statistics with STATA[END_REF] [START_REF] Chatterjee | Regression Analysis by Example[END_REF]) suggest that there is multicollinearity when "the largest VIF is greater than 10; or the mean VIF is considerably larger than 1."

With our largest VIF higher than 10000 and the mean VIF higher than 1600 in both the elections, we meet both criteria.

In order to fix the problem of the multicollinearity, there are two main possible solutions. We can either try to fix the multicollinearity by employing the "centring" procedure (see Hamilton 2009, p. 226), or we can decide to eliminate some variables.

To get some indications on which variables need to be eliminated, we decided to use the step-wise procedure. Table 5 shows the variables that the procedure suggests to maintain and to drop. Table 5 summarizes the results of the step-wise regressions on both the elections (PE stands for 2007 Presidential Election and RE stands for 2010 Regional Election). Let's read the step-wise regression carefully: what this procedure suggests is basically taking out two out of three variables from one category (age and housing in the case of the Presidential Election, and national origin and housing for the Regional Election).

This result was predictable. In fact, we have three different variables for each category (except for "work", which only has the variable "Unemployment"). Since these are ecological data expressed in percentages, to keep just one of the variables in a given category means that the other two variables make up the remaining proportion that the chosen variable does not cover. That is, if in a given polling station the percentage of the owners in the population is 60%, it means that the remaining 40% is composed by renters and free housing.

Hence, when deciding which variable to maintain for each category, it seems more convenient to choose either the more significant or one which was an outer.

Then, for the national origin we keep French, since in this way we'll have kind of a dichotomy French/not French (and we would not have this advantage by using one of the two other categories)1 since the immigrants are part of the French: mostly, they were born in another country and then became French citizen.. The two outer variables for "education" are "no diploma" or "BAC(s)" (end of the secondary education /university), and we decided it was indeed better to take "BAC(s)" for two reasons: 1) when looking at our regression model, the β coefficients for this variable is the highest, 2) the percentage of population with a BAC (or more) was higher than the percentage without diploma (about 38% versus 12%) 2 .

With regard to the "housing", we decided to keep the variable "owners" over "renters" and "free housing", because of the aforementioned advantage: by keeping the variable "owners", we have a net distinction in any given area between the percentage of people who own their house and the ones who do not 3 .

Finally, we decided to take all the three variables for the category "age", since in this case the middle category "Working age" is the more inclusive in terms of years, and also the most numerous. Furthermore, this is the only demographic variable we have 4 .

Tables 6 and7 show the second regression model, which is a more parsimonious one, and it includes only one variable per category, which the only exception of the "age". This model includes only seven variables, while the first one included thirteen variables. As you can notice, with six variables less, the adjusted R 2 does not drop that much. In fat it is only 0.0022 less with regard to the Presidential Election, and 0.0149 less with regard to the Regional Election. Furthermore, when applying the model to the Regional Election, all the variables are significant. When considering the Presidential Election, the only non significant variable is "owners".

In order to definitively accept this model and test it on the other two aggregation levels, we have to check whether there is still a multicollinearity problem.

Table 11 shows the multicollinearity diagnostic for the second regression model. As the Table 11 demonstrates, by selecting one variable for each category (with the exception of "age") we meet the criteria suggested by Chatterjee et al (2000 [START_REF] Chatterjee | Regression Analysis by Example[END_REF]), since even if our mean VIF is slightly larger than 1, the largest VIF is not greater than 10.

City and Electoral Districts aggregation levels

In order to test whether the aggregation level affects the quality of the analysis, we test the second regression model on the same data, but using two different aggregation levels: cities and electoral districts.

With regard to the city level, as showed in Table 1, we have 124 unities for the Presidential Election and 137 for the Regional one. For both elections, 20 of those units represent the twenty smaller areas (arrondissements) in which Paris is divided.

Tables 9 and10 show the linear regression model applied to the city level. As you can see, for both elections, the model has a high adjusted R 2 value. However, only three variables are significant for the Presidential Election ("French" and "BAC", and one variable from the "age" category). With regard to the Regional Election five variables are significant, but the adjusted R 2 is even higher. These are usually signs of the possible presence of multicollinearity. In fact, it was shown that when multicollinearity occurs, the variances are large and thus far from the true value (Morris 1982[9]; Pagel and Lunneberg 1985 [START_REF] Pagel | Empirical evaluation of ridge regression[END_REF]; Mooney and Duval 1993[8]).

Table 11 shows the diagnostic for the multicollinearity applying the regression model at the city level. The diagnostic shows that at the city level there is a problem with multicollinearity. The problem is not that severe, but in both cases the biggest VIF is larger than 10 and the mean VIF is almost doubled with respect to the polling station level.

This outcome was indeed predictable, since it is well known that a remedy to the multicollinearity problem is the increase of the number of the observations. However, this usually applies to individual data and specifically refers to the increasing of the sample. In this case, we are not dealing with a sample but with the entire population available for the area considered (the petite couronne), and the data analysed are exactly the same than at the previous level: the only difference is that they have been aggregated differently, that is at a higher level (cities).

Finally, we apply the linear regression model at the electoral district level (tables 12 and 13). As you can see, running the regression model at the electoral district level, the adjusted R 2 increases even more, as the significant variables decrease: only three per elections (the same in both elections), and with lower levels of significance. This suggests that the multicollinearity could be even increased.

Furthermore, the model looses more explanatory power. Table 14 shows the multicollinearity diagnostic at the electoral district level. With respect to the city level, the multicollinearity is increased, both with regard to the larger VIF and the mean VIF, showing that the more aggregated are the data the more the collinearity problem increases.

Conclusions

The main goal of this paper was to show that when using ecological data, the higher the aggregation level employed, the lower the explanatory power of the regression model.

In order to demonstrate this, we took two elections (first round of the 2007 Presidential Election and first round of the 2010 Regional Election) and we applied the same model at the petite couronne (Departments 75, 92, 93, 94), on the same data aggregated at three different levels.

The higher levels (cities and electoral districts) are obtained by aggregating the polling station data. Upon the successful design of a linear regression model at the polling station level (high explanatory power and absence of multicollinearity), we applied that same model on the city and electoral district levels.

We found out the the higher the aggregation level, the more the multicollinearity and the less the explanatory power of the model. Hence, we can conclude that the aggregation level has an impact on the quality of the estimates and on the explanatory power of the model.

Figure 1 :

 1 Figure 1: Cities -top -and polling stations in la petite couronne

Table 1 :

 1 Number of observations per aggregation level.

		Polling stations Cities Electoral districts
	Presidential Election 2007	3089	124	59
	Regional Election 2010	3268	136	59

Table 2 :

 2 First model: observations and adjusted R 2 .

	Election	Obs. Adj. R 2
	Presidential Election 2007 3089	0.5011
	Regional Election 2010	3268	0.2645

Table 3 :

 3 First model: standardized β coefficients.

	Variables	Std. β coeff.	Std. Err.
	(ID: Abstention)	PE	RE	PE	RE
	Young	-.0261346 -.0617661**	.0236746 .0831016
	Working age	-.0080677 -.3182642***	.0219815 .0794442
	Old	.0935859** -.2617109***	.0195017 .0696243
	French	.0935859	-.4335976	.6670382 2.403678
	Strangers	.-.6190991	-.2436914	.6675262 2.405822
	Immigrants	-.0302352	-.0581329	.0242692 .0868779
	Unemployed	.1655555***	.0841649*	.0417229 .1514527
	No diploma	.1507885*	.2397054**	.0280403 .1006813
	CAP/BEP	.1150329* .3921551***	.0354464 .1275432
	BAC (all types)	.1150329*	.3941463**	.0226946 .0818271
	Owners	.1438963	.0389029	.0136184 .0490349
	Renters	.1735346*	.0319292	.0138451 .0497786
	Free housing	0528821*	0585268*	.0213972 .0760094
	***= t ≥ 4; **= t ≥ 3; *= t ≥ 1.96
	weights: inscrits

Table 4 :

 4 explains: "1/VIF (variance inflation factor) (or 1/R 2 ) tells us what proportion of an x variable's variance is independent of all the other x variables. A low proportions [...] indicates potential trouble". First model: multicollinearity

	Variables	VIF		1/VIF
		PE	RE	PE	RE
	Strangers	10844.34 11417.53 0.000092 0.000088
	French	10819.49 11391.43 0.000092 0.000088
	BAC (all types)	46.42	47.57 0.021540 0.021020
	Renters	43.16	43.35 0.023168 0.023068
	Owners	41.69	42.01 0.023986 0.023803
	Immigrants	18.30	19.33 0.054636 0.051746
	No diploma	17.82	18.59 0.056113 0.053798
	CAP/BEP	11.21	11.17 0.089190 0.089527
	Working age	4.69	4.74 0.213166 0.210754
	Unemployed	4.37	4.54 0.229086 0.220261
	Old	3.70	3.67 0.270191 0.272833
	Free housing	2.40	2.34 0.416697 0.426763
	Young	1.66	1.70 0.600729 0.588033
	Mean VIF	1681.48 1769.84
	Chatterjee et al. (2000

Table 5 :

 5 Step-wise regression

	Removing.	PE RE	YES NO	YES NO	NO NO	NO YES	YES NO	YES YES	NO NO	NO NO	NO NO	NO NO	NO YES	NO YES	NO NO
	P≥t	PE RE	0.050 0.000	0.050 0.000	0.000 0.000	0.050 0.791	0.050 0.000	0.050 0.372	0.000 0.000	0.000 0.000	0.000 0.000	0.000 0.000	0.038 0.050	0.010 0.050	0.003 0.002
	Std. Err.	PE RE	-.06564	-.0683442	.0130077 .0589315	.0138716 -	-.0416564	--	.0418275 .1157916	.0244067 .0870589	.0256073 .1069413	.0134074 .0703127	.0128846 -	.0132596 -	.0214196 0505448
	Variables Coefficients	(ID: Abstention) PE RE	Young --.246053	Working age --.7187066	Old .0813572 -.6389446	French -.0699387 -	Strangers -.1523258	Immigrants --	Unemployed .2550422 .4125384	No diploma .0945562 .4533225	CAP/BEP .099413 .7758603	BAC (all types) -.0791843 .2587817	Owners .0267003 -	Renters .034293 -	Free housing .0632967 .1535383	PE: obs = 3089; Adj R 2 = 0.5136	RE: obs = 3268; Adj R 2 = 0.3530

Table 6 :

 6 Second model: observations and adjusted R 2

	Election	Obs Adj R 2
	Presidential Election 2007 3089	0.4989
	Regional Election 2010	3268	0.2496

Table 7 :

 7 Second model: standardized β coefficients and standard error

	Variables	Coefficients	Std. Err.
	(ID: Abstention)	PE	RE	PE	RE
	Young	-.0424543* -.0685306**	.0216045 .0753659
	Working age	.0409709* -.1482356***	.0158999 .0580375
	Old	.131871*** -.1428987***	.0154836 .0556107
	French	-.1524588*** -.0910477**	.0097089 .0346454
	Unemployed	.1834103***	.0649041*	.0383116 .1399367
	BAC (all types)	-.5578792*** -.2789387***	.0063275 .0228775
	Owners	-.017734	.0400247*	.0027952 .0101924

***= t ≥ 4; **= t ≥ 3; *= t ≥ 1.96 weights: inscrits

Table 8 :

 8 Second model: multicollinearity

	Variables	VIF	1/VIF
		PE RE	PE	RE
	Unemployed	3.68 3.84 0.271595 0.260290
	BAC (all types) 3.51 3.54 0.284638 0.282233
	French	2.48 2.56 0.403516 0.390102
	Working age	2.43 2.46 0.411882 0.406643
	Old	2.31 2.27 0.433036 0.440716
	Owners	1.69 1.71 0.591151 0.585970
	Young	1.39 1.39 0.717100 0.719406
	Mean VIF	2.50 2.54

Table 9 :

 9 Second model, cities: observations and adjusted R 2

	Election	Obs Adj R 2
	Presidential Election 2007 124	0.7000
	Regional Election 2010	136	0.8003

Table 10 :

 10 City level: standardized β coefficients and standard error.

	Variables	Coefficients	Std. Err.
	(ID: Abstention)	PE	RE	PE	RE
	Young	-.1205111* -.1508704*	.1064167 .178327
	Working age	-.0209643 -.3545781**	.0838745 .1406654
	Old	.0609746 -.2686493**	.0962709 .1584229
	French	-.313109* -.2835766**	.054503 .0902889
	Unemployed	.0522401	.1456201	.2432053 .4178619
	BAC (all types)	-.6178728**	-.237419*	.0329872 .0551368
	Owners	-.0572686	.044503	.0161465 .0278532

***= t ≥ 4; **= t ≥ 3; *= t ≥ 1.96 weights: inscrits

Table 11 :

 11 City level: multicollinearity

	Variables	VIF		1/VIF
		PE	RE	PE	RE
	Unemployed	10.49 11.30 0.095365 0.088493
	BAC (all types) 6.74 5.33 0.148320 0.187692
	French	5.35 5.50 0.186881 0.181974
	Working age	3.39 2.61 0.294693 0.382832
	Owners	3.36 3.49 0.297212 0.286929
	Old	3.11 2.70 0.321340 0.370344
	Young	1.50 1.46 0.668535 0.684407
	Mean VIF	4.85 4.63

Table 12 :

 12 Second model, electoral districts: observations and adjusted R 2

	Election	Obs Adj R 2
	Presidential Election 2007	59	0.8240
	Regional Election 2010	59	0.8490

Table 13 :

 13 Electoral districts level: standardized β coefficients and standard error

	Variables	Coefficients	Std. Err.
	(ID: Abstention)	PE	RE	PE	RE
	Young	-.1522416* -.1471945*	.1243341 .2556732
	Working age	-.0130935 -.2113426	.1115796 .235215
	Old	.1017578 -.1214201	.1285218 .2628305
	French	-.3678302* -.3531195*	.0669929 .1368121
	Unemployed	.0904198	.0671718	.293413 .6330592
	BAC (all types)	-.6500518** -.4319581*	.0426144 .0898281
	Owners	-.0520452	.0376641	.0193139 .0417978

***= t ≥ 4; **= t ≥ 3; *= t ≥ 1.96 weights: inscrits

Table 14 :

 14 Electoral district level: multicollinearity

	Variables	VIF		1/VIF
		PE	RE	PE	RE
	BAC (all types) 14.32 15.40 0.069812 0.064933
	Unemployed	13.07 14.95 0.076504 0.066888
	Old	7.34 6.96 0.136323 0.143775
	Working age	7.25 7.64 0.137885 0.130924
	French	6.73 6.97 0.148630 0.143407
	Owners	2.71 2.84 0.368551 0.351632
	Young	1.73 1.70 0.578233 0.589295
	Mean VIF	7.59 8.07

Table 15 :

 15 Correlation table PE 2007 (obs=3089)

	Free	housing
	BAC Owners Renters
	Old French Strangers Immigrants Unemployed No diploma CAP/BEP
	Young Working	age	1.00
			Young	Working age

Table 16 :

 16 Correlation table REG 2010 (obs=3268)

	Free	housing
	BAC Owners Renters
	Old French Strangers Immigrants Unemployed No diploma CAP/BEP
	Young Working	age	1.00	-0.16 1.00
			Young	Working age	Old

The abstention, measured in percentage of the inscrits was 16.22% (national level) in

vs 53.64 in 2010. For the space studied here, the gap is even higher (14% in 2007 and 55.79 in 2010).

Being inscrit on electoral list is mandatory to be able to vote. The inscription is theoretically an obligation for all French citizens, but the constraint has never being used, and the non inscription rate is estimated around 7 to 9% of the French adult population(Braconnier and Dormagen,

2007a[START_REF] Braconnier | Ségrégation sociale et ségrégation politique, Non inscrits, mal inscrits et abstentionnistes[END_REF]). 2 Source: INSEE National census, 2007.

For people knowing the French educational system: low level means diploma from CAP to BEP and high level -the name being a little abusive for the youngest generations -means Baccalauréat and above.

In

2008, Paris counts 2 590 000 inhabitants and the second largest city of la petite couronne is Boulogne-Billancourt with 109 000 inhabitants only. The largest Parisian arrondissement regarding population is the 15 th with 238 000 persons.

In France, if we consider all elections on a national scale, excluding referendums, the highest rate of abstention was in the

Européennes election with 59.37% of the inscrits.

The percentage of French in the population is 84.84, the strangers are the 15.16%, and the immigrants the

20.48%. The total is more than 100%[START_REF] Braconnier | La démocratie de l'abstention. Aux origines de la démobilisation électorale en milieu populaire[END_REF] The percentage of no diploma in the population is 13.01, the percentage of people holding a title form a CAP to a BEP are 18.91, and the ones olding a BAC or more are the 38.20%. All the remaining population out of this percentages is still studying.[START_REF] Chatterjee | Regression Analysis by Example[END_REF] The percentage of owners in the population is 17.45, the renters are 25.73%, and the people benefiting of the free housing the 1.81%.[START_REF] Hamilton | Statistics with STATA[END_REF] The percentage of young in the population is 9.86, the people in their working age are the 55.60%, and the seniors, that is the old the 12.84%.


Alexey Ozerov

Mathieu Lagrange

Emmanuel Vincent

Uncertainty-based learning of Gaussian mixture models from noisy data

come


Christophe Chalons
email: christophe_chalons@ljll.univ-paris-diderot.fr.

Mathieu Girardin
email: mathieu.girardin@cea.fr.

Samuel Kokh
email: samuel.kokh@cea.fr.

Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms

We propose a large time step and asymptotic preserving scheme for the gas dynamics equations with external forces and friction terms. By asymptotic preserving, we mean that the numerical scheme is able to reproduce at the discrete level the parabolic-type asymptotic behaviour satisfied by the continuous equations. By large time-step, we mean that the scheme is stable under a CFL stability condition driven by the (slow) material waves, and not by the (fast) acoustic waves as it is customary in Godunov-type schemes. Numerical evidences are proposed and show a gain of several orders of magnitude in both accuracy and efficiency.

Introduction

We are interested in the simulation of subsonic compressible flows in a specific regime where the (main) driving phenomena are stiff source terms and material transport. More precisely, we consider the system of gas dynamics with external body forces and friction. Such flow configuration may be encountered in several industrial processes like the flows involved within the core of a nuclear power plant.

We propose here a method that fulfils our task by addressing three issues. First we require our method to enable the use of large time steps in order to avoid classic Courant-Friedrichs-Lewy (CFL) restriction based on the (fast) acoustic waves of the model. Second, we want our method to accurately approximate (slow) waves that account for material transport. Third, the discretization of the stiff source terms with large time steps may severely affect the accuracy of the method. We propose to overcome this difficulty by imposing an even stronger property on our numerical scheme. Indeed, when one considers the asymptotic regime obtained for both long time and large friction coefficients, the solution of the system is formally expected to behave like the solution of a typical parabolic system. We aim at deriving a scheme that preserves this property for the discrete approximate of the solution. Such property is referred to as an asymptotic preserving (AP) property. Since its introduction in the pioneer work of [START_REF] Jin | Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms[END_REF][START_REF] Jin | Numerical Schemes for Hyperbolic Conservation Laws with Stiff Relaxation Terms[END_REF], the notion of AP numerical schemes has been investigated and implemented in the past years in a wide range of context stemming from hydrodynamics with radiative transfer [START_REF] Gosse | Asymptotic-preserving and well-balanced schemes for radiative transfer and the Rosseland approximation[END_REF][START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF][START_REF] Buet | An asymptotic preserving scheme for hydrodynamics radiative transfer models: numerics for radiative transfer[END_REF][START_REF] Berthon | An HLLC scheme to solve the M1 Model of radiative transfer in two space dimensions[END_REF][START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF], multiscale kinetics [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF], diffusive limit of the transport equation [START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF][START_REF] Naldi | Numerical schemes for hyperbolic systems of conservation laws with stiff diffusive relaxation[END_REF] to problems similar to our model [START_REF] Bouchut | Upwinding of the source term at interfaces for Euler equations with high friction[END_REF][START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF][START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF]. Without being exhaustive, we also refer the reader to [START_REF] Caflisch | Uniformly accurate schemes for hyperbolic systems with relaxation[END_REF][START_REF] Dumbser | Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws[END_REF][START_REF] Berthon | Late-time relaxation limits of nonlinear hyperbolic systems. A general framework[END_REF][START_REF] Abgrall | Special issue : Numerical approximations of hyperbolic systems with source terms and applications[END_REF][START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF].

Classic means to fulfil our first requirement for avoiding CFL based on the acoustic waves, consists in deriving an implicit in time discretization. Unfortunately this usually induces more numerical diffusion, including for the approximation of the material waves. In order to meet both first and second need, we propose a mixed implicit-explicit strategy: the terms responsible for the acoustic waves receive a time implicit treatment while the ones responsible for the transport waves are treated by an explicit update. This task is achieved by means of a Lagrange-Projection [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] algorithm as in [START_REF] Coquel | Entropy-satisfying relaxation method with large timesteps for Euler IBVPs[END_REF]. This approach provides a natural decoupling of the acoustic waves and the material waves. An approximation based on a relaxation strategy [START_REF] Suliciu | Ont the thermodynamics of fluids with relaxation and phase transitions[END_REF][START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF][START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF][START_REF] Chalons | Navier-stokes equations with several independant pressure laws and explicit predictor-corrector schemes[END_REF][START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF][START_REF] Bouchut | A reduced stability condition for nonlinear relaxation to conservative laws[END_REF] provides a simple mean to circumvent the nonlinearities involved with the equation of state of the fluid.

The effects of both gravity and friction source terms are then incorporated into the solver thanks to the concept of simple approximate Riemann solver and consistency with the integral form introduced by Gallice [START_REF] Gallice | Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source[END_REF][START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF]. This powerful method allows to account for both source terms and convective fluxes at the same time. It has been applied in [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF] for deriving an explicit scheme for the same flow model. We shall show that this approach can be followed within a mixed implicit-explicit framework for designing an AP scheme.

The outline of the paper is as follows. In the next section, we give the model under consideration and its parabolic-type asymptotic limit. In section 3, we first briefly recall the Lagrange-Projection decomposition and the pressure relaxation strategy. We then recall the concepts of simple approximate Riemann solver and consistency in the integral sense in section 4. The last part of this section gives the explicit in time numerical scheme and section 5 focuses on the Lagrangian system. At last, the implicit in time scheme for the Lagrangian system is given in section 6 and the overall mixed implicit-explicit scheme is described in section 7. Finally, section 8 gives the main properties of the mixed implicit-explicit scheme and the last section provides some numerical illustrations.

Governing equations and asymptotic behaviour

The gas dynamics equations with gravity and friction terms in Eulerian coordinates are given by

   ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p) = ρ(g -αu), ∂ t (ρE) + ∂ x ((ρE + p)u) = ρu(g -αu), (1)
where ρ, u and E denote the density, the velocity and the total energy of the fluid, g the gravitational acceleration and α the friction parameter. The pressure law p = p(ρ, e) is assumed to be a given function of the density ρ and the internal energy e defined by e = E -u 2 2 , satisfying the usual Weyl assumptions [START_REF] Weyl | Shock waves in arbitrary fluids[END_REF]. Under these assumptions and when the source terms are omitted, ( 1) is shown to be strictly hyperbolic over the phase space Ω given by Ω = {(ρ, ρu, ρE)

T ∈ R 3 , ρ > 0, e > 0},
with eigenvalues given by

λ 1 = u -c < λ 2 = u < λ 3 = u + c,
where c = [(∂p/∂ρ) + (p/ρ 2 )(∂p/∂ ) ρ ] 1/2 is the sound speed. Moreover, the characteristic fields associated with λ 1 and λ 3 are genuinely non linear while the characteristic field associated with λ 2 is linearly degenerate. Let us also recall that λ 1 and λ 3 give rise to the so-called acoustic waves, while λ 2 is associated to the transport phenomenon. We refer for instance the reader to [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] for more details. Let s = s(ρ, e) be the strictly convex mathematical specific entropy which satisfies

-T ds = de + pd 1 ρ , (2)
where T > 0 is the temperature. We obtain for a smooth solution of (1) the following equation

∂ t (ρs) + ∂ x (ρsu) = 0.
We are particularly interested in studying the long time behavior of (1) when the friction parameter goes to infinity. Let be a small positive parameter. We model this flow regime by replacing α with α/ with a slight abuse of notation and by performing the change of variable t = t in the system [START_REF] Abgrall | Special issue : Numerical approximations of hyperbolic systems with source terms and applications[END_REF]. We obtain the system

         ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p) = ρ g - α u , ∂ t (ρE) + ∂ x ((ρE + p)u) = ρu g - α u . (3)
Let us assume that the velocity u admits an asymptotic expansion in powers of of the following form

u = u 0 + u 1 + O( 2 ).
Multiplying the second equation by and letting go to 0 first gives u 0 = 0. Then inserting u = u 1 + O( 2 ) in the first equation, dividing by and letting go to 0 gives

∂ t ρ + ∂ x ρu 1 = 0.
If we now insert u = u 1 + O( 2 ) in the second equation and let go to 0, we get

∂ x p = ρg -ραu 1 .
At last, inserting u = u 1 + O( 2 ) in the third equation, dividing by and letting go to 0 gives

∂ t (ρe) + ∂ x (ρeu 1 + pu 1 ) = ρu 1 (g -αu 1 ).
The long time behaviour of the solutions of (1) for large friction coefficients is then given by the following system of partial differential equations

   ∂ t ρ + ∂ x (ρu 1 ) = 0, ∂ x p = ρ(g -αu 1 ), ∂ t (ρe) + ∂ x ((ρe + p)u 1 ) = ρu 1 (g -αu 1 ), (4)
where we note that in comparison to (1) the flow speed u has been replaced by its first order corrector u 1 in both the mass flux of the first equation and the friction term of the second equation. This observation will play a crucial role in the forthcoming developments. In addition, we note that this model has to be understood as the diffusive or parabolic limit of the hyperbolic model (1) since the second derivative of the pressure p naturally appears in the first equation of the limit system using its second equation :

∂ t ρ + ∂ x ρg -∂ x p α = 0.
From a numerical point of view and as already said in the previous section, one of our objectives is to preserve this asymptotic behaviour at the discrete level. In other words, we aim at proposing a consistent numerical scheme for (1) leading to a consistent numerical scheme for (4) when goes to zero and up to the expected changes of variables. Before proceeding and to conclude this section, let us formally rephrase this property in terms of limits with respect to the small parameter and to the time and space steps ∆t, ∆x used in the numerical approximation. Let us denote M the initial model (1), M 0 the limit model (4), S ∆t,∆x a consistent numerical scheme for (1) and S 0 ∆t,∆x its asymptotic limit. Recall that consistency of S ∆t,∆x means that lim ∆t,∆x→0

S ∆t,∆x = M , for all > 0. By definition and with a little abuse in the notations, S ∆t,∆x is said to be asymptotic preserving if S 0 ∆t,∆x is consistent with M 0 , that is lim

∆t,∆x→0 S 0 ∆t,∆x = M 0 , or equivalently lim ∆t,∆x→0 lim →0 S ∆t,∆x = lim →0 M .
In other words, the asymptotic preserving property is formally equivalent to the following order of limits interchange property : lim

∆t,∆x→0 lim →0 S ∆t,∆x = lim →0 lim ∆t,∆x→0
S ∆t,∆x .

From a practical point of view this equality formally means that for large friction coefficients, or equivalently for small values of , an asymptotic preserving scheme is expected to give good numerical results even for reasonable mesh sizes (with respect to ). This will be observed in the last section devoted to the numerical experiments.

3 Lagrange-Projection approach and relaxation procedure

In this section we briefly recall the so-called Lagrange-Projection strategy applied to [START_REF] Abgrall | Special issue : Numerical approximations of hyperbolic systems with source terms and applications[END_REF] and propose a relaxation procedure for approximating the solutions of the underlying Lagrangian system. As motivated in the introduction, these are two key ingredients of the method we propose, together with the notion of consistency in the integral sense that will be recalled in the next section. Let us begin with the Lagrange-Projection decomposition.

Lagrange-Projection decomposition

We describe here a procedure that allows to approximate the evolution of the system (1) over a time interval [t 0 , t 0 +∆t]. The guideline of the method consists in decoupling the terms responsible for the acoustic waves and the transport waves. By using the chain rule for the space derivatives we split up the operators of system [START_REF] Abgrall | Special issue : Numerical approximations of hyperbolic systems with source terms and applications[END_REF] and obtain two subsystems. The first subsystem describes the transport process and reads

   ∂ t ρ + u∂ x ρ = 0, ∂ t (ρu) + u∂ x (ρu) = 0, ∂ t (ρE) + u∂ x (ρE) = 0. (5)
The second subsystem accounts for acoustic, gravity and friction effects, namely

   ∂ t ρ + ρ∂ x u = 0, ∂ t (ρu) + ρu∂ x u + ∂ x p = ρ(g -αu), ∂ t (ρE) + ρE∂ x u + ∂ x (pu) = ρu(g -αu).
If we note τ = 1 ρ the specific volume, the above system also reads

     ∂ t τ -τ ∂ x u = 0, ∂ t u + τ ∂ x p = g -αu, ∂ t E + τ ∂ x (pu) = u(g -αu). (6)
Then for t ∈ [t 0 , t 0 + ∆t], we propose to approximate τ (x, t)∂ x • by τ (x, t 0 )∂ x • in [START_REF] Bouchut | Upwinding of the source term at interfaces for Euler equations with high friction[END_REF]. If one introduces the mass variable m defined by dm = τ (x, t 0 ) -1 dx, we obtain

   ∂ t τ -∂ m u = 0, ∂ t u + ∂ m p = g -αu, ∂ t E + ∂ m (pu) = u(g -αu), (7)
that will be referred to as the Lagrangian system. Let us note that system [START_REF] Buet | An asymptotic preserving scheme for hydrodynamics radiative transfer models: numerics for radiative transfer[END_REF] is consistent the usual form of the gas dynamics equations in Lagrangian coordinates with friction and gravity terms. We refer for instance the reader to [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] for more details. It is worth noticing that [START_REF] Buet | An asymptotic preserving scheme for hydrodynamics radiative transfer models: numerics for radiative transfer[END_REF] is easily shown to be hyperbolic over the phase space Ω Lag given by

Ω Lag = {(τ, u, E) T ∈ R 3 , τ > 0, e > 0},
with eigenvalues given by λ Lag is linearly degenerate. Importantly, we note that the sound speed only appears in the characteristic speeds of this Lagrangian model. The flow speed u is no longer present but is on the other hand the unique characteristic velocity of the first subsystem [START_REF] Bouchut | A reduced stability condition for nonlinear relaxation to conservative laws[END_REF]. System (7) (respectively (5)) then clearly encompasses to the so-called acoustic (resp. material or transport) waves. From a numerical point of view, the scheme associated with this decomposition simply consists of an usual two-step splitting strategy where [START_REF] Buet | An asymptotic preserving scheme for hydrodynamics radiative transfer models: numerics for radiative transfer[END_REF] is solve in the first step and (5) in the second one. Recall that the idea will be here to propose a time implicit treatment of the Lagrangian system [START_REF] Buet | An asymptotic preserving scheme for hydrodynamics radiative transfer models: numerics for radiative transfer[END_REF] to avoid a too restrictive CFL condition involving the sound speed c and an explicit treatment of the transport part [START_REF] Bouchut | A reduced stability condition for nonlinear relaxation to conservative laws[END_REF] to keep accuracy on the material waves. A key ingredient to get a low cost implicit scheme for the Lagrangian system will rely on a relevant pressure relaxation approximation described in the next section.

Relaxation approximation

We propose in this section a relaxation approximation of the Lagragian system [START_REF] Buet | An asymptotic preserving scheme for hydrodynamics radiative transfer models: numerics for radiative transfer[END_REF]. The main objective is to overcome the non linearities that make difficult the resolution of this system. From a numerical point of view, this strategy will be used to design a low cost time implicit treatment. We first consider the model neglecting the source terms and then extend the approach to the full model.

Relaxation approximation of the homogeneous model

The design principle of the so-called pressure relaxation methods is to introduce a larger system than the original one but easier to solve. More precisely, the objective is to discard the non linearities induced by the pressure law p = p(ρ, e). Such a strategy is now well known in the literature and we refer for instance the reader to [START_REF] Suliciu | Ont the thermodynamics of fluids with relaxation and phase transitions[END_REF][START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF][START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF][START_REF] Chalons | Navier-stokes equations with several independant pressure laws and explicit predictor-corrector schemes[END_REF][START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF][START_REF] Bouchut | A reduced stability condition for nonlinear relaxation to conservative laws[END_REF] and the references therein. To do so, we introduce a new variable Π that can be seen as a linearization of the pressure p and that is considered as a new unknown. In particular, it evolves according to its own partial differential equation. More precisely, we propose the following relaxation system for [START_REF] Buet | An asymptotic preserving scheme for hydrodynamics radiative transfer models: numerics for radiative transfer[END_REF] when the source terms are omitted :

       ∂ t τ -∂ m u = 0, ∂ t u + ∂ m Π = 0, ∂ t Π + a 2 ∂ m u = λ(p -Π), ∂ t E + ∂ m (Πu) = 0, ( 8
)
where a is a constant to be precised and λ the relaxation parameter. At least formally, we observe that in the asymptotic regime λ → +∞ we have Π → p and we recover the initial system (7) without the source terms. In order to prevent this relaxation procedure form instabilities, it is now well established that a must be chosen sufficiently large and according to the subcharacteristic condition

a > max (ρc) , (9)
for all the states under consideration (see for instance [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF] for a rigorous proof). In addition, it can be easily proved that ( 8) with λ = 0 is strictly hyperbolic with three eigenvalues given by -a, 0 and a which are nothing but approximations of the exact eigenvalues -ρc, 0 and ρc for system [START_REF] Buet | An asymptotic preserving scheme for hydrodynamics radiative transfer models: numerics for radiative transfer[END_REF]. Then and in particular, the subcharacteristic condition means that information propagates faster in the relaxation model. More importantly, the characteristic fields associated with these new eigenvalues are shown to be linearly degenerate. This property allows to solve analytically the Riemann problem associated with ( 8) when λ = 0, that is when considering an initial data made of two constant states separated by an initial discontinuity. This property justifies by itself the introduction of the proposed relaxation model and its simplicity. If we go further into the details and as it is customary, the exact Riemann solutions are self-similar and made of three contact discontinuities propagating with velocities -a, a and 0 and separating two intermediates states. From a numerical point of view, the numerical strategy for approximating the solutions of ( 7) using (8) consists in first solving [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF] with λ = 0, that is [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF] and then to take into account the source term

       ∂ t τ -∂ m u = 0, ∂ t u + ∂ m Π = 0, ∂ t Π + a 2 ∂ m u = 0, ∂ t E + ∂ m (Πu) = 0,
       ∂ t τ = 0, ∂ t u = 0, ∂ t Π = λ(p -Π), ∂ t E = 0,
in the asymptotic regime λ → +∞. Which amounts to set Π = p(ρ, E), before solving again [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF] as the time goes on. The new variable Π is said to be at equilibrium.

To conclude this section, notice that the self-similar Riemann solutions associated with [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF] being explicitly known, it is natural to use an exact Godunov scheme to numerically solve [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF]. See for instance the references above for more details.

Relaxation approximation of the Lagrangian system with source terms

Let a be a real parameter chosen in agreement with the subcharacteristic condition [START_REF] Caflisch | Uniformly accurate schemes for hyperbolic systems with relaxation[END_REF]. In order to approximate the solution of ( 7), we propose to supplement [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF] with the friction and gravity terms so that we now have to solve

       ∂ t τ -∂ m u = 0, ∂ t u + ∂ m Π = g -αu, ∂ t Π + a 2 ∂ m u = 0, ∂ t E + ∂ m (Πu) = u(g -αu). (11)
Let us underline that the solutions of the Riemann problem associated with [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF] are neither self-similar nor explicitly known anymore. This makes the use of the exact Godunov method quite a complex task. In what follows, we decide nevertheless to approximate the non self-similar Riemann solutions to [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF] by self-similar approximate Riemann solutions and to use an approximate Godunov-type method for solving [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF]. In order to guarantee the consistency of the proposed self-similar approximate Riemann solutions to [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF] with the exact ones, we will impose a generalized notion of consistency in the integral sense due to Gallice [START_REF] Gallice | Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source[END_REF][START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF] and adapted to systems with source terms. This is recalled in the next section.

To conclude this section, let us observe that ( 11) can be given the following equivalent form

       ∂ t τ -∂ m u = 0, ∂ t - → w + a∂ m - → w = a(g -αu), ∂ t ← - w -a∂ m ← - w = -a(g -αu), ∂ t E + ∂ m (Πu) = u(g -αu), (12)
where the new variables -→ w and ←w are defined by

- → w = Π + au, ← - w = Π -au. (13)
These quantities are nothing but the strong Riemann invariants associated with the characteristic speeds ±a of the relaxation system [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF] when the source terms are omitted. The closure relations for [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF] are naturally given by

u = - → w -← - w 2a , Π = - → w + ← - w 2 . (14)
This new formulation will be used hereafter to define the proposed implicit in time numerical strategy.

4 Consistency in the integral sense and explicit in time Godunovtype scheme

We briefly recall in this section the notion of consistency in the integral sense of a self-similar approximate Riemann solver for a given set of hyperbolic equations with source terms that we write in the following condensed form

∂ t U + ∂ x F(U) = S(U), (15)
supplemented with the validity of an entropy inequality

∂ t η + ∂ x q ≤ 0, (16)
where (η, q) is a strictly convex entropy-entropy flux pair. We also derive the corresponding explicit in time Godunov-type scheme for approximating the solutions to [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] and refer to [START_REF] Gallice | Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source[END_REF][START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF] for the details.

Solving the Riemann problem amounts to find the solution to [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] with the following piecewise constant initial data

U(x, t = 0) = U L if x < 0, U R if x > 0,
for any given U L and U R in the phase space. Unlike the homogeneous case corresponding to the choice S(U) = 0, the exact Riemann solution that we denote U(x, t; U L , U R ) is not self-similar. Notice however that an approximate Riemann solver W( x t ; U L , U R ) may be self-similar as in the homogeneous case provided that some consistency relations are imposed. More precisely, let us consider a simple approximate Riemann solver

U l-1 U l U l+1 = U R λ l λ l-1 λ l-2 λ 1 λ 2 ∆x L ∆x R ∆t U 2 U 1 = U L x t
W x t ; U L , U R =                U 1 = U L , x t < λ 1 , . . . U k , λ k-1 < x t < λ k , . . . U l+1 = U R , x t > λ l . (17)
From Gallice [START_REF] Gallice | Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source[END_REF][START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF], if ∆x = 1 2 (∆x L + ∆x R ) with ∆x L > 0, ∆x R > 0 and ∆t > 0 are respectively space and time steps that verify the CFL condition

max 1≤k≤l |λ k | ∆t min(∆x L , ∆x R ) ≤ 1 2 , ( 18
)
the approximate Riemann solver is said to be consistent with the integral form of (15) over the interval

[ -∆x L 2 , ∆x R 2 ]
if the integral of ( 17) approximates correctly the integral of the exact solution in the sense that there exists a function S such that

F(U R ) -F(U L ) -∆x S(∆x, ∆t; U L , U R ) = l k=1 λ k (U k+1 -U k ), (19)
where S(∆x, ∆t; U L , U R ) is consistent with the source terms S(U) in the sense that lim

U L ,U R →U ∆t,∆x→0 S(∆x, ∆t; U L , U R ) = S(U).
Hereafter and using very classical notations, (∆x j ) j∈Z and ∆t represent the constant time and variable space steps of the mesh under consideration for defining the approximate solutions. More precisely and in order to define the Godunov-type scheme associated with this approximate Riemann solver, we define the mesh interfaces x j+1/2 = x j-1/2 + ∆x j for j ∈ Z, and the intermediate times t n = n∆t for n ∈ N. Note that ∆x in [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] then plays the role of ∆x j+1/2 = 1 2 (∆x j + ∆x j+1 ). In the sequel, U n j denotes the approximate value of U at time t n and on the cell [x j-1/2 , x j+1/2 ). For n = 0 and j ∈ Z, we set U 0 j = 1 ∆x x j+1/2

x j-1/2 U 0 (x) dx where U 0 (x) is the initial condition. Then, the explicit in time Godunov-type scheme reads

                 U n+1 j = U n j - ∆t ∆x j (F n j+ 1 2 -F n j-1 2 ) + ∆t 2 ( ∆x j+1/2 ∆x j S n j+ 1 2 + ∆x j-1/2 ∆x j S n j-1 2
),

F n j+ 1 2 = F(U n j , U n j+1 ), S n j+ 1 2 = S(∆x j+1/2 , ∆t; U n j , U n j+1 ), (20)
with

F(U L , U R ) = 1 2 F(U L ) + F(U R ) - l k=1 |λ k |(U k+1 -U k ) .
As far as the consistency with the entropy inequality ( 16) is concerned, the simple approximate Riemann solver is said to be consistent with the integral form of ( 16) if and only if there exists a function σ such that under the CFL condition [START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF] we have

q(U R ) -q(U L ) -∆x σ(∆x, ∆t; U L , U R ) ≤ l k=1 λ k η(U k+1 ) -η(U k ) , (21)
with lim

U L ,U R →U ∆t,∆x→0 σ(∆x, ∆t; U L , U R ) = 0. ( 22
)
Then, the numerical scheme defined by [START_REF] Gosse | An asymptotic-preserving well-balanced schemes for the hyperbolic heat equations[END_REF] satisfies the following discrete entropy inequality

                 η(U n+1 j ) ≤ η(U n j ) - ∆t ∆x j (q n j+ 1 2 -q n j-1 2 ) + ∆t 2 ( ∆x j-1/2 ∆x j σ n j-1 2 + ∆x j+1/2 ∆x j σ n j+ 1 2
),

q n j+ 1 2 = q(U n j , U n j+1 ), σ n j+ 1 2 = σ(∆x j+1/2 , ∆t; U n j , U n j+1 ), (23)
with

q(U L , U R ) = 1 2 q(U L ) + q(U R ) - l k=1 |λ k | S(U k+1 ) -S(U k ) . (24)
The CFL condition associated with this explicit in time Godunov-type scheme naturally reads

max 1≤k≤l | λ k (U n j , U n j+1 ) | ∆t min(∆x j , ∆x j+1 ) ≤ 1 2 ,
for all j. Again, we refer to [START_REF] Gallice | Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source[END_REF][START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF][START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF] for more details.

5 Application to the Lagrangian system and explicit in time Godunovtype scheme

We suppose again that a is a parameter that complies with the subcharacteristic constraint [START_REF] Caflisch | Uniformly accurate schemes for hyperbolic systems with relaxation[END_REF]. We consider a step ∆m of the space variable expressed through the mass variable and a time step ∆t. The objective of this section is to define a consistent simple approximate Riemann solver for [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF]. We have in this case

U =     τ - → w ← - w E     , F(U) =     -u a - → w -a ← - w Πu     , S(U) =     0 a(g -αu) -a(g -αu) u(g -αu)     , F(U L , U R ) =     F τ (U L , U R ) F - → w (U L , U R ) F ← - w (U L , U R ) F E (U L , U R )     . m/t = -a m/t = 0 m/t = +a ∆m L ∆m R ∆t U * L U L U * R U R m t Figure 2:
Wave structure of the approximate Riemann solver for the Lagrangian system [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF] in the (m, t)-plane.

Note in particular that introducing the notation

F(U) = (F τ , F - → w , F ← -
w , F E ) T , the energy flux satisfies the relation

F E = - F - → w + F ← - w 2a F τ . ( 25
)
This relation will be used in the calculations below.

In order to mimic the self-similar solution to (12) when the source terms are omitted, we propose to consider a simple approximate Riemann solver made of 3 waves (see figure 2), namely a stationary wave and two waves propagating with velocities ±a :

W m t ; U L , U R =        U L , m t < -a, U * L , -a < m t < 0, U * R , 0 < m t < a, U R , m t > a.
Following [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF], we define S as follows :

S(∆m, ∆t; U L , U R ) =     0 a(g -αũ) -a(g -αũ) ũ(g -αũ)     ,
where ũ represents a consistent approximation of the velocity u, in the sense that lim

U L ,U R →U ∆t,∆x→0 ũ(∆m, ∆t; U L , U R ) = u
The definition of ũ will be specified later on.

We now turn to the definition the intermediate states U * L and U * R . Each state containing four components, eight relations are expected. As motivated in the previous section, we first impose the consistency relations [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] which gives here

         (u L -u R ) = -a(τ * L -τ L ) + a(τ R -τ * R ), a( - → w R -- → w L ) -∆m a(g -αũ) = -a( - → w * L -- → w L ) + a( - → w R -- → w * R ), -a( ← - w R -← - w L ) + ∆m a(g -αũ) = -a( ← - w * L -← - w L ) + a( ← - w R -← - w * R ), (Π R u R -Π L u L ) -∆m ũ(g -αũ) = -a(E * L -E L ) + a(E R -E * R ). (26)
Then, we make the natural choice of imposing the Rankine Hugoniot relations associated with the mass conservation across each wave of the approximate Riemann solver. We get

   u L -aτ L = u * L -aτ * L , u R + aτ R = u * R + aτ * R , u * L = u * R . (27)
Note however that (27) provides only two independent relations since the first equation of ( 26) is a linear combination of the three equations in [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF]. Then, two equations are still missing. In the sequel we shall note u * = u * L = u * R . In order to account for the source terms, we propose to impose a generalized jump condition across the stationnary wave and associated with the momentum equation in [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF]. This amounts to take into account the source term at the interface of the initial condition. More precisely, we impose

Π * R -Π * L = ∆m(g -αũ). (28)
So that only one equation is now missing. At last and in regards to the energy equation, we then propose to mimic the relation ( 25) at the discrete level by imposing the following relation on the numerical flux of the Godunov-type method:

F E (U L , U R ) = -F τ (U L , U R ) × F - → w (U L , U R ) + F ← - w (U L , U R ) 2a . ( 29
)
It remains to define ũ and following [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF] we set ũ = u * .

This choice aims at proposing the same approximation of the velocity u in the mass flux at the interface (that is u * ) as in the friction term (that is ũ). As mentioned in section 2, this is also true at the continuous level for the parabolic system (4). It turns out to be essential in order to obtain the asymptotic preserving property.

At this point, we can now define the intermediate states U * L and U * R . We have after easy calculations

                                   u * = 1 2a + α∆m a(u R + u L ) -(Π R -Π L ) + g∆m , τ * L = τ L + u * -u L a , τ * R = τ R + u R -u * a , Π * R = Π R + Π L 2 -a u R -u L 2 + (g -αu * )∆m 2 , Π * L = Π * R -(g -αu * )∆m, E * L = E L + 1 a p L u L -u * (p * - ∆m 2 g -αu * ) , E * R = E R -1 a p R u R -u * (p * + ∆m 2 g -αu * ) , (30)
where we have set

p * = Π R + Π L 2 -a u R -u L 2 = - → w L + ← - w R 2 .
Regarding the numerical flux of the Godunov-type scheme, we find

             F τ (U L , U R ) = -u * , F - → w (U L , U R ) = a - → w L + a(g -αu * )∆m 2 = a 2 u * + ap * , F ← - w (U L , U R ) = -a ← - w R + a(g -αu * )∆m 2 = a 2 u * -ap * , F E (U L , U R ) = p * u * .
We now use the above flux definition to derive an explicit in time Godunov-type scheme for [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF]. We consider the approximate variable U n j known for j ∈ Z and we set ∆m j = ρ n j ∆x, ∆m j+1/2 = ∆m j + ∆m j+1 2 .

The superscript Lag will denote the updated values after the approximation of [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF]. Following [START_REF] Gosse | An asymptotic-preserving well-balanced schemes for the hyperbolic heat equations[END_REF], we are led to the following scheme

                                     τ Lag j = τ n j + ∆t ∆m j (u * j+ 1 2 -u * j-1 2 ), - → w Lag j = - → w n j -a ∆t ∆m j ( - → w n j -- → w n j-1 ) + ∆t a ∆m j-1/2 ∆m j (g -αu * j-1 2 ), ← - w Lag j = ← - w n j + a ∆t ∆m j ( ← - w n j+1 -← - w n j ) -∆t a ∆m j+1/2 ∆m j (g -αu * j+ 1 2
),

E Lag j = E n j - ∆t ∆m j ((up) * j+ 1 2 -(up) * j-1 2 ) + g ∆t ∆m j ∆m j+1/2 u * j+ 1 2 + ∆m j-1/2 u * j-1 2 2 -α ∆t ∆m j ∆m j+1/2 (u * j+ 1 2 ) 2 + ∆m j-1/2 (u * j-1 2 ) 2 2 , (31)
where

u * j+ 1 2 = 1 2a + α∆m j+1/2 ( - → w n j -← - w n j+1 + g∆m j+1/2 ), p * j+ 1 2 = - → w n j + ← - w n j+1 2 . ( 32
)
Notice that the source terms clearly appear to receive an upwind treatment in [START_REF] Suliciu | Ont the thermodynamics of fluids with relaxation and phase transitions[END_REF]. This scheme is easily shown to be stable under the CFL condition given by max

j∈Z ∆t ∆m j a ≤ 1 2 , (33)
for all j. Let us note that as a satisfies the subcharacteristic condition [START_REF] Caflisch | Uniformly accurate schemes for hyperbolic systems with relaxation[END_REF] this CFL condition naturally involves the sound speed c.

6 Implicit in time Godunov-type scheme for the Lagrangian system

We follow here a standard approach for deriving an implicit in time Godunov-type scheme from (31) by simply replacing the terms evaluated at time t n with the terms noted with the superscript Lag . We get

       - → w Lag j = - → w n j -a ∆t ∆m j ( - → w Lag j -- → w Lag j-1 ) + ∆t a ∆m j-1/2 ∆m j (g -αu * j-1 2 ), ← - w Lag j = ← - w n j + a ∆t ∆m j ( ← - w Lag j+1 -← - w Lag j ) -∆t a ∆m j+1/2 ∆m j (g -αu * j+ 1 2 ), ( 34
) with u * j+ 1 2 = 1 2a + α∆m j+1/2 ( - → w Lag j -← - w Lag j+1 + g∆m j+1/2 ), (35)
and

                 τ Lag j = τ n j + ∆t ∆m j (u * j+ 1 2 -u * j-1 2
),

E Lag j = E n j - ∆t ∆m j ((up) * j+ 1 2 -(up) * j-1 2 ) + g ∆t ∆m j ∆m j+1/2 u * j+ 1 2 + ∆m j-1/2 u * j-1 2 2 -α ∆t ∆m j ∆m j+1/2 (u * j+ 1 2 ) 2 + ∆m j-1/2 (u * j-1 2 ) 2 2 , (36)
with

p * j+ 1 2 = - → w Lag j + ← - w Lag j+1 2 . ( 37
)
It is important to notice that (34) is independent from (36). More precisely, once (34) is solved, the update values (36) for τ and E follow explicitly. As far as -→ w and ←w are concerned, the update formulas (34) are coupled and require the resolution of a linear system. The corresponding matrix is shown to be pentadiagonal and strictly diagonally dominant. Therefore, it is invertible and (34) can be solved for any ∆t > 0. The proposed implicit in time numerical scheme for solving the Lagrangian system is then actually cheap thanks to the relaxation strategy.

7 Implicit-explicit in time Godunov-type scheme for the Eulerian system (1)

In order to complete the definition of our numerical scheme, it remains to define the second step of the operator splitting associated with the Lagrange-Projection decomposition, which corresponds to the Eulerian projection [START_REF] Bouchut | A reduced stability condition for nonlinear relaxation to conservative laws[END_REF]. Following [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF], we consider a very classic upwind and explicit in time numerical scheme given by

X n+1 j = X Lag j + ∆t ∆x (u * j-1 2 ) + X Lag j-1 + (u * j+ 1 2 ) --(u * j-1 2 ) + X Lag j -(u * j+ 1 2 ) -X Lag j+1 , (38)
where X ∈ {ρ, ρu, ρE} and

u + = u + |u| 2 , u -= u -|u| 2 ,
for all u. The update formula (38) is shown to be stable under the CFL condition given by ∆t ∆x (u * j-1

2

) + -(u * j+ 1 2 ) -< 1, (39)
see again [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] for instance. Note that this CFL condition involves only the flow speed and is not based on the acoustic waves.

For the sake of clarity, let us briefly recall the different steps of the overall method called LP-IMEX, suppose that at the instant n we know (ρ n j , (ρu) n j , (ρE) n j ) for j ∈ Z, we perform the following steps:

(i) compute (τ n j , - → w n j , ← - w n j , E n j ) at equilibrium by evaluating Π n j = p(ρ n j , e n j ), (ii) compute (τ Lag j , - → w Lag j , ← - w Lag j , E Lag j
) thanks to the implicit scheme defined by ( 34)-( 36),

(iii) evaluate (ρ Lag j , (ρu) Lag j , (ρE) Lag j ) thanks to (τ Lag j , - → w Lag j , ← - w Lag j , E Lag j ), (iv) compute (ρ n+1 j , (ρu) n+1 j , (ρE) n+1 j
) thanks to the explicit scheme defined by ( 38) and (35).

In the numerical experiments and for the sake of comparison, we will also consider the following explicitexplicit numerical scheme that will be referred to as the LP-EXEX scheme. Suppose that at the instant n we know (ρ n j , (ρu) n j , (ρE) n j ) for j ∈ Z, we perform the following steps: (i) compute (τ n j , -→ w n j , ←w n j , E n j ) at equilibrium by evaluating Π n j = p(ρ n j , e n j ), (ii) compute (τ Lag j , -→ w Lag j , ←w Lag j , E Lag j ) thanks to the explicit scheme defined by ( 31),

(iii) evaluate (ρ Lag j , (ρu) Lag j , (ρE) Lag j ) thanks to (τ Lag j , -→ w Lag j , ←w Lag j , E Lag j ), (iv) compute (ρ n+1 j , (ρu) n+1 j , (ρE) n+1 j ) thanks to the explicit scheme defined by ( 38) and [START_REF] Weyl | Shock waves in arbitrary fluids[END_REF].

Main properties

We give in this section the main properties of the proposed numerical scheme.

Theorem 1. Under the CFL condition (39), the implicit-explicit in time numerical scheme LP-IMEX is well defined and satisfies the following stability properties:

(i) it is a conservative scheme for the density ρ. It is also a conservative scheme for ρu and ρE when the source terms are omitted, (ii) the density ρ n j is positive for all j and n > 0 provided that ρ 0 j is positive for all j, (iii) it is asymptotic preserving. In addition, under the CFL condition (33), the explicit-explicit overall numerical scheme LP-EXEX (iv) satisfies an entropy inequality.

This result is worth a few comments. We first note that the overall numerical scheme is stable under the CFL condition (39) which only involves the flow speed and not the acoustic waves. It is then less restrictive than the classical CFL restrictions of the usual explicit Godunov-type numerical schemes. The CFL condition (39) involves however the solution computed at the end of the Lagrangian step and then is not determined explicitly. Finding an explicit formula for the time step restriction that guarantees (39) and then the stability of the scheme is an open question at the moment.

Properties (i) and (ii) are obtained from standard manipulations, see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] and [START_REF] Coquel | Entropy-satisfying relaxation method with large timesteps for Euler IBVPs[END_REF]. We first prove (iii) and then (iv).

Proof of (iii)

This section aims at proving that the proposed LP-IMEX numerical scheme is asymptotic preserving. More precisely and similarly to the continuous analysis proposed in section 2, we aim at proving that if we perform the scale change

∆t → ∆t , α → α , (40)
which corresponds to the study of the large time behaviour for a large friction coefficient, then we get a consistent approximation of (4) when goes to zero. We propose a two-part proof: we first prove the property in Lagrangian coordinates and then consider the Eulerian framework.

Asymptotic preserving property in Lagrangian coordinates

In Lagrangian coordinates, the limit system (4) reads

   ∂ t τ -∂ m u 1 = 0, ∂ m p = g -αu 1 , ∂ t e + ∂ m pu 1 = u 1 (g -αu 1 ). ( 41
)
Let us perform the scale change (40) in ( 34)-( 36). We get

                                     - → w Lag j = - → w n j -a ∆t ∆m j ( - → w Lag j -- → w Lag j-1 ) + ∆ta ∆m j-1/2 ∆m j (g - α u * j-1 2 ), ← - w Lag j = ← - w n j + a ∆t ∆m j ( ← - w Lag j+1 -← - w Lag j ) - ∆ta ∆m j+1/2 ∆m j (g - α u * j+ 1 2 ), τ Lag j = τ n j + ∆t ∆m j (u * j+ 1 2 -u * j-1 2
),

E Lag j = E n j - ∆t ∆m j (up) * j+ 1 2 -(up) * j-1 2 + g∆t ∆m j ∆m j+1/2 u * j+ 1 2 + ∆m j-1/2 u * j-1 2 2 - α∆t 2 ∆m j ∆m j+1/2 (u * j+ 1 2 ) 2 + ∆m j-1/2 (u * j-1 2 ) 2 2 , ( 42
) with u * j+ 1 2 = 2a + α∆m j+1/2 ( - → w Lag j -← - w Lag j+1 + g∆m j+1/2 ) and p * j+ 1 2 = - → w Lag j + ← - w Lag j+1 2 .
Multiplying the first equation of ( 42) by , we obtain -→ w Lag j -←w Lag j = 0 + O( ). This implies by definition of -→ w and ← -

w that for j ∈ Z u Lag j = 0 + O( ), (43)
which also gives p *

j+ 1 2 = P * j+1/2 + O( ), where
P * j+1/2 = Π Lag j + Π Lag j+1 2
We also have thanks to (43) that

u * j+ 1 2 = 2a + α∆m j+1/2 ( - → w Lag j -← - w Lag j+1 + g∆m j+1/2 ) = v * j+ 1 2 + O(ε 2 ), (44)
where

v * j+ 1 2 = 1 α g - Π Lag j+1 -Π Lag j ∆m j+1/2 ,
and thus provide a definition of v * j+ 1 2 that is consistent with the second equation in (41). By reinjecting the expression of u * j+ 1 2 in the last two equations of (42), we obtain finally

                                       τ Lag j =τ n j + ∆t ∆m j (v * j+ 1 2 -v * j-1 2 ) + O( ), v * j+ 1 2 = 1 α (g - Π Lag j+1 -Π Lag j ∆m j ), P * j+ 1 2 = Π Lag j + Π Lag j+1 2 , e Lag j =e n j - ∆t ∆m j ((vP ) * j+ 1 2 -(vP ) * j-1 2 ) + g ∆t ∆m j ∆m j+1/2 v * j+ 1 2 + ∆m j-1/2 v * j-1 2 2 -α ∆t ∆m j ∆m j+1/2 (v * j+ 1 2 ) 2 + ∆m j-1/2 (v * j-1 2 ) 2 2 + O( ), (45)
which is consistent with (41) when tends to 0.

Asymptotic Preserving property in Eulerian coordinates

It remains to prove that after the Eulerian projection, the overall scheme is consistent with (4). In (38) we perform the scale change (40). We have for X ∈ {ρ, ρu, ρE}:

X n+1 j = X Lag j + ∆t ∆x (v * j-1 2 ) + X Lag j-1 + (v * j+ 1 2 ) --(v * j-1 2 ) + X Lag j -(v * j+ 1 2 ) -X Lag j+1 + O( ). (46)
Since we have u Lag

j = 0 + O( ), then (ρu) Lag j
= 0 + O( ) for all j ∈ Z and considering X = ρu in the previous equality gives

u n+1 j = 0 + O( ).
Let us remark that the first equation in (45) reads

ρ n j = ρ Lag j 1 + ∆t ∆x (v * j+ 1 2 -v * j-1 2 ) + O( ).
Then, for X ∈ {ρ, ρE} the Eulerian projection (46) may be recast into

(X) n+1 j = (X) Lag j 1 + ∆t ∆x (v * j+ 1 2 -v * j-1 2 ) -L a ((X) Lag , v * ) + O( ) = ρ n j X Lag j ρ Lag j -L a ((X) Lag , v * ) + O( ), (47)
where the advection operator L a ((X) Lag , v * ) is consistent with ∂ x (Xv) and is defined by

L a ((X) Lag , v * ) = ∆t ∆x (X) Lag j (v * j+ 1 2 ) + + (X) Lag j+1 (v * j+ 1 2 ) --(X) Lag j (v * j-1 2 ) -+ (X) Lag j-1 (v * j-1 2 ) + .
Injecting the last equation of ( 45) into (47) yields in the limit → 0

                                   ρ n+1 j =ρ n j -L a (ρ Lag , v * ) + O( ), v * j+ 1 2 = 1 α g -2 Π Lag j+1 -Π Lag j (ρ n j + ρ n j+1 )∆x , P * j+ 1 2 = Π Lag j + Π Lag j+1 2 , (ρe) n+1 j
=(ρe) n j - ∆t ∆x (vP ) * j+ 1 2 -(vP ) * j-1 2 +ρ n j ∆t g v * j+ 1 2 + v * j-1 2 2 -α (v * j+ 1 2 ) 2 + (v * j-1 2 ) 2 2 -L a ((ρe) Lag , v * ) + O( ),
which is clearly consistent with (4). Hence, the overall scheme is asymptotic preserving is the proof is completed.

Proof of (iv)

We now prove that the explicit scheme LP-EXEX satisfies a discrete entropy inequality. As above, we first prove such an inequality in Lagrangian coordinates considering the scheme (31) that may be recast into

                             τ Lag j =τ n j + ∆t ∆m j (u * j+ 1 2 -u * j-1 2 ), u Lag j =u n j - ∆t ∆m j (p * j+ 1 2 -p * j-1 2 ) + ∆t ∆m j ∆m j-1/2 g -αu j-1 2 + ∆t ∆m j ∆m j+1/2 g -αu j+ 1 2 , E Lag j =E n j - ∆t ∆m j (up) * j+ 1 2 -(up) * j-1 2 + g ∆t ∆m j ∆m j+1/2 u * j+ 1 2 + ∆m j-1/2 u * j-1 2 2 -α ∆t ∆m j ∆m j+1/2 (u * j+ 1 2 ) 2 + ∆m j-1/2 (u * j-1 2 ) 2 2 , (48) with u
* j+ 1 2 = 1 2a + α∆m j+1/2 a(u n j+1 + u n j ) -p(ρ n j+1 , e n j+1 ) -p(ρ n j , e n j ) + g∆m j+1/2 and p * j+ 1 2 = p(ρ n j+1 , e n j+1 ) + p(ρ n j , e n j ) 2 - a(u n j+1 -u n j ) 2 .
Then we take the Eulerian projection (38) into account to obtain a discrete entropy inequality for the overall scheme in Eulerian coordinates. Let s = s(ρ, e) be the strictly convex mathematical entropy. The entropy inequality associated with (1) writes

∂ t (ρs) + ∂ x (ρus) ≤ 0. (49)
In Lagrangian coordinates, ( 7) is associated with the following entropy inequality

∂ t s ≤ 0. ( 50
)
In the sequel and with a little abuse in the notations, we will consider the pressure as a function of the density and the entropy p = p(ρ, s).

Entropy inequality in Lagrangian coordinates

The scheme (48) is a Godunov-type scheme. We prove here that the associated approximate Riemann solver is consistent with the entropy inequality (50). Let us check that

0 ≤ -a(s * L -s L ) + a(s R -s * R ), (51)
where s * L = s(ρ * L , e * L ) and s * R = s(ρ * R , e * R ) so that ( 21) is verified with η = s, q = 0 and σ(∆m, ∆t; U L , U R ) = 0. Let us first prove the following result.

Proposition 1. If a > 0 is such that    ρ * L > 0, ρ * R > 0, ρ 2 ∂ ρ p(ρ, s L ) ≤ a 2 , ∀ρ ∈ I(ρ L , ρ * L ), ρ 2 ∂ ρ p(ρ, s R ) ≤ a 2 , ∀ρ ∈ I(ρ R , ρ * R ), ( 52
)
then we have e * L ≥ e(ρ * L , s L ), e * R ≥ e(ρ * R , s R ). (53)
Proof. Recall that we have

     E * L = E L + 1 a p L u L -u * [p * - ∆m 2 (g -αu * )] , E * R = E R + 1 a u * [p * + ∆m 2 (g -αu * )] -p R u R , and
     Π * R = p * + ∆m 2 (g -αu * ), Π * L = p * - ∆m 2 (g -αu * ), which gives      E * L = E L - 1 a (Π * L u * -Π L u L ), E * R = E R + 1 a (Π * R u * -Π R u R ), so that e * R = E * R - 1 2 u * 2 = e R + 1 2 (u 2 R -u * 2 ) + 1 a (Π * R u * -Π R u R ).
Using the definitions of u * and Π * R given in [START_REF] Natalini | Convergence to equilibrium for the relaxation approximations of conservation laws[END_REF], straightforward calculations then lead to

e * R = e R + 1 2a 2 (Π * 2 R -Π 2 R ) and Π * R = Π R + a(u * -u R
). The latter equality, together with the last two equalities in [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF], gives

Π * R a 2 = 1 ρ R + p R a 2 - 1 ρ * R .
It is then not difficult to check that

       e * R -e(ρ * R , s R ) = 1 2a 2 (p(ρ * R , s R ) -Π * R ) 2 + φ(ρ R ), φ(ρ) = e(ρ, s R ) - p(ρ, s R ) 2 2a 2 -e(ρ * R , s R ) + p(ρ * R , s R ) 2 2a 2 + p(ρ * R , s R ) 1 ρ + p(ρ, s R ) a 2 - 1 ρ * R - p(ρ * R , s R ) a 2 ,
with, using the well-known relation [START_REF] Berthon | An HLLC scheme to solve the M1 Model of radiative transfer in two space dimensions[END_REF],

   φ (ρ) = (p(ρ, s R ) -p(ρ * R , s R ))( 1 ρ 2 - 1 a 2 ∂ ρ p(ρ, s R )), φ(ρ * R ) = 0.
Therefore φ(ρ) ≥ 0 for all ρ ∈ I(ρ R , ρ * R ) since ∂ ρ p ≥ 0 under the assumptions of Weyl and by [START_REF] Caflisch | Uniformly accurate schemes for hyperbolic systems with relaxation[END_REF]. We have thus proved in particular that e * R ≥ e(ρ * R , s R ). The proof of the second inequality follows the same idea and the proof of the proposition is then completed.

In order to get (51), we then note that the temperature T being positive by assumption, we have ∂ e s(ρ, e) < 0 where s = s(ρ, e(ρ, s)). So that

s R -s * R = s(ρ * R , e(ρ * R , s R )) -s(ρ * R , e * R ) = ∂ e s(ρ * R , e * L )(e(ρ * R , s R ) -e * R ) ≥ 0, that is s R -s * R ≥ 0.
We get in the same way s Ls * L ≥ 0 which gives (51) since a > 0. Then, under the CFL condition (33), the scheme in Lagrangian coordinates satisfies the discrete entropy inequality [START_REF] Jin | Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms[END_REF] which reads

s Lag j ≤ s n j - ∆t ∆m j (q n j+ 1 2 -q n j-1 2 ), (54)
where q n j+ 1 2 naturally follows from [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF].

Entropy inequality in Eulerian coordinates

Let us now prove an entropy inequality for the whole Lagrange-Projection scheme. We first recall that the function ρs = ρs(ρ, ρu, ρE) is strictly convex by assumption, see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]. The Eulerian projection (38) being a convex combination for each X = ρ, ρu, ρE, we thus have

(ρs) n+1 j ≤ (ρs) Lag j + ∆t ∆x (u * j-1 2 ) + (ρs) Lag j-1 + (u * j+ 1 2 ) --(u * j-1 2 ) + (ρs) Lag j -(u * j+ 1 2
) -(ρs) Lag j+1 .

Then using the first equation of ( 48) in ( 54), together with the definition of ∆m j , we have

(ρs) Lag j ≤ (ρs) n j - ∆t ∆x (q n j+ 1 2 -q n j-1 2 ) - ∆t ∆x (ρs) Lag j (u * j+ 1 2 -u * j- 1 2
).

Then using these two inequalities together with the simple relation u = u + + u -, we get

(ρs) n+1 j ≤ (ρs) n j - ∆t ∆x (g n j+ 1 2 -g n j-1 2 ), (55)
where

g n j+ 1 2 = (u * j+ 1 2 ) + (ρs) Lag j + (u * j+ 1 2 ) -(ρs) Lag j+1 + q n j+ 1 2 .
Which is nothing but the expected discrete entropy inequality for the whole Lagrange-Projection scheme (note indeed that g n j+ 1

2 is clearly consistent with the entropy flux ρsu).

Numerical results

For the sake of stressing the importance of the source term discretization, we propose to also consider in the following a more simple approximation strategy that will be referred to as the LP-EXEX SP scheme. In this scheme, the source terms will be treated by means of a separate operator splitting. Suppose that for some instant n we know (ρ n j , (ρu) n j , (ρE) n j ) for j ∈ Z, the LP-EXEX SP numerical scheme reads: (i) compute (τ n j , -→ w n j , ←w n j , E n j ) at equilibrium by evaluating Π n j = p(ρ n j , e n j ), (ii) compute (τ Lag ) thanks to the explicit scheme defined by (38) and ( 32), (v) account for the gravity and friction source terms by integrating the system of ordinary differential equations

d dt   ρ ρu ρE   =   0 ρ(g -αu) ρu(g -αu)   .
The update in term of the variable ρ, u, and e then reads

     ρ n+1 j = ρ n+1, j , u n+1 j = u n+1, j e -α∆t + g α (1 -e -α∆t ), e n+1 j = e n+1, j . ( 56
) (vi) evaluate (ρ n+1 j , (ρu) n+1 j , (ρE) n+1 j
). We propose to test both LP-EXEX SP and LP-IMEX scheme against a test case that has been proposed in [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF]. In the sequel, we shall consider that the fluid is equipped with a perfect gas equation of state p = (γ-1)ρe and we set the gravity acceleration, the friction coefficient and the specific heat ratio to the following values

g = 9.81 m • s -2 , α = 10 6 s -1 , γ = 1.4.
The initial condition is defined by At the boundaries, we impose periodic boundary conditions thanks to a fictious cell at each end of the domain.

In the sequel, both LP-EXEX and LP-EXEX SP computations will be performed with a time step defined by

∆t = min ρ n j ∆x 2a , (57)
in order to agree with the classic acoustic CFL (33). The choice of the time step for the LP-IMEX scheme will be specified case by case. For each test, we will compute a reference solution thanks to the LP-EXEX scheme over a 10 000-cell grid. If we refer to this solution thanks to the superscript ref and if Y denote a fluid variable, for the sake of comparison we shall consider in the sequel the L 1 relative error with respect to the reference solution at the instant t defined by

err(Y, t) = Y (•, t) -Y ref (•, t) L 1 ([0,1]) Y ref (•, t) L 1 ([0,1]) .
9.1 Test case 1: sensitivity with respect to the space step for large friction

We run our numerical tests with the LP-EXEX SP scheme using a spatial discretization over 100 cells, 1000 cells and 10 000 cells. This leads to time step values ∆t that are respectively of magnitude 10 α , 1 α and 1 10α . In figure 3, we display the result obtained at t = 0.01 s and we can see that there is a large amount of numerical diffusion due to the discretization of the source term with large values of α. It is necessary to choose small values of ∆t relatively to 1 α in order to preserve the accuracy of the solution. We now consider the same test performed with the LP-IMEX scheme. We now choose ∆t in agreement with the CFL condition (39) by setting

∆t = min ∆x 2 max u n j , 1 α , (58)
so that we always have ∆t ≤ 1 α . The results obtained with LP-IMEX scheme at instant t = 0.01 s are presented in figure 4 for discretization grids of 100 cells and 1000 cells. It is clear that the approximate solution is much more accurate than the one computed with the LP-EXEX SP scheme, even for coarse mesh. Let us note that for the relatively large space steps that have been used here, choice (58) imposes ∆t = 1 α . Table 1 displays the relative error obtained with both the LP-EXEX SP and the LP-IMEX scheme. It shows that the asymptotic preserving property significantly lessen the numerical diffusion and then improve accuracy by several orders of magnitude.

Test case 2: sensitivity with respect to the time step

We are now interested in testing the LP-IMEX scheme in situations where ∆t is much bigger than 1 α . In order to procede, we relax the previous time step choice (58) by suppressing the control of with respect to 1 α . We simply choose to define ∆t in agreement with the CFL condition (39) based on the material velocity by setting ∆t = min ∆x 2 max u n j .

(59) The test is performed with a 1000-cell mesh. The graph of the approximate solution at t = 0.01 s are displayed in figure 5. For this grid choice, we obtain that the magnitude of ∆t is 1000 α . Let us underline that these time steps are 1000 times larger than the time steps used in section 9.1. It appears that even for such relatively large time step and space step the approximate solution remains very accurate. This results is the consequence of the good behaviour of the numerical scheme for large friction coefficients. We compare in table 2 the relative errors of the fluid variables for this choice of time step values. We can verify here that the LP-IMEX scheme enables the use of time steps that are much larger than the acoustic based time step and also much larger than 1 α while preserving accurate simulation results.

Conclusion

We designed a large time step and asymptotic preserving scheme for the gas dynamics equations with gravity and linear friction. The stability is proved under a time step CFL restriction based on the velocity u only and not on the sound speed c. The scheme gives very good results and combines accuracy thanks to the asymptotic preserving property and efficiency thanks to the large time step stability condition. Future developments include an extension to several space dimensions and high-order accuracy, to low-Mach number flows and to more complicated systems of partial differential equations used in the modeling of two phase flows.

1 = -ρc < λ Lag 2 = 0 < λ Lag 3 =

 123 ρc, where c still denotes the sound speed. Here again, the extreme characteristic fields associated with λ Lag 1 and λ Lag 3 are genuinely non linear while the intermediate characteristic field associated with λ Lag 2

,

  the use the explicit scheme defined by[START_REF] Suliciu | Ont the thermodynamics of fluids with relaxation and phase transitions[END_REF] with α = 0 and g = 0, so that the gravity and friction terms are not taken into account yet (iii) evaluate (ρ Lag j

  (ρ, u, p) = (1.0, 0, 10000.0), if x ∈[0, 0.35] ∩ [0.65, 1], (ρ, u, p) = (2.0, 0, 26390.2), if x ∈[0.35, 0.65].

Figure 3 :

 3 Figure 3: Profile at time t = 0.01 s of the velocity (top left), the pressure (top right), the density (bottom)obtained for a 100-cell, 1000-cell and 10 000-cell grid with the LP-EXEX SP scheme and the reference solution (LP-EXEX scheme with 10 000-cell mesh).

Figure 4 :

 4 Figure 4: Profile at time t = 0.01 s of the velocity (top left), the pressure (top right), the density (bottom) obtained for a 100-cell and 1000-cell grid with the LP-IMEX scheme and the reference solution (LP-EXEX scheme with 10 000-cell mesh).

α 2 . 3 -

 23 607495 × 10 -3 1.099137 × 10 -1 3.288768 × 10 -

Figure 5 :

 5 Figure 5: Profile at time t = 0.01 s of the velocity (top left), the pressure (top right), the density (bottom) obtained for a 1000-cell grid with the LP-IMEX scheme and the reference solution (LP-EXEX scheme with 10 000-cell mesh).

  Figure 1: Wave structure of the approximate Riemann solver in the (x, t)-plane.

	W( x

t ; U L , U R ) made of l + 1 intermediate states U k separated by discontinuities propagating with velocities λ k (see figure

1

), namely

Table 1 :

 1 Comparison of the relative errors between the approximated solutions obtained with both LP-EXEX SP and LP-IMEX schemes. The space domain is discretized with a 1 000-cell space discretization and ∆t = 1 EXEX SP 1.686931 × 10 -2 6.858335 × 10 -1 2.539820 × 10 -2 LP-IMEX 3.959560 × 10 -4 1.195630 × 10 -2 5.635518 × 10 -4

	α

Table 2 :

 2 Comparison of the relative L 1 -errors obtained with the LP-IMEX scheme for a 1 000-cell space discretization and two different ∆t values. × 10 -4 1.195630 × 10 -2 5.635518 × 10 -4 LP-IMEX 1000

	numerical scheme	∆t	err(ρ, t = 0.01)	err(u, t = 0.01)	err(P, t = 0.01)
	LP-IMEX	1 α	3.959560

Acknowledgement. This present work has been achieved within the framework of the research project LAT-SAP. This project was funded by the CEA/DEN/DANS/DM2S for the 2011 edition of the CEMRACS.