Cyril Agut

Julien Diaz

Mathematical Modelling and Numerical Analysis Will be set by the publisher Modélisation Mathématique et Analyse Numérique STABILITY ANALYSIS OF THE INTERIOR PENALTY DISCONTINUOUS GALERKIN METHOD FOR THE WAVE EQUATION

Keywords: 1991 Mathematics Subject Classication. 35L05, 65M12, 65M60. Discontinuous Galerkin, penalization coecient, CFL condition, wave equation

We consider here the Interior Penalty Discontinuous Galerkin (IPDG) discretization of the wave equation. We show how to derive the optimal penalization parameter involved in this method in the case of regular meshes. Moreover, we provide necessary stability conditions of the global scheme when IPDG is coupled with the classical Leap-Frog scheme for the time discretization. Numerical experiments illustrate the fact that these conditions are also sucient.

Introduction

Nowadays, in many domains as for instance medical imaging, seismic imaging, radar, earthquakes simulations, one requires the accurate numerical solution to acoustic, elastodynamic or electromagnetic wave equations. In many cases, the equations have to be solved in very large domains with strong heterogeneities. Therefore, one has to use sophisticated discretization methods to compute the most accurate solution for the smallest computational cost.

Finite Dierences Methods, which are widely used because of their small computational cost and the simplicity of their implementation, are not adapted to deal with strong heterogeneities. Indeed, they rely on structured grids, which can not accurately approximate the shape of the various layers of the domain. Finite Elements Methods (FEM) are more adapted to this kind of problems since they allow for the use of unstructured grids. Among all existing FEM, the Spectral Element Method [START_REF] Cohen | Higher-Order Numerical Methods for Transient Wave Equations[END_REF][START_REF] Cohen | Higher-order nite elements with mass-lumping for the 1d wave equation[END_REF][START_REF] Komatitsch | Introduction to the spectral element method for three-dimensional seismic wave propagation[END_REF][START_REF] Seriani | Spectral element method for acoustic wave simulation in heterogeneous media[END_REF], is probably one of the most ecient to solve the wave equation since the resulting mass matrix is diagonal. Hence, it can be easily coupled to explicit time-schemes, which requires the inversion of the mass matrix at each time step. However, SEM requires the use of quadrilateral (in 2D) or hexahedral (in 3D) meshes which can be dicult to generate for realistic applications. Let us however mention that SEM can be extended to handle triangular meshes [START_REF] Cohen | Higher-order triangular nite elements with mass-lumping for the wave equation[END_REF], but it requires additional degrees of freedom and the implementation is more complex. Moreover, as far as we know, the extension to tetrahedral meshes has not been proposed yet.

Discontinuous Galerkin Methods are more and more popular for solving the wave equation since they lead to block-diagonal mass matrices without the help of quadrature formula. Moreover, they can be used with any type of meshes and even allow for the variation of the physical parameters inside the cells of the mesh (provided that the variation is polynomial). DGM are also naturally adapted to parallel computing since all volume integrals are computed locally and the communications between the cells are ensured by integrals over the faces of the elements. In [START_REF] Arnold | Unied analysis of disconitnuous galerkin methods for elliptic problems[END_REF], the authors provide a detailed review of the various Discontinuous Galerkin approximations of the Laplacian operator. They show that the so-called Interior Penalty Discontinuous Galerkin Method (IPDGM), also known as Symmetric Interior Penalty (SIP) [START_REF] Arnold | An interior penalty nite element method with discontinuous elements[END_REF][START_REF] Bassi | A high-order accurate discontinuous nite element method for the numerical solution of the compressible navier-stokes equations[END_REF], is one of the most suitable since it is stable and adjoint consistent, which guarantees the optimal order of convergence of the scheme. This explains why this method has been succesfully used to solve Helmholtz equation [START_REF] Ainsworth | Dispersive and dissipative properties of discontinuous galerkin nite element methods for the second-order wave equation[END_REF][START_REF] Alvarez | A discontinuous nite element formulation for helmholtz equation[END_REF] and the wave equation [START_REF] Ainsworth | Dispersive and dissipative properties of discontinuous galerkin nite element methods for the second-order wave equation[END_REF][START_REF] Grote | Discontinuous galerkin nite element method for the wave equation[END_REF][START_REF] Grote | Convergence analysis of a fully discrete dicontinuous galerkin method for the wave equation[END_REF]. Comparisons of the performances of IPDGM and SEM can be found in [START_REF] De Basabe | Stability of the high-order nite elements for acoustic or elastic wave propagation with high-order time stepping[END_REF] and [START_REF] Baldassari | Modélisation et simulation numérique pour la migration terrestre par équation d'ondes[END_REF]. It is worth noting that the rst paper concluded that the performances of SEM are better than IDPGM when IPDGM is applied to structured grids composed of squares, while the second paper show that the performances of IPDGM are better when it is applied to triangular meshes.

Nevertheless, in spite of all its interesting properties, IPDGM still suers from two diculties. The rst one is the determination of the penalization parameter, which penalizes the discontinuities of the solution through the faces. The accurate determination of the optimal parameter is crucial, since a too small value leads to instabilities while a too large value could (strongly) hamper the CFL (Courant-Friedrichs-Lewy) condition, which gives the maximal time step that can be used to ensure the stability of the scheme. In [START_REF] Ainsworth | Dispersive and dissipative properties of discontinuous galerkin nite element methods for the second-order wave equation[END_REF], the authors conjectured a minimal value of the penalization parameters, depending on p, the polynomial degree of the basis functions and on the size of the elements. They proved their conjecture up to p = 3. The extension of this result to p > 3 and to unstructured meshes, is still to be done. The second diculty is the determination of the CFL condition. It is well-known that this condition decreases when the penalization parameter increases, but no analytical formula has been proposed yet. The aim of this paper is a) to prove the conjecture of Ainsworth, Monk and Muniz up to p = 5; b) to propose a solution methodology to prove it for a given p; and c) to provide an analytic formula linking the CFL condition to the penalization parameter. We restrict ourselves to the cases of structured meshes composed of segments (in 1D), squares (in 2D) or cubes (in 3D). In section 1, we recall the IPDG discretization of the wave equation. In section 2, we propose two theorems, the rst one provides explicit necessary stability conditions on the penalization parameter and the time step while the second one provides a more restrictive but implicit necessary stability condition. The proof of these theorems in the one dimensional case is given in section 3. Section 4 is devoted to the proof in the three dimensional case and contains a discussion on the adaptation of the theorems to structured meshes composed of rectangles or parallelepipeds. We do not present the proof in the two dimensional case, but it can be adapted without any diculties from the three dimensional case. Finally, we present numerical results in section 5 that illustrate the fact that the stability condition is actually necessary and sucient for numerical applications.

Interior Penalty Discontinuous Galerkin Discretization of the wave equation

In this section, we recall the so called Interior Penalty Discontinuous Galerkin method applied to the acoustic wave equation in homogeneous bounded media Ω ⊂ R d , d = 1, 2, 3. For the sake of simplicity, we impose homogeneous Dirichlet boundary conditions on the boundary Γ := ∂Ω but this study can be extended to Neumann boundary conditions without major diculties.

                     Find u : Ω × [0, T ] → R such that : 1 µ ∂ 2 u ∂t 2 -div 1 ρ ∇u = f in Ω × ]0, T ] , u (x, 0) = u 0 , ∂u ∂t (x, 0) = u 1 in Ω, u = 0 on ∂Ω. (1)
where u stands for the displacement, µ is the compressibility modulus, ρ is the density and f is the source term. We introduce a triangulation T h of Ω and the following space of approximation with piecewise discontinuous polynomial functions :

V h := v ∈ L 2 (Ω) : v |K ∈ P p (K) , ∀K ∈ T h .
The set of the mesh faces is denoted F h which is partitionned into two subsets F i h and F b h corresponding respectively to the interior faces and those located on the boundary. For F ∈ F i h , we note arbitrarily K + and K -the two elements sharing F and we dene ν as the unit outward normal vector pointing from K + to K -. Moreover, v ± represents the restriction of a function v to the element K ± and we dene the jump and the average of a piecewise smooth function

v ∈ V h over F ∈ F i h such that [[v]] = v + -v -, {{v}} = v + + v - 2 .
(2)

For F ∈ F b h , we dene [[v]] = v and {{v}} = v.
The IPDG discretization of (1) reads as

       Find u h ∈ V h such that, ∀v h ∈ V h : K∈T h K 1 µ ∂ 2 t u h v h dx = -a h (u h , v h ) + K∈T h K f v h dx. (3)
where a h is a bilinear form dened by

a h (u h , v h ) = B T h (u h , v h ) -I (u h , v h ) -I (v h , u h ) + B S (u h , v h ) , with B T h (u h , v) = K∈T h K 1 ρ ∇u h ∇v, I (u h , v) = F ∈F i h F [[v]] 1 ρ ∇u h • ν , B S (u h , v) = F ∈F i h F γ [[u h ]] [[v]] .
The bilinear form B S is devoted to enforce the coercivity of a h and the penalization function γ is dened on each interior face F by

γ = α ξ F
where α is a positive parameter. There are many denitions of the function ξ F in the litterature. The most commonly used are:

• ξ F = h (F )
where h (F ) denotes the diameter of F . See for instance [START_REF] Ainsworth | Dispersive and dissipative properties of discontinuous galerkin nite element methods for the second-order wave equation[END_REF][START_REF] Arnold | Unied analysis of disconitnuous galerkin methods for elliptic problems[END_REF][START_REF] Grote | Convergence analysis of a fully discrete dicontinuous galerkin method for the wave equation[END_REF]. It is worth noting that this denition does not make sense in 1D.

• ξ F = min(h (K + ) , h (K -))
where h (K ± ) is the diameter of K ± . See for instance [START_REF] Grote | Discontinuous galerkin nite element method for the wave equation[END_REF].

• ξ F = min(ρ (K + ) , ρ (K -)) where ρ (K ± ) is the diameter of the inscribed circle (or sphere) of K ± . See for instance [START_REF] Shahbazi | An explicit expression for the penalty parameter of the interior penalty method[END_REF]. Whatever the denition of ξ F , there exists α p 0 > 0 such that the coercivity of a h is ensured for all α ≥ α p 0 . Obviously, the optimal parameter α 0 depends on the choice of the basis functions of V h , but also on ξ F . It has been shown by Shabazi in [START_REF] Shahbazi | An explicit expression for the penalty parameter of the interior penalty method[END_REF], using inverse inequalities proposed in [START_REF] Warburton | On the constants in hp-nite element trace inverse inequalities[END_REF], that the third denition was the most appropriate for triangular meshes. In [START_REF] Ainsworth | Dispersive and dissipative properties of discontinuous galerkin nite element methods for the second-order wave equation[END_REF], the authors showed that, on square and cubic meshes, using the rst denition of ξ F , the optimal parameter is α p 0 = p (p + 1) 2

for p = 0, . . . , 3 and they conjectured this relation for p ≥ 4. It is worth noting that, for such meshes, the rst and the third denition of ξ F are equivalent. At this point, we choose not to explicit the expression of ξ F . This will be done in the next section.

We refer to [START_REF] Ainsworth | Dispersive and dissipative properties of discontinuous galerkin nite element methods for the second-order wave equation[END_REF][START_REF] Arnold | Unied analysis of disconitnuous galerkin methods for elliptic problems[END_REF][START_REF] Grote | Discontinuous galerkin nite element method for the wave equation[END_REF] for more details on the properties of the bilinear form a h .

Considering {ϕ i } i=1,...,m the classical discontinuous Lagrange basis functions of degree p of V h , where m denotes the number of degrees of freedom of the problem, we obtain the following linear system:

∂ 2 t U = M -1 KU + M -1 F (4)
where

(M ) i,j = K∈T h K ϕ i ϕ j , (K) i,j = a h (ϕ i , ϕ j ) , (F ) i = K∈T h K f ϕ i .
Now, we have to discretize in time. Using the well known Leap-Frog scheme, we obtain the following fully discretized scheme:

U n+1 -2U n + U n-1 ∆t 2 = -M -1 KU n + M -1 F n . (5)
Since ( 5) is an explicit scheme, its L 2 stability is constrained by a CFL condition. It is well known [START_REF] Grote | Convergence analysis of a fully discrete dicontinuous galerkin method for the wave equation[END_REF] that this CFL condition decreases when α increases and behaves as for large α. However, no explicit formula of the CFL condition has been proposed yet. This is the objective of the next section.

Stability analysis

In this section, we rst propose necessary conditions over γ and ∆t ensuring the L 2 -stability of scheme [START_REF] Baldassari | Modélisation et simulation numérique pour la migration terrestre par équation d'ondes[END_REF]. This theorem provides an explicit dependance of ∆t with respect to γ and h. Next we propose a more restrictive necessary stability condition. In this second theorem, the dependance of ∆t with respect to γ is no longer explicit. However the condition can be numerically computed using the roots of a polynomial of degree 2p. We assume here that the domain Ω is unbounded Ω = R d and uniformly meshed by segments (if d = 1), squares (if d = 2) or cubes (if d = 3). The length of the edges of the elements is denoted by h.

The necessary stability conditions are given by the following theorem.

Theorem 3.1. The scheme ( 5) is L 2 -stable only if, for p ≤ 5,

γ ≥ p(p + 1) 2h = α 0 p h ; (6)
and, denoting γ = α /h,

√ d c∆t h ≤    C 1,p if α ≤ α 1,p C 2,p (α) if α > α 1,p . (7)
where α 1,p , C 1,p and C 2,p (α) are dened with respect to the polynomial degree p such that:

p α 1,p C 1,p C 2,p (α) 1 2 √ 3 3 0.577 1 3 (α -1) 2 27 5 = 5.4 1 √ 15 0.258 2
-15 + 6α + (405 -240α + 36α 2 )
1 /2 where, for the case p = 4, we have

             g 4,1 (α) = 518 -98α + 5α 2 1 2 , g 4,2 (α) = cos 1 3 arccos 1 10 g 4,3 (α) √ 5 g 3 4,1 (α) , g 4,3 (α) = -47705 + 14574α -1470α 2 + 50α 3
and for the polynomials of degree 5,

                     g 5 = 1 3 arccos 10447 126350 √ 133 g 5,1 (α) = 1555 -200α + 7α 2 1 2 , g 5,2 (α) = cos 1 3 arccos 1 14 g 5,3 (α) √ 7 g 3 5,1 (α) , g 5,3 (α) = -299825 + 61440α -4200α 2 + 98α 3 .

Remark 3.2. a

• As it is was noted in [START_REF] Grote | Convergence analysis of a fully discrete dicontinuous galerkin method for the wave equation[END_REF], the stability condition on ∆t behaves as C / √ α for α large enough. More

precisely, C = 2 (p + 1) (p + 2)
.

• This stability condition is constant for

p (p + 1) 2 ≤ α ≤ α 1
p . This shows that it is not necessary to choose α too close from α 0 p to improve the CFL condition.

• In [START_REF] Ainsworth | Dispersive and dissipative properties of discontinuous galerkin nite element methods for the second-order wave equation[END_REF], they authors proved (6) for p = 0, . . . , 3 and conjectured this relation for any p. Theorem (3.3) extends its validity until p = 5. • The condition [START_REF] Bassi | A high-order accurate discontinuous nite element method for the numerical solution of the compressible navier-stokes equations[END_REF] does not depend on the dimension d. This would not have been the case if we had expressed γ as a function of the circumcircle (or circumsphere) diameter which is √ dh. Since h is the diameter of the inscribed circle or sphere, we conjecture that the third denition of ξ F is the most appropriate. We will strengthen this conjecture when we discuss the extension of this theorem to meshes composed of rectangles or parallelepipeds. We have performed numerical experiments to compute the numerical CFL condition on nite meshes. We have observed (see section 6) that this condition was equivalent to condition [START_REF] Alvarez | A discontinuous nite element formulation for helmholtz equation[END_REF], except for a small range of value of α. Therefore we need a more restrictive necessary condition provided by the following theorem.

Theorem 3.3. Let V p,α = λ ∈ R : Q p,α (λ) = 0 and | Qp,α (λ) | ≤ 1 where Q p,α (λ) is a polynomial of degree
2p and Qp,α (λ) is a rational function whose expressions are given in appendix B. Then, scheme ( 5) is L 2 -stable only if, for p ≤ 5, the conditions ( 6) and ( 7) are satised and

             V p,α = ∅ or √ d c∆t h ≤ C 3,p (α) = 1 2 min V p,α . (8)
Remark 3.4. a

• This theorem does not provide an explicit CFL condition. However, it can be computed numerically by the following algorithm:

(1) Compute all the roots of

Q p,α , (2) Select the real roots such that | Qp,α (λ) | ≤ 1, (3)
Choose the maximum of these roots. • The numerical results in section 6 show that this theorem gives in practical cases necessary and sucient conditions.

• The numerical study of condition (8) that we present in section 6 shows that the set V p,α is actually empty except when α belongs to a small segment around α p 1 . This means that theorem 3.1 provides a sucient and necessary stability condition when α is not in this segment. Moreover the remarks 3.2 are still valid. We were unfortunately unable to establish this theorem for any p and we have restricted ourselves to p ≤ 5.

The proofs in the one dimensional case are given in section 4 while the extension to d = 3 is the subject of section 5. The proof for d = 2 can be easily deduced from the case d = 3.

Study in the 1-Dimensional case

This section contains the proofs of theorems 3.1 and 3.3 in the one dimensional case. It consists of three steps. The rst step is a Fourier analysis presented in section 4.1; the second step is devoted to the proof of conditions ( 6) and is presented in section 4.2; the last step concerns the proof of ( 7) and ( 8) in section 4.3. The proofs are detailed for p = 3 and are easily extendable to the cases p = 1, 2, 4 and 5.

Here, we assume that the domain is Ω = R and is meshed by segments of length h. We consider a velocity c 2 = µ /ρ = 1 but we can extend the proof to other velocities by setting ∆t = ∆t /c. We consider the scheme [START_REF] Baldassari | Modélisation et simulation numérique pour la migration terrestre par équation d'ondes[END_REF] without source term that is to say

M U n+1 -2U n + U n-1 ∆t 2 + KU n = 0. (9)
Considering the equation on one element J of the mesh, we have

∀J ∈ T h M 1,p U n+1 J -2U n J + U n-1 J ∆t 2 + K W 1,p T U n J-1 + K 1,p U n J + K W 1,p U n J+1 = 0 (10)
where U J corresponds to the vector of unknowns U restricted to the element J and M 1,p , K 1,p and K W 1,p are respectively the mass and stiness matrices in dimension 1 considering polynomials of degree p:

M 1,p (i, j) = h [0,1] φi (x) φj (x) dx, K 1,p (i, j) = 1 h [0,1] ∂ φi ∂ x (x) ∂ φj ∂ x (x) dx + 1 2h φi (1) ∂ φj ∂ x (1) + 1 2h φj (1) ∂ φi ∂ x (1) +γ φi (1) φj (1) - 1 2h φi (0) ∂ φj ∂ x (0) - 1 2h φj (0) ∂ φi ∂ x (0) +γ φi (0) φj (0) , K W 1,p (i, j) = - 1 2h φi (1) ∂ φj ∂ x (0) + 1 2h φj (0) ∂ φi ∂ x (1) -γ φi (1) φj (0) , (11)
where { φi } i=1,...,p+1 are the classical discontinuous Lagrange basis functions on the reference element [0, 1].

Fourier Analysis of the IPDG scheme in 1D

In order to study the stability of the IPDG scheme, we have to introduce the discrete Fourier transform

F h : L 2 h → L 2 (K h ) U → Ũ = F h (U ) (k) = h 2π J∈Z U J e -ikJh with K h = -π h , π h and L 2 h = U = (U J ) J∈Z , J∈Z U J 2 < +∞ .
Now, applying this discrete Fourier transform to [START_REF] Cohen | Higher-order nite elements with mass-lumping for the 1d wave equation[END_REF], we obtain

, ∀β ∈ [-π, π] M 1,p Ũ n+1 J (β) -2 Ũ n J (β) + Ũ n-1 J (β) ∆t 2 + K β Ũ n J (β) = 0 (12)
where β = hk and [START_REF] Fauqueux | Eléments nis mixtes spectraux et couches absorbantes parfaitement adaptées pour la propagation d'ondes élastiques en régime transitoire[END_REF] for all β ∈ [-π, π], is equivalent to the L 2 -stability of (9), thanks to the Parseval equalities. Since M 1,p is positive denite and K β is hermitian, all the eigenvalues of N β = M -1 1,p K β are real. Moreover, a classical stability analysis shows that ( 12) is stable if and only if

K β = K W 1,p T e -iβ + K 1,p + K W 1,p e iβ . The L 2 -stability of
0 ≤ λ ≤ 4 ∆t 2
for all λ ∈ Λ (β), where Λ (β) denotes the set of the eigenvalues of N β . Then, a necessary and sucient condition of the stability of ( 9) is

λ min ≥ 0 and ∆t ≤ 2 √ λ max with λ min = min β∈[-π,π] [min (Λ (β))] and λ max = max β∈[-π,π] [max (Λ (β))].
In section 4.2, we show that the condition λ min ≥ 0 is equivalent to [START_REF] Bassi | A high-order accurate discontinuous nite element method for the numerical solution of the compressible navier-stokes equations[END_REF] and in section 4.3, we show that the condition λ max ≤ 4 ∆t 2 implies ( 7) and (8). In the following, we consider the change of variable α = hγ to simplify the presentation. To show the equivalence between (3.1) and λ min ≥ 0, we have to consider the characteristic polynomial of N β :

q α (λ, β) = (-1) p+1 λ p+1 + p i=0 c i (α, β) λ i .
The coecients c i (α, β) can be computed by a symbolic calculus software such as Maple. We present them in appendix A for 1 ≤ p ≤ 5.

In order to study the sign of the eigenvalues of N β , we will use the following lemma.

Lemma 4.1. Let P be a polynomial of degree n with n real roots such that P (X) = n i=0 c i X i . All the roots of P are non negative if and only if

(-1) i c i ≥ 0.
Proof. The proof of this lemma can be found in [START_REF] Agut | Stability analysis of the interior penalty discontinuous galerkin method for the wave equation[END_REF].

Hence, we have to nd a condition over α such that, ∀i ∈ {0, . . . , p}, ∀β ∈ [-π, π],

(-1) i c i (α, β) ≥ 0.
Here we only detail the computations in the case p = 3 and we give the expression of the characteristic polynomials for p = 3 in the appendix A. For p = 3, we have

                       c 3 (α, β) = 8 h 2 ((15 -α) cos (β) -4α) c 2 (α, β) = 240 h 4 cos 2 (β) -(23 + α) cos (β) + (18α -65) c 1 (α, β) = 2880 h 6 4 cos 2 (β) + (65 -3α) cos (β) + (141 -32α) c 0 (α, β) = 100800 h 8 3 cos 2 (β) + 2 (3 -α) cos (β) + (2α -9) .
• Let us rst study the condition over c 3 . We have

, ∀β ∈ [-π, π], -c 3 (α, β) ≥ 0 ⇔ (α -15) cos (β) + 4α ≥ 0.
It is clear that this condition is satised for all β if and only if

(α -15) + 4α ≥ 0, ( 15
-α) + 4α ≥ 0 (13)
which implies that

α ≥ 3, α ≥ -5. (14)
Consequently, -c 3 (α, β) ≥ 0, whatever the choice of β, if and only if

α ≥ 3 (15)
• Let us now consider the condition over c 2 . Setting X = cos (β), this condition is equivalent to

f α (X) := X 2 -(23 + α) X + (18α -65) ≥ 0, ∀X ∈ [-1, 1]. (16)
This second order polynomial admits two roots:

       X 1 = 1 2 23 + α + (α -13) 2 + 620 1 /2 , X 2 = 1 2 23 + α -(α -13) 2 + 620 1 /2 .
We know that f α (X) is a second-order polynomial on the variable X and its head coecient is positive. Thus, to have f α non negative ∀X ∈ [-1; 1], we need one of the following conditions:

(1) the two roots are in ] -∞; -1] i.e.

X 1 ≤ -1, ( 2
) the two roots are in ]1; +∞[ i.e. X 2 ≥ 1, ( 3
) X 1 = X 2 .
Since (α -13) 2 + 620 > 0, X 1 > X 2 and 3. is impossible.

The case X 1 ≤ -1 is also impossible since X 1 ≥ 0 when α ≥ 0 so we just have to consider the case X 2 ≥ 1, which leads to the inequality

23 + α -(α -13) 2 + 620 1 /2 ≥ 2,
which is equivalent to α ≥ 87 17 .

(

) 17
Finally, c 2 (α, β) ≥ 0, whatever the choice of β, if and only if [START_REF] Shahbazi | An explicit expression for the penalty parameter of the interior penalty method[END_REF] holds. • Now, let us study the sign of c 1 (α, β).

Using the change of variable X = cos (β), the condition

-c 1 (α, β) ≥ 0, ∀β ∈ [-π, π] is equivalent to f α (X) := -4X 2 + (3α -65) X + (32α -141) ≥ 0, ∀X ∈ [-1, 1].
The polynomial f α admits the two following roots:

               X 1 = -1 8
  65 -3α + 3α + 61 3 2 + 14000 9 1 /2   , X 2 = -1 8   65 -3α - 3α + 61 3 2 + 14000 9 1 /2   .
Since, the head coecient of the polynom f α is negative, we need X 1 ≤ -1 and X 2 ≥ 1.

The condition X 1 ≤ -1 implies that

In the same way, the condition X 2 ≥ 1 is equivalent to

α ≥ 6. (19)
Consequently, -c 1 (α, β) ≥ 0, ∀β ∈ [-π, π] if and only if α ≥ 6. • Finally, let us look at the positivity of c 0 (α, β) , ∀β ∈ [-π, π].

Once again, using the change of variable X = cos (β), we have

f α (X) := 3X 2 + 2 (3 -α) X + 2α -9 ≥ 0.
This polynomial function f α admits the two following roots:

     X 1 = 1 3 (2α -9) , X 2 = 1.
In the same way than previously, we need X 1 ≥ 1 which leads to the condition

α ≥ 6. (20)
In conclusion, taking into account the conditions ( 15), ( 17), ( 18), ( 19) and (20), we have

λ min ≥ 0 ⇔ α ≥ 6. (21)
We used the same technique to derive a condition over α for all polynomial degrees p from p = 1 to p = 5. Since the calculations are very similar, we do not detail them here and we just present the conditions in Tab. 1.

p c 0 c 1 c 2 c 3 c 4 c 5 1 α ≥ 1 α ≥ 1 2 α ≥ 3 α ≥ 30 11 α ≥ 2 3 α ≥ 6 α ≥ 6 α ≥ 87 17 α ≥ 3 4 α ≥
p α (λ (β 0 ) , β 0 ) = 0 and ∂p α ∂λ (λ (β 0 ) , β 0 ) = 0.
However, the analysis of these points is extremely complicated and we have restricted ourselves to the study of the value β 0 such that λ (β 0 ) = 0. Therefore, we cannot arm that our condition is a sucient condition of stability. Nevertheless, the numerical results presented in section 6 show that this condition is sucient for numerical experiments. Using the implicit functions theorem, we have

∂p α ∂λ λ (β) + ∂p α ∂β = 0.
Hence β 0 satises ∂p α ∂β (λ (β 0 ) , β 0 ) = 0 and the steps of the proof are, for a given α:

(1) Find all the β 0 satisfying

∂p α ∂β (λ (β 0 ) , β 0 ) = 0;
(2) For each β 0 , compute λ max,β0 = max Λ (β 0 );

(3) Find the maximum value among all the λ max,β0 . In the following, we only detail the case p = 3, since the proofs are similar for all p = 1, . . . , 5.

We have

∂p α ∂β (λ, β) = sin (β) 8 h 2 (α -15) λ 3 + 240 h 4 (23 + α -2 cos (β)) λ 2 + 2880 h 6 (3α -65 -8 cos (β)) λ + 100800 h 8 (-6 cos (β) + 2 (α -3)) := sin (β) pα (λ, β) .
Consequently,

∂p α ∂β (β 0 ) = 0 ⇔        sin (β 0 ) = 0 or pα (λ (β 0 ) , β 0 ) = 0. (22)
First, we consider only the condition sin (β 0 ) = 0 and we show that it is equivalent to [START_REF] Alvarez | A discontinuous nite element formulation for helmholtz equation[END_REF]. Then, to obtain the necessary and sucient condition (8), we take into account the condition pα (λ (β 0 ) , β 0 ) = 0.

• If β 0 = 0, we obtain the following eigenvalues:

0; 60 h 2 ; 90 + 20α + 2g 1 (α) h 2 ; 90 + 20α -2g 1 (α) h 2
where g 1 (α) = 4545 -1320α + 100α 2 1 2 . It is clear that, for α ≥ 0, the two greatest eigenvalues are λ 1 = 60

h 2 and λ 2 = 1 h 2 (90 + 20α + 2g 1 (α)).
Studying the sign of the quantity

h 2 (λ 1 -λ 2 ) = -150 + 20α + 2g 1 (α)
we can easily obtain

λ max,0 = λ 2 if α ≥ 6 λ max,0 = λ 1 if α < 6.
(

) 23
We have proved in section 4.2 that the condition α ≥ 6 is a necessary stability condition. Therefore, we only have to consider λ max,0 = λ 2 .

• If β 0 = π, we obtain the following eigenvalues:

2 h 2 45 + √ 1605 ; 2 h 2 45 - √ 1605 ; 2 h 2 (-15 + 6α + g 2 (α)) ; 2 h 2 (-15 + 6α -g 2 (α))
where g 2 (α) = 405 -240α + 36α (

) 24
Now we have to compare λ max,0 and λ max,π . We easily verify that

       λ max,0 ≥ λ max,π if α ≥ 2 √ 1605 + 393 49 9.66, λ max,π > λ max,0 if 2 √ 1605 + 393 49 < α. (25)
Consequently, considering ( 23) and ( 25), a necessary condition of stability is

∆t h ≤        2 λ max,π if 6 ≤ α ≤ α 1,p , 2 λ max,0 if α 1,p < α, ( 26
)
with α 1,p = 2 √ 1605 + 393 49
, which corresponds to the necessary condition (6). • Let us now nd β 0 such that pα (λ (β 0 ) , β 0 ) = 0. In fact, we don't have to compute β 0 , we interest ourselves only to max Λ (β 0 ). We can easily obtain that

cos (β 0 ) = (α -15) h 6 λ 3 (β 0 ) + 30 (23 + α) h 4 λ 2 (β 0 ) + 360 (3α -65) h 2 λ (β 0 ) + 25200 (α -3) 60 (h 4 λ 2 (β 0 ) + 48h 2 λ (β 0 ) + 1260) := Qα (λ (β 0 )) .
Using this expression of cos (β 0 ) in the characteristic polynomial (35) we obtain that λ = λ (β 0 ) is solution to

q α (λ, β 0 ) = - 1 15h 8 (h 4 λ 2 + 48h 2 λ + 1260) 6 i=0 λ i h 2i ci (α) = 0
or, equivalently, solution to

Q p,α (λ) = 6 i=0 λ i h 2i ci (α) = 0 with                                      c 0 (α) = 635040000 α 2 -12α + 36 , c 1 (α) = 3628800 15α 2 + 70α -96 , c 2 (α) = 86400 31α 2 -447α + 5316 , c 3 (α) = 14400 8α 2 -135α -1728 , c 4 (α) = 180 17α 2 -442α + 7740 , c 5 (α) = 60 α 2 + 16α -357 , c 6 (α) = α 2 -30α + 210.
After having computed the roots of Q p,α , we have to verify that β 0 is well dened, that is to say

| cos (β 0 ) | = | Qα (λ) | ≤ 1.
That is why we are interested uniquely in the eigenvalues λ verifying Q α (λ) = 0 and | Qα (λ) | ≤ 1.

The d-dimensional case

In this section, we propose to adapt the technique proposed in [START_REF] Fauqueux | Eléments nis mixtes spectraux et couches absorbantes parfaitement adaptées pour la propagation d'ondes élastiques en régime transitoire[END_REF] to extend the analysis from the 1D case to the dD case. Here, we only detail the three-dimensional case, since the technique is exactly the same for the two dimensional case. First of all, we consider an innite homogeneous 3D domain Ω uniformly meshed by cubes of edge h. We introduce the following notations.

• Ω = K J ∈T h K J where K J = 3 k=1 S J k = 3 k=1
[J k h, (J k + 1) h] and J = (J k ) k=1,...,3 .

• On K [0, 1] d , we dene the Lagrange basis functions ( φl ) l∈{1,...,p+1} 3 by

φl (x) = 3 k=1 φl k (x k )
where φl k are the 1D Lagrange basis functions.

• Since the mesh is uniform, the basis functions are dened thanks to the functions ( φl ) l∈{1,...,p+1} 3 by

ϕ J m (x) = φm x -Jh h 1 K j (x)
where 1 K J is the indicator function of K J . These functions can be written as a product of d 1D basis functions:

ϕ J m (x) = 3 k=1 φm k x -J k h h 1 S J k (x k ) .
• The dierent faces of the reference element K are denoted by an exponent C corresponding to the orientation of the face: North (N), South (S), East (E), West (W), in Front of (F) and in the Back (B) (cf. Since the mesh is uniform, we can rewrite the problem on an element I = {I 1 , I 2 , I 3 }:

M 3,p δ n U I1,I2,I3 = K 3,p U I1,I2,I3 + K E 3,p U I1+1,I2,I3 + K W 3,p U I1-1,I2,I3 +K F 3,p U I1,I2+1,I3 + K B 3,p U I1,I2-1,I3 + K N 3,p U I1,I2,I3+1 +K S 3,p U I1,I2,I3-1 (27)
where

• U I1,I2,I3 corresponds to the restriction of U on the element I.

• δ n U I1,I2,I3 = U n+1 I1,I2,I3 -2U n I1,I2,I3 + U n-1 I1,I2,I3 ∆t 2 • M 3,p is a block of the mass matrix M , M 3,p (i, j) = h 3 K φi φj dx, i, j ∈ {1, . . . , p + 1} 3 • K 3,p is a diagonal block of the matrix K K 3,p (i, j) = h K ∇ φi • ∇ φj dx - C∈{N,S,E,W,B,F } h 2 Γ C ( φi ∇ φj + φj ∇ φi ) ν C dσ + C∈{N,S,E,W,B,F } h 2 Γ C γ φi φj dσ, i, j ∈ {1, . . . , p + 1} 3
where ν C is the outward unit normal vector to the face Γ C .

• K C 3,p is a block of the matrix K corresponding to the interactions between an element I and its neighbour on the face Γ C :

K E 3,p (i, j) = [0,1] 2 h 2 ( φi (1, x 2 , x 3 ) ∇ φj (0, x 2 , x 3 ) + φj (0, x 2 , x 3 ) ∇ φi (1, x 2 , x 3 )) ν E -h 2 φi (1, x 2 , x 3 ) φj (0, x 2 , x 3 ) dx 2 dx 3 K F 3,p (i, j) = [0,1] 2 h 2 ( φi (x 1 , 1, x 3 ) ∇ φj (x 1 , 0, x 3 ) + φj (x 1 , 0, x 3 ) ∇ φi (x 1 , 1, x 3 )) ν F -h 2 φi (x 1 , 1, x 3 ) φj (x 1 , 0, x 3 ) dx 1 dx 3 K N 3,p (i, j) = [0,1] 2 h 2 ( φi (x 1 , x 2 , 1) ∇ φj (x 1 , x 2 , 0) + φj (x 1 , x 2 , 0) ∇ φi (x 1 , x 2 , 1)) ν N -h 2 φi (x 1 , x 2 , 1) φj (x 1 , x 2 , 0) dx 1 dx 2 K W 3,p ((i 1 , i 2 , i 3 ) , (j 1 , j 2 , j 3 )) = K E 3,p ((j 1 , i 2 , i 3 ) , (i 1 , j 2 , j 3 )) K B 3,p ((i 1 , i 2 , i 3 ) , (j 1 , j 2 , j 3 )) = K F 3,p ((i 1 , j 2 , i 3 ) , (j 1 , i 2 , j 3 )) K S 3,p ((i 1 , i 2 , i 3 ) , (j 1 , j 2 , j 3 )) = K N 3,p ((i 1 , i 2 , j 3 ) , (j 1 , j 2 , i 3 ))
Then, multiplying the equation ( 27) by the inverse of the mass matrix M 3,p , we obtain

δ n U I1,I2,I3 = N 3,p U I1,I2,I3 + N E 3,p U I1+1,I2,I3 + N W 3,p U I1-1,I2,I3 +N F 3,p U I1,I2+1,I3 + N B 3,p U I1,I2-1,I3 + N N 3,p U I1,I2,I3+1 +N S 3,p U I1,I2,I3-1 (28)
where

N 3,p = M -1 3,p K 3,p and N C 3,p = M -1 3,p K C 3,p .
Now, we are interested in rewriting the matrices N 3,p and N C 3,p with respect to the matrices we have obtained for the one dimensional case.

From the 3 dimensional case to the one dimensional case

The coecients of M 3,p , K 3,p , K C 3,p , N 3,p and N C 3,p can be deduced from the coecients of M 1,p , K 1,p , K W 1,p , N 1,p and N W 1,p thanks to the following theorem.

TITLE WILL BE SET BY THE PUBLISHER Theorem 5.1. For all m = (m k ) k=1,...,3 ∈ {1, . . . , p + 1} 3 and n = (n k ) k=1,...,3 ∈ {1, . . . , p + 1} 3 , we have

1. M 3,p (m, n) = 3 i=1 M 1,p (m i , n i ) , 2. K 3,p (m, n) = 3 i=1   K 1,p (m i , n i ) 3 k=1, k =i M 1,p (m k , n k )   , 3. K C 3,p (m, n) = K W 1,p (m p C , n p C ) 3 k=1, k =p C M 1,p (m k , n k ) , 4. N 3,p (m, n) = 3 p=1   N 1,p (m p , n p ) 3 k=1, k =p δ m k ,n k   , 5. N C 3,p (m, n) = N W 1,p (m p C , n p C ) 3 k=1, k =p C δ m k ,n k , (29)
where

p C =          1 if C ∈ {E, W } , 2 if C ∈ {N, S} , 3 if C ∈ {B, F } ,
and

N 1,p = M -1 1,p K 1,p and N C 1,p = M -1 1,p K C 1,p .
The proof of this theorem is given in appendix C.

Consequences on the stability analysis

Let us now apply a Fourier transform in the three directions to (28) to obtain, for β = [-π, π] 3 ,

δ n Ũβ1,β2,β3 = N β Ũβ1,β2,β3 (30)
where the matrix N β is dened by

N β (m, n) = 3 j=1   N 1,p (m j , n j ) 3 q=1, q =j δ mq,nq + e iβj N W 1,p (m j , n j ) 3 q=1, q =j δ mq,nq + e -iβj N W 1,p (n j , m j ) 3 q=1, q =j δ mq,nq  
which can be rewritten as

N β (m, n) = 3 j=1   N 1,p (m j , n j ) + e iβj N W 1,p (m j , n j ) + e -iβj N W 1,p (n j , m j ) 3 q=1, q =j δ mq,nq   = 3 j=1 N βj (m j , n j ) 3 q=1, q =j δ mq,nq .
Using the stability analysis as in section 4.1, the stability of the scheme is ensured if and only if To compute these values, we use the following lemma. Lemma 5.2. Let λ β i i=1,...,p+1 denote the eigenvalues of N β and v β i i=1,...,p+1 the associated eigenvectors. Then, the eigenvalues of N β are given by

λ β i = 3 k=1 λ β k i k
and the associated eigenvectors by

v β i (m) = 3 k=1 v β k i k (m k ) . (31)
Proof. Let v β i dened by (31), then

N β v β i (m) = p+1 n1,...,n3=1   3 q=1 N βq (m q , n q ) 3 k=1, k =q δ m k ,n k 3 k=1 v β k i k (n k )   = 3 q=1     p+1 nq=1 N βq (m q , n q ) v βq iq (n q )   3 k=1, k =q v β k i k (m k )   = 3 q=1   λ βq iq v βq iq (m q ) 3 k=1, k =q v β k i k (m k )   = 3 q=1 λ βq iq v β i (m)
It is then clear that λ min,3 = 3λ min and λ max,3 = 3λ max .

Hence, the scheme (31) is stable, if and only if

λ min ≥ 0 and c∆t h ≤ 1 √ 3 2 λ max .
The rst condition is equivalent to condition [START_REF] Bassi | A high-order accurate discontinuous nite element method for the numerical solution of the compressible navier-stokes equations[END_REF], while the second one implies ( 7) and (8).

Extension to rectangular or parallelepiped mesh

For the sake of simplicity, we restricted our theorem to the case of squared or cubic mesh. However, one can extend the proof to the case of rectangular or parallelepipeds meshes to show that a necessary stability condition is in 2D:

γ ≥ p (p + 1) 2 min (h x , h y )
and in 3D:

γ ≥ p (p + 1) 2 min (h x , h y , h z ) .
Here, h x , h y and h z denote respectively the length of the edges of the elements in the x, y and z direction. The minimal value of h x , h y and h z is actually the diameter of the inscribed sphere of each element. This remark conrms that the third denition of ξ F using the diameter of the inscribed sphere or circle in 2D is the most appropriate.

The proof could also be extended to obtain a CFL condition, but its expression is complicated and does not add much insight.

Numerical results

In this section, we rst represent the behaviour of the CFL condition with respect to α and we show that the set V p,α is empty for almost all the values of α (section 6.1). This illustrate the fact that theorem 3.1 is actually necessary and sucient for most of the values of α. Then, we compare our analytical CFL condition in innite domain to the CFL condition computed numerically on nite meshes in order to illustrate the validity of theorem 3.3 (section 6.2). We computed numerically the greatest eigenvalue λ max of the matrix M -1 K and we deduced the CFL condition of the scheme using the formula In this paper, we have proposed necessary conditions of L 2 -stability of an IPDG method using regular meshes and the numerical results showed that these conditions are actually sucient in practice. It also conrm the conjecture of Ainsworth, Monk and Muniz up to p = 5. Moreover, we have observed that the CFL condition is constant with respect to α on a segment

p (p + 1) 2
, α and is decreasing as α -1 /2 for α > α. This means that it is not necessary to choose α too close to p (p + 1) 2

to improve the CFL condition. Finally, we have observed that a good choice for ξ F should be the diameter of the inscribed circle (or sphere). This should be conrmed by an analysis on triangular meshes, which will be the topic of a future work.

                               c 4 (α, β) = - 10 h 2 ((α -24) cos (β) -5α) c 3 (α, β) = - 120 h 4 5 cos 2 (β) + (4α + 287) cos (β) + 10 (15α -88) α c 2 (α, β) = - 10080 h 6 5 cos 2 (β) + (3α -305) cos (β) + 990 -133α c 1 (α, β) = - 201600 h 8 15 cos 2 (β) + (8α + 165) cos (β) + 2 (59α -468) c 0 (α, β) = - 50803200 h 10 2 cos 2 (β) + (α -10) cos (β) + 8 -α .
• In the case of discontinuous basis functions of degree 5, the characteristic polyomial associated to the matrix A β is

p α (λ, β) = λ 6 + c 5 (α, β) λ 5 + c 4 (α, β) λ 4 + c 3 (α, β) λ 3 + c 2 (α, β) λ 2 + c 1 (α, β) λ + c 0 (α, β) (35)
with We present here the expressions of the polynomial Q p,α and the rational function Qp,α for 1 ≤ p ≤ 5.

                                       c 5 (α, β) = - 12 h 2 ((α -35) cos (β) + 6α) c 4 (α, β) = - 420 h 4 -3 cos 2 (β) + (2α + 336) cos (β) + 1155 -134α c 3 (α, β) = - 40320
• For the polynomials of degree 1, Qp,α is dened by

Qp,α (λ) = h 2 λ 2 α 3 -1 + α.
We have, • Proof of 1.

Considering the notations and the results of the section 5, we have • Proof of 2. We rst have the following lemma for the volumic term.

Lemma C.2. For all m = (m k ) k=1,...,3 ∈ {1, . . . , p + 1} 3 and n = (n k ) k=1,...,3 ∈ {1, . . . , p + 1} 3 , we have

h K ∇ φm • ∇ φn dx = 3 i=1   1 h [0,1] ∂ φmi ∂x i (x i ) ∂ φni ∂x i (x i ) dx i 3 k=1, k =i M 1,p (m k , n k )   .
Proof. We know that, ∀m ∈ {1, . . . , p + 1} d φm (x) = which implies that

∂ φm ∂x k (x) = ∂ φm k ∂x k (x k ) 3 i=1, i =k φmi (x i ) .
Then, using the same reasoning as previously, Proof. First of all, using (37), we have In the same way, for the penalization term, we have: which ends the proof.

h K ∇ φm • ∇ φn dx = 3 i=1   1 h [0,1] ∂ φmi ∂x i (x i ) ∂ φni ∂x i (x i ) dx i 3 k=1, k =i h [0,1] φm k (x k ) φn k (x k ) dx k   = 3 i=1   1 h [0,1] ∂ φmi ∂x i (x i ) ∂ φni ∂x i (x i ) dx i 3 k=1, k =i

4. 2 .

 2 Study of the condition λ min ≥ 0

  Fig 1).

Figure 1 .

 1 Figure 1. Notations of the faces in 3D

λ

  min,3 ≥ 0 and λ max,3 ≤ 4 ∆t 2 where λ min,3 = min β∈[-π,π] 3 (min Λ (N β )), λ max,3 = max β∈[-π,π] 3 (max Λ (N β )) and Λ (N β ) is the set of eigenvalues of N β .

6. 1 .

 1 Behaviour of the CFL condition with respect to α In Fig. 2, 4, 6, 8 and 10 we plot the functions C 1,p (blue line with diamonds), C 2,p (α) (red line with circles) and C 3,p (α) (black line) respectively for p = 1, 2, 3, 4 and 5. The function C 3,p (α) only modies the CFL condition in a small segment around α 1,p . The behaviour is conrmed in Fig. 3, 5, 7, 9 and 11 which represent a zoom around α 1,p . Theses numerical results conrm the fact that theorem 3.1 provides actually a necessary and sucient condition except for a small range of α. Moreover, the CFL condition remains constant for α from p (p + 1) 2 to a close value of α 1 p , which means that it is not necessary to choose α = p (p + 1) 2to optimize the time step.

Figure 2 .

 2 Figure 2. The 3 conditions for p = 1

Figure 3 .

 3 Figure 3. Zoom on the 3 cond. for p = 1

Figure 4 .

 4 Figure 4. The 3 conditions for p = 2

Figure 5 .

 5 Figure 5. Zoom on the 3 cond. for p = 2

Figure 6 .

 6 Figure 6. The 3 conditions for p = 3

Figure 7 .

 7 Figure 7. Zoom on the 3 cond. for p = 3

Figure 8 .

 8 Figure 8. The 3 conditions for p = 4

Figure 9 .

 9 Figure 9. Zoom on the 3 cond. for p = 4

Figure 10 .

 10 Figure 10. The 3 conditions for p = 5

Figure 11 . 5 6. 2 .

 1152 Figure 11. Zoom on the 3 cond. for p = 5

.

  In Fig.12, 13, 14, 15 and 16 we compare the analytical CFL (red line) obtained by theorem 3.3 to the numerical CFL (triangles), respectively for p = 1, 2, 3, 4 and 5. All gures show a very good agreement between the analytical and the numerical CFL.

Figure 12 .

 12 Figure 12. Numerical comparison in P 1

Figure 13 .

 13 Figure 13. Numerical comparison in P 2

Figure 14 . 3 Figure 15 .

 14315 Figure 14. Numerical comparison in P 3

Figure 16 .

 16 Figure 16. Numerical comparison in P 5

  with

h 6 - 4 8 - 9

 6489 cos 2 (β) + (2α -702) cos (β) + 256α -2849α c 2 (α, β) = -1814400 h cos 2 (β) + (4α770) cos (β) + 9343 -326α c 1 (α, β) = -101606400 h 10 -12 cos 2 (β) + (5α -303) cos (β) + 94α -1170 c 0 (α, β) = -10059033600 h 12 -5 cos 2 (β) + (2α -20) cos (β) + 25 -2α . Appendix B. Definition of Q p,α and Qp,α

α 2 -

 2 2α + 1 , c1 (α) = 12 α 2 -α , c2 (α) = α 2 -6α + 6,

M 3 , 3 i=1M 1

 331 p (m, n) = h 3 K φm φn dx i ) φni (x i ) dx i = ,p (m i , n i )

M 1 ,

 1 p (m k , n k )  which ends the proof. Now, we have to deal with the surface terms. Let us rst remark that, over all the faces Γ C ,∇ φm • ν| Γ C = ∂ φmp C ∂x p C (x p C ) ν 1,C 3 k=1, k =p C φm k (x k )(37)where ν 1,C is the outward unit normal vector in the one dimensional case dened byν 1,C =    1 if C ∈ {E, N, F } , -1 if C ∈ {W, S, B}and x p C is dened byx p C =    1 if C ∈ {E, N, F } , 0 if C ∈ {W,S, B} . Then, we can propose the following lemma. Lemma C.3. For all C ∈ {E, W, N, S, F, B}, we have ∀m, n ∈ {1, . . . , p + 1} d

2 Γ

 2 h Γ C φm (∇ φn • ν) dσ = 1 h φmp C (x p C ) ∂ φnp C ∂x p C (x p C ) ν 1,C 3 k=1, k =p C M 1,p (m k , n k ) h C γ φm φn dσ = γ φmp C (x p C ) φnp C (x p C ) 3 k=1, k =p C M 1,p (m k , n k ) .

3 k=1, k =p C h 2 [0, 1 ]

 321 x p C ) ∂ φnp C ∂x p C (x p C ) ν 1,C φm k (x k ) φn k (x k ) dx kwhich can be rewritten asΓ C φm (∇ φn • ν) dσ = φmp C (x p C ) ∂ φnp C ∂x p C (x p C ) ν 1,C 3 k=1, k =p C M 1,p (m k , n k ) .

γ 2 Γ C γ 3 k=1, k =p C φm k (x k ) 3 k=1, k =p C φn k (x k ) 3 k=1, k =p C dσ k = 3 k=1, k =p C h 2 [0, 1 ]

 2333k321 φm φn dσ = φmp C (x p C ) φnp C (x p C ) h γ φmp C (x p C ) φnp C (x p C ) φm k (x k ) φn k (x k ) dx k which clearly implies that Γ C γ φm φn dσ = γ φmp C (x p C ) φnp C (x p C ) 3 k=1, k =p C M 1,p (m k , n k )

Table 1 .

 1 Conditions over α for each coecient c i and each polynomial degree p From these results, we can easily deduce the smallest penalization parameters ensuring the stability of the scheme (see Tab. 2). It is clear that, for 1 ≤ p ≤ 5, the stability is guaranteed if and only if

	10 α ≥	543 55	α ≥	325 34	α ≥	581 73	α ≥ 4
	5 α ≥ 15 α ≥ 15 α ≥	336 23	α ≥	1185 86	α ≥	124 11	α ≥ 5
	α ≥	p (p + 1) 2	, or, equivalently, γ ≥	p (p + 1) 2h	.

α ≥ 1 α ≥ 3 α ≥ 6 α ≥ 10 α ≥ 15

Table 2 .

 2 Stability condition over α for each polynomial degree p

	4.3. The CFL condition
	Now, we propose to prove that the condition
	∆t ≤	√	2 λ max

implies

[START_REF] Alvarez | A discontinuous nite element formulation for helmholtz equation[END_REF]

and

[START_REF] Cohen | Higher-Order Numerical Methods for Transient Wave Equations[END_REF]

. To compute λ max , we use the implicit function theorem. Indeed, each eigenvalue λ ∈ Λ (β) is an implicit function λ (β) dened by p α (λ (β) , β) = 0. Since p α (λ, β) is a continuous and periodic function of β, so is λ (β) and we only need to nd all the values β 0 ∈ [-π; π] such that λ (β 0 ) = 0 or such that λ (β 0 ) is not dierentiable in β = β 0 . These latter points check

√α

/2 = 1 ms -1 . We impose also Dirichlet boundary conditions at the both ends of the domain and the length of the space step is h = 0.1.

Appendix A. Expression of the polynomial p α • In the case of discontinuous basis functions of degree 1, we can easily obtain the following characteristic polyomial associated to the matrix

with

c 0 (α, β) = 12 h 4 cos 2 (β) -2α cos (β) + 2α -1 .

• In the case of discontinuous basis functions of degree 2, the characteristic polyomial associated to the matrix A β is

with

• In the case of discontinuous basis functions of degree 4, the characteristic polyomial associated to the matrix A β is

• In the case p = 2, the denition of Qp,α is

The polynomial Q p,α is such that

• For p = 3, we have

and Qp,α is dened by

• For the polynomials of degree 4, Qp,α is dened by

And, we have

with

• For p = 5, Qp,α is such that

and we have

This section is devoted to the proof of the following theorem Theorem C.1. For all m = (m k ) k=1,...,3 ∈ {1, . . . , p + 1} 3 and n = (n k ) k=1,...,3 ∈ {1, . . . , p + 1} 3 , we have

where

and

Proof. a Finally, using the two lemmas, we obtain

• Proof of 3.

To rewrite the terms K C 3,p (m, n) for all m, n, we use a similar reasoning as for K 3,p . • Proof of 4.

To prove 4. and 5., we need the following lemma.

Lemma C.4. Let m = (m p ) p=1,...,3 and n = (n p ) p=1,...,3 .

We have

Proof. Let A be the matrix dened ∀m, n ∈ {1, . . . , p + 1} d by

We have

where I is the identity matrix, which ends the proof.

Let us now consider the matrix N 3,p = M -1 3,p K 3,p . First, we rewrite K 3,p as

Then, M -1

3,p T q and, using lemma C.4

So that

Performing the same calculations for T 2 and T 3 , we obtain

• Proof of 5.

We apply the technique used to prove 4. to show that


Andrea G B Tettamanzi
email: andrea.tettamanzi@unice.fr

Célia Da Costa Pereira

Testing Carlo Cipolla's Laws of Human Stupidity with Agent-Based Modeling

published or not. The documents may come

L'archive ouverte

I. INTRODUCTION

In a tongue-in-cheek essay [START_REF] Cipolla | The basic laws of human stupidity[END_REF] written in English, tainted with Swiftian humor, printed in 1976 in a hundred of copies and distributed to friends as a Christmas gift, the late Italian historical economist Carlo M. Cipolla put forward a theory of human stupidity, articulated in five fundamental laws, which is regarded by many as a real work of genius.

One should not be misled by the humorous tone of Cipolla's essay into thinking his theory cannot be taken seriously. In fact, in most cultures, humor and jokes are a way to tell truths that hurt without breaking social norms or sounding irrespectful.

If taken seriously, Cipolla's laws should enable one to make falsifiable claims. Therefore, it should be possible to test his theory experimentally. In particular, if one accepts Darwin's theory of evolution, then for Cipolla's theory of human stupidity to be accepted it must be compatible with it. In this paper, we investigate under which circumstances and based on which assumptions the two theories do not contradict each other. We use agent-based simulation as an experimental tool to this aim.

In the last decade, agent-based modeling has become a widely accepted tool for studying the dynamics of complex systems in the social sciences [START_REF] Bonabeau | Agent-based modeling: Methods and techniques for simulating human systems[END_REF], including economical phenomena like financial markets [START_REF] Ehrentreich | Agent-Based Modeling[END_REF]. Agent-based modeling is the computational study, carried out by means of simulation, of economies that are modeled as evolving systems of autonomous interacting agents [START_REF] Axelrod | The Complexity of Cooperation: Agent-Based Models of Competition and Collaboration[END_REF], [START_REF] Tesfatsion | Agent-based computational economics: modeling economies as complex adaptive systems[END_REF]. The main advantage of agent-based over analytical models is that they allow the removal of restrictive assumptions and can, thus, lead to a more realistic description of the system under investigation.

The paper is organized as follows: a summary of the main claims of Cipolla's theory of human stupidity is given in Section II; the theory is then critically questioned in Section III and some conjectures are put forward to solve the arising issues. Section IV describes an experimental setup using agent-based simulation designed to test these conjectures; the outcome of the experiments are reported in Section V and conclusions are drawn in Section VI.

I b I h H i H s S h S b B s B i

II. CIPOLLA'S THEORY OF HUMAN STUPIDITY

Cipolla identifies in stupid people one of the main obstacles to welfare in human societies and offers a coherent account of why stupidity is so powerful and hard to act against. To this end, he builds an abstract model of a human agent's social behavior, which may be summarized by two coordinates:

x the average gain (or loss) that an agent obtains as a result of his or her actions; y the average gain (or loss) that an agent produces to other agents or groups of agents.

Two things have to be stressed: first of all, this notion of gain or loss is subjective and includes psychological and emotional factors-it must be thus understood as a utility as used in economics or game theory [START_REF] Neumann | Theory of Games and Economic Behavior[END_REF]; secondly, x and y are statistical expectations, for the actions of a given agent may yield different outcomes in different circumstances and an agent's behavior is, therefore, to be described by means of a probablity distribution over the Cartesian plane.

As a result, human agents can be plotted as points on a diagram like the one shown in Figure 1 based on their x, y behavior. Figure 1 divides the two-dimensional plane into four quadrants or eight sectors, corresponding to different categories of agents.

Cipolla's definition of a stupid person (his Third Law) is "a person who causes losses to another person or to a group of persons while himself deriving no gain and even possibly incurring a loss", thus, technically, an agent with x ≤ 0 and y < 0. Stupid persons thus fall in quadrant S of Figure 1. The persons falling in quadrant I (x > 0, y ≥ 0) may be defined as intelligent, those falling in quadrant H (x ≤ 0, y ≥ 0) as helpless, and those falling in quadrant B (x > 0, y < 0) as bandits. Each quadrant may be further divided into two sectors: for instance, quadrant B may be divided into sector B i , containing bandits with undertones of intelligence (for they cause fewer losses to others than the gains they derive) and B s , containing bandits with undertones of stupidity (for their actions yield to them gains inferior to the losses inflicted to other people), and so on. Now, Cipolla observes that any numerical estimate of the fraction σ of stupid people always and inevitably turns out to be an underestimate (First Law) and that the probability that a given person be stupid is independent of any other characteristic of that person (Second Law). From these observations, we can infer that stupid people must be an overwhelming majority of any sample we may draw from the general population.

Cipolla's Fourth and Fifth Laws basically say that stupid persons are the most dangerous type of persons, because it is impossible for non-stupid people to organize defenses against them, due to their irrationality and unpredictability.

III. CRITIQUE

There is something counterintuitive in Cipolla's theory: individuals who do not derive any gain and even possibly incur a loss from their actions, like stupid and helpless people, should have a competitive disadvantage vis-à-vis more opportunistic individuals. One would thus expect that, in the long run, they would lose ground and eventually get extinct. This, at least, is what would be suggested by the Darwinian law of natural selection or "survival of the fittest", which is generally accepted as a well-corroborated explanation of the evolution of the species, but also the basis of recent economic models, like the adaptive market hypothesis [START_REF] Lo | The adaptive markets hypothesis[END_REF]. On the other hand, for the same reason, rational individuals like intelligent people and bandits should be expected to be more successful than stupid and helpless people, thus being able in the long run to take over the entire population.

The fact that evolution must be taken into account is indirectly confirmed by Cipolla himself, who writes that the potential for a stupid person to cause damages depends, first of all, on the genetic factor. He even goes as far as to postulate a gene for stupidity, which is inherited by stupid individuals from their parents. Why, then, despite their apparent evolutionary handicap, can stupid people make up an overwhelming fraction σ of the human race? They must have some sort of competitive advantage which makes them particularly fit for reproduction, but how can it be?

There are several conjectures we could make to attempt to answer this question: 1) stupid people, causing damages to all other members of the society, neutralize, as it were, the forces of natural selection; 2) stupid people are more resilient to damages inflicted by others;

3) the observed σ fraction is an effect of the initial distribution of agents in the population: although it may not always be the case that σ is overwhelmingly high, once for some reason the stupid become a large majority, the situation cannot be reversed and evolution remains trapped in that state; 4) etc.

In the next section, we describe an experiemental setup which should allow us to test these conjectures and suggest a solution to the enigma.

IV. EXPERIMENTAL SETUP

Members of the society are modeled as agents, whose behavior, stochastic in nature, is governed by a probability distribution which is innate in each agent and genetically determined.

The agents are evolutionary, in the sense that they are individuals of an evolutionary algorithm, whose details are illustrated below. Therefore, every agent has a genome, a set of parameters which determine its behavior. These genetic parameters do not change during the agent's entire lifetime but are subject to random mutation when passed on from a parent to its offspring. An agent's genome consists of the parameters of a twodimensional normal joint probability distribution,

(µ x , µ y , σ x , σ y , θ), (1)
where µ x , µ y is the mean, σ x is the standard deviation along the x-axis, σ y is the standard deviation along the y-axis, and θ is the rotation angle, which is used to randomly determine the outcome of its "actions".

Agents are subject to an evolutionary algorithm like the one described in [START_REF] Mauri | Investigating the role of social behavior in financial markets through agent-based simulation[END_REF].

Agents interact as follows:

1) an agent, whom we will call "active" because it will be the one that takes an action, is randomly selected from the population; 2) another agent, whom we will call "passive" because it will undergo the effects of the action taken by the active agent, is randomly selected from the remaining agents; 3) a pair x, y is randomly extracted from the active agent's normal distribution; 4) the active agent's wealth is updated according to x; 5) the passive agent's wealth is updated according to y.

The above cycle, called a period, which is the clock that measures time in a simulation, is repeated over and over. We will denote by t the current period in a simulation.

Notice that, in general, this is not a zero-sum game. If most agents act intelligently, the population will enjoy an overall wealth increase; if, on the other hand, most agents act stupidly, the overall welfare of the population will decrease and nothing prevents it from eventually becoming extinct.

Optionally, one might enforce a zero-sum game by redistributing net wealth surplus or loss proportionally to all the agents in a population to model the fact that in the real-world resources are finite (or, at least, quite unelastic) and the gain of some always determines a loss for the rest of the society.

Another possible variation of the model has to do with the way wealth is updated: let W a t and W p t be the wealth of the active and of the passive agent, respectively, at time t; then,

• a linear wealth trasfer is when W a t+1 = W a t + x and W p t+1 = W p t + y; • a logarithmic wealth transfer may be defined, whereby

W a t+1 = W a t + x, if x ≤ 0; W a t + log(x + 1), otherwise; (2)
and similarly for W p t+1 ; • a hyperbolic wealth transfer may be defined, whereby

W a t+1 = W a t + x, if x ≤ 0; W a t + x x+1 , otherwise; (3)
and similarly for W p t+1 .

Notice that both the logarithmic and the hyperbolic transfer functions introduce an asymmetry between losses, which remain linear in all cases, and gains, whose transfer gets harder as they increase. The key intuition these two definitions try to capture is that it is exceedingly easy to destroy wealth, but creating it is hard, and gets harder the more wealth one wants to create.

To model the fact that rational agents know better, we may then introduce the notion of defense, by stating that rational agents (i.e., agents such that µ x > 0) are able to build defenses against bandits, but not, according to Cipolla's Fourth and Fifth laws, against stupid agents, with whom any countermeasure would turn out to be vain. This is how we propose to model such defense: in an interaction, if the passive agent is rational and the active agent is behaving like a bandit (i.e., x > 0 and y < 0), both x and y are discounted by multiplying them by a "defense factor" 1 -δ, where

δ = µ p x µ p x + 1 (4)
represents a hyperbolic discount factor which grows with the rationality of the passive agent; an infinitely rational agent will be able to completely neutralize a bandit's actions; an agent who is just slightly rational will discount but a small percentage of the damage caused by a bandit.

Finally, to model the hypothesis that stupid agents are more resilient than others to damages inflicted by their peers, we introduce the possibility that the x and y effects of an interaction be "relativized" with respect to the µ x of the receiving agent (be it active or passive). Accordingly,

• the active agent's wealth would be updated according to x -µ a x ; • the passive agent's wealth would be updated according to y -µ p x .

A. Evolutionary Algorithm

Evolutionary algorithms (EAs) [START_REF] De | Evolutionary Computation: A unified approach[END_REF], [START_REF] Eiben | Introduction to Evolutionary Computing[END_REF] are a broad class of stochastic optimization algorithms, inspired by biology and in particular by those biological processes that allow populations of organisms to adapt to their surrounding environment: genetic inheritance and survival of the fittest. Each individual of the population represents a point in the space of the potential solutions for the considered problem. The evolution is obtained by iteratively applying a (usually quite small) set of stochastic operators, known as mutation, recombination, and selection. Mutation randomly perturbs a candidate solution; recombination decomposes two distinct solutions and then randomly mixes their parts to form novel solutions; and selection replicates the most successful solutions found in a population at a rate proportional to their relative quality. The initial population may be either a random sample of the solution space or may be seeded with solutions found by simple local search procedures, if these are available. The resulting process tends to find, given enough time, globally optimal solutions to the problem much in the same way as in nature populations of organisms tend to adapt to their surrounding environment.

An evolutionary algorithm is used by the simulator to make the agents evolve according to their interactions with their peers.

1) Representation:

The parameters listed in Equation 1, which make up the genome of an agent, are encoded as an array of five floating-point numbers.

The bivariate normal probability density function may be written

f (x, y) = 1 2π |C -1 | e (x-µx,y-µy)C(x-µx,y-µy) T , (5)
where the matrix

C = a b b c
is positive-definite. Equation 5 may thus be rewritten as

f (x, y) = √ ac -b 2 2π • e -ar 2 x -2brxry-cr 2 y , (6)
where r x = x -µ x and r y = y -µ y .

Instead of encoding parameters a, b, and c directly in the genome, we encode the angle θ and the standard deviations σ x and σ y along the two axes (x and y), from which parameters a, b, and c may be computed based on the following equations:

a = cos 2 θ 2σ 2 x + sin 2 θ 2σ 2 y , (7) b
= sin 2θ 4σ 2 y - sin 2θ 4σ 2 x , (8)
c = sin 2 θ 2σ 2 x + cos 2 θ 2σ 2 y . (9)
This choice has the advantage that all values for θ, σ x , and σ y , even negative ones, yield a positive-definite matrix C. By the way, this representation, based on rotation angles, is wellknown and used in evolution strategies [START_REF] Schwefel | Numerical optimization of computer models[END_REF]. 2) Initialization: The genomes of the agents in the initial population are generated at random, according to the following rules:

• the µ x , µ y position of the agent on Cipolla's diagram is randomly generated from a standard normal bivariate distribution centered at the origin and with variance 1 along both axes;

• σ x and σ y are independently extracted from a uniform distribution on [1, 2);

• θ is extracted from a uniform distribution on [0, 2π).

All individuals are assigned an initial wealth of 100.

To test the hypothesis that Cipolla's laws are sensitive to the initial distribution of the population, such distribution may be modified by filtering out generated agents that do not respect a given constraint. In particular, we are interested in performing simulations with an initial population comprising only

• stupid agents (µ x < 0 and µ y < 0);

• agents that we might call deleterious, bacause they cause on average lower gains or greater damages to others than the gains or losses they derive from their actions (µ y < -µ x ).

Although other initial distributions could be proposed, these are the two most "promising" initial distributions, besides the unfiltered one, to investigate in the light of Cipolla's observations. The unfiltered and the two filtered initial distribution are shown in Figure 2.

3) Fitness, Selection, and Reproduction: The fitness of an agent is given by its wealth.

If, at any time, an agent's wealth becomes negative, the agent dies. This is an implicit form of selection, which captures the fact that if an agent keeps on being harmed when interacting with its peers, eventually it will get extinct.

An agent whose wealth exceeds the predefined threshold of twice the initial wealth reproduces by asexual division. A new agent is created, which inherits (with mutation) the genome of its parent. The wealth is equally divided among the parent and offspring. This, too, may be regarded as an implicit form of selection, for the more an agent's behavior is capable of producing wealth, the more that agent will reproduce and pass on its genes.

As new agents may be generated by reproduction at any time as other agents may became extinct, in general the population size will fluctuate and in some cases may greatly increase or shrink until only one agents is left, depending on dynamics determined by the average behavior of agents in the population.

V. EXPERIMENTS AND RESULTS

We perform experiments by trying all possible combinations of the following parameters:

• i ∈ {all, stupid, deleterious} (initial distribution):
filter the initial distribution of the agents so that, respectively, all types of agents are generated, only stupid agents are generated, only deleterious agents (i.e., agents such that µ y < -µ x ): see Figure 2;

• f ∈ {linear, logarithmic, hyperbolic} (transfer function): the function used to modify wealth of agents participating in an interaction is, respectively, linear, logarithmic, or hyperbolic;

• d ∈ {false, true} (defense): rational agents (i.e., agents such that µ x > 0) know how to defend themselves from bandits;

• r ∈ {false, true} (relative): the effects of the interactions between agents are felt by either participant relative to their µ x ;

• z ∈ {false, true} (zero-sum): a zero-sum game is enforced, so that the total wealth of the population remains constant.

This gives a total number of 72 combinations, which will be encoded as a string of parameters as follows:

• Boolean parameters are omitted if their value is false, included otherwise;

• the remaining parameters will be followed by an abbreviation of their value;

• the parameters will be separated by hyphens.

For example, ia-flin-d-r-z will denote simulations with defense, unfiltered initial distribution, linear transfer function, relative effects, and zero-sum game enforcement.

All simulations are run with an initial population size of 1,000 agents and for one million time steps (i.e., agent interactions) or until the population size shrinks to one, in which case there is no point in continuing and we consider that the population has become extinct. In order to save computational resources, we set the maximum population size at 10,000. If that limit is reached, reproduction can happen only if another agent dies thus leaving one empty slot. However, this limitation does not change the results for the size of 10,000 is attained late in the simulation, if ever, when the final distribution type has already clearly emerged. to get a first idea about which settings are most compatible with Cipolla's laws.

is-flog is-flog-z is-flog-r is-flog-r-z is-flog-d is-flog-d-z is-flog-d-r is-flog-d-
-5 -4 -3 -2 -1 -4 -3 -2 -1 0 1 2 3 x y is-fhyp is-fhyp-z is-fhyp-r is-fhyp-r-z is-fhyp-d is-fhyp-d-z is-fhyp-d-r is-fhyp-d-r-z 5
y 0.8 1.0 1.2 1. 8 -6 -4 -2 0 -4 -2 0 2 4 x y -6 -5 -4 -3 -2 -1 -4 -2 0 2 4 x
There are some general observations we can make based on these distributions:

• it is clear that relativization of the effects of inter-actions (parameter r) is critical to the survival and proliferation of stupid agents; when r is turned off, the stupid and the helpless are mercilessly wiped out from the population, to the point that, even when they dominate the initial population, they manage to drive themselves into extinction;

• if we restrict our attention to runs with parameter r turned on, the most promising distributions may be observed when a zero-sum game is enforced (parameter z turned on): a clear preponderance of stupid agents shows up in simulations ia-flin-d-r-z, is-flog-r-z, is-flog-d-r-z, and id-flin-r-z;

• the only setting of simulation parameters which obtains a preponderance of stupid agents starting from a "neutral" initial distribution is ia-flin-d-r-z, with parameter d turned on: this suggests that the ability of rational agents to defend themselves from bandits may play a role in favoring a drift toward the evolutionary success of stupidity;

• an initial distribution biased toward stupid agents appears to have an impact on a prevalence of stupidity in the final population: a majority of the simulations that end up having a larger proportion of stupid agents than any other type has an initial population of stupid only (i = stupid).

Overall, eight parameter settings (six out of which assume an initial population of all stupid agents) achieved a final population distribution featuring a majority of stupid agents, though in some cases not an overwhelming one, as one would expect based on Cipolla's Laws. These are: ia-flin-d-r-z, is-flog-r-z, is-flog-d-r-z, is-fhyp-r, is-fhyp-r-z, is-fhyp-d-r, is-fhyp-d-r-z, and id-flin-r-z. We are now going to examine these simulations in deeper detail.

Figure 4 shows the evolution of the composition of the population during typical simulations with each of the eight selected parameter settings. This gives us an idea of the evolutionary dynamics at work during these simulation and how the observed final distribution was reached.

We may notice that in simulation starting from a stupidonly initial distribution, the takeover of the population by the stupid is accompanied by a relative increase of the helpless; this is not the case with the two settings ia-flin-d-r-z and id-flin-r-z, where the proliferation of the stupid is virulent and leaves little hope to the other types; nevertheless, even in these latter cases, it is the stupid with undertones of helplessness that are responsible for the takeover and a few of them seemingly continue to spill over into the H s sector, preventing it from becoming empty.

Figure 5 shows how the total wealth of the population changes during simulations for the two parameter settings where a zero-sum game is not enforced: in both cases, we observe the same pattern, whereby an initial crisis, whose cause is to be found in the great number of stupid agents with undertones of banditism (S b ), is overcome as their number decreases in favor of stupid agents with undertones of helplessness (S h ), some of which evolve into helpless with undertones of stupidity (H s ), as the population seemingly performs a 45 • clockwise rotation away from the B quadrant. The total wealth evolution for zero-sum game parameter settings is not shown because, by definition, total wealth remains constant in those simulations. Evolution of the total wealth of the population during the two simulations using parameter settings where a zero-sum game is not enforced.

Figure 6 shows the final wealth distribution for typical simulations with each of the eight selected parameter settings. An interesting remark is that simulations starting from an initial population of stupid agents only end up with a flat or right-skewed wealth distribution, which is highly unlikely and contrary to what is usually observed in real human populations, where the wealth distribution tends to follow a power law. In this sense, we might say that parameter settings ia-flin-d-r-z and id-flin-r-z look much more in line with what one would expect.

Figure 7 shows how the population size evolves in simulations using the eight selected parameter settings. It is interesting to observe that all simulations with i = stupid exhibit a similar pattern, whereby the population size starts to increase, culminates after ca. 200,000 periods (except in is-fhyp-r and is-fhyp-d-r, where the population size is initially steady), and suddenly plummets to a minimum around the 400,000th period; such crisis is then overcome, and the population size reaches what looks like an equilibrium in three out of the six settings, whereas it starts to increase again in the other cases. Surprisingly, the two remaining settings, namely ia-flin-d-r-z and id-flin-r-z, show almost identical dynamics despite their differences: the population size increases with two initial bumps followed by temporary stagnations, after which the increase becomes steady.

VI. CONCLUSION

Some of the parameter settings we have tried led to emergent behaviors that are quite in line with Carlo Cipolla's theory. One parameter setting in particular, namely ia-flin-d-r-z, looks like a very promising first approximation of Cipolla's laws. This is all the more remarkable as the initial population is unfiltered. One hypothesis which appears to be critical to the corroboration of the theory is that the interactions among the agents be a zero-sum game: this is by no means an obvious assumption and, therefore, it calls for an explanation. The one we can offer is that the subjective utility of the agents, i.e., their perception of their own welfare, is somehow relative to the welfare of their peers. In a sense, it would be as if an agent compared its welfare to that of their peers and considered itself happy to the extent that it were in a better situation and unhappy to the extent that it were in a worse situation. If this were the case, then obviously any increase of an agent's wealth would determine a corresponding reduction of the wealth of the rest of the population.

Our results should be regarded as a first step toward understanding which additional hypothesis must be made in order to reconcile Cipolla's theory with other well-corroborated theories, like the theory of evolution. Possible future work might include enriching the evolutionary algorithm with sexual reproduction (i.e., recombination) and testing other conjectures that may help reproducing a state of affairs compatible with Cipolla's laws, such as a different propensity of agents to reproduce depending their µ x and µ y .

Possible applications of our approach might be to model certain phenomena relevant to the behavior of human systems, like traffic jams, crowd movements, market bubbles and crashes, and the like. One might argue that stupid behaviors (in the sense of Cipolla) can actually be observed and probably play a significant role in those phenomena.

Fig. 1 .

 1 Fig. 1. Cipolla's diagram with its eight sectors. The uppercase letter identifies the quadrant and, therefore, the main character of an agent falling in the corrisponding sector. The lowercase subscript refers to the undertones of the agent, deriving from its closeness to the neighboring quadrant.

Fig. 2 .

 2 Fig. 2. The three initial distributions of the population used in the simulations: (a) the unfiltered standard normal bivariate distribution centered at the origin and with variance 1 along both axes; (s) the same distribution filtered to comprise stupid agents only; (d) the same distribution filtered to comprise deleterious agents only. The xand y-axis are drawn in red.

Figure 3

 3 Figure3gives a comprehensive picture of the final populadistribution for each parameter setting, which we can use

  id-flin-z id-flin-r id-flin-r-z id-flin-d id-flin-d-z id-flin-d-r id-flin-d-r-z

Fig. 3 .

 3 Fig. 3. The distribution of the population after 10 6 simulation steps or after the population size has shrinked to one (these cases are clearly visible, with a single agent in the center of the plot). The x = 0 and y = 0 lines are shown in red, to help the reader locate the individuals with respect to the four quadrant of Cipolla's diagram.

Fig. 4 .

 4 Fig. 4. Evolution of the composition of the population during a simulation with each of the eight selected parameter settings. The eight sub-populations located in the eight quadrants of Figure 1 are stacked in the diagrams starting from I b and proceeding counterclockwise up to B i and labeled on the lefthand and right-hand side of the diagram, if they make up for at least 5% of the population.

  Fig.5.Evolution of the total wealth of the population during the two simulations using parameter settings where a zero-sum game is not enforced.

Fig. 6 .

 6 Fig.6. Final wealth distribution at the end of simulations with the eight selected parameter settings.

Fig. 7 .

 7 Fig. 7. Evolution of the population size in simulations using the eight selected parameter settings.