FRANC ¸ OIS BOUCHUT Tom Ás

Morales De Luna

A SUBSONIC-WELL-BALANCED RECONST RUCTION SCHEME FOR SHALLOW WATER FLOWS

Keywords:

We consider the Saint-Venant system for shallow water flows with nonflat bottom. In past years, efficient well-balanced methods have been proposed in order to well resolve solutions close to steady states at rest. Here we describe a strategy based on a local subsonic steady state reconstruction that allows one to derive a subsonic-well-balanced scheme, preserving exactly all the subsonic steady states. It generalizes the now well-known hydrostatic solver, and like the latter it preserves the nonnegativity of the water height and satisfies a semidiscrete entropy inequality. An application to the Euler-Poisson system is proposed. Key words. shallow water, subsonic reconstruction, subsonic steady states

before (entropy satisfying and the conservativity and nonnegativity of ρ). The object of this paper is to go further in well-balanced schemes by building a solver with all the above requirements and, overall, the property of preserving exactly all subsonic steady states.

2. Saint-Venant system and well-balanced schemes. The Saint-Venant system describes the evolution of the water height ρ(t, x) and the velocity u(t, x) in the horizontal direction, of a thin layer of water flowing over a slowly varying topography. In one space dimension, the system writes as (2.1)

∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p(ρ)) + ρgz x = 0,
where g > 0 is the gravitational constant and z(x) is the topography. We shall denote (2.2)

Z = gz.
The physically relevant case is p(ρ) = gρ 2 /2, but we shall deal with the general case p(ρ). We shall assume as usual that p > 0, and we suppose that

(2.3) ρ 2 p (ρ) is strictly increasing, ρ 2 p (ρ) → ∞ as ρ → ∞, (2.4)
1 0 p (ρ) ρ dρ < ∞, p (ρ) → 0 as ρ → 0, ∞ 1 p (ρ) ρ dρ = ∞.
These assumptions are satisfied in particular for the pressure law of isentropic gas dynamics p(ρ) = κρ γ , with γ > 1 and κ > 0. We define as usual the internal energy e(ρ) by (2.5) e (ρ) = p(ρ) ρ 2 .

Note that the integrability conditions in (2.4) imply that e(ρ) + p(ρ)/ρ has a finite limit as ρ → 0 and tends to ∞ as ρ → ∞. For future reference we denote the flux by (2.6) F (U ) = ρu, ρu 2 + p(ρ) , U = (ρ, ρu).

The Saint-Venant model is very robust, being hyperbolic and admitting an entropy inequality related to the physical energy, (2.7)

∂ t η(U, Z) + ∂ x G(U, Z) ≤ 0, where (2.8) η(U ) = ρu 2 /2 + ρe(ρ), G(U ) = ρu 2 /2 + ρe(ρ) + p(ρ) u, η(U, Z) = η(U ) + ρZ, G(U, Z) = G(U ) + ρuZ.
The steady states of (2.1) can be described as follows. We subtract u times the first equation in (2.1) from the second and divide the result by ρ. We get (2.9)

∂ t u + ∂ x u 2 2 + e(ρ) + p(ρ) ρ + Z = 0.
Therefore, the steady states are exactly the functions ρ(x), u(x) satisfying (2.10) ρu = Cst, u 2 2 + e + p ρ (ρ) + Z = Cst.

In particular, we have the so-called steady state at rest (2.11) u = 0, e+ p ρ + Z = Cst.

As exposed in [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF], a first-order finite volume method for solving (2.1) writes generically with U = (ρ, ρu) as (2.12)

U n+1 i -U i + Δt Δx i F i+1/2--F i-1/2+ = 0, (2.13) F i+1/2-= F l (U i , U i+1 , ΔZ i+1/2 ), F i+1/2+ = F r (U i , U i+1 , ΔZ i+1/2 )
for some left/right numerical fluxes F l (U l , U r , ΔZ), F r (U l , U r , ΔZ), with

(2.14) ΔZ i+1/2 = Z i+1 -Z i ,
and where Δx i denotes a possibly variable mesh size, Δx i = x i+1/2x i-1/2 . We then have the following characterizations (see [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF]).

The conservativity or density writes, with F l/r = (F ρ l/r , F ρu l/r ),

(2.15)

F ρ l = F ρ r ≡ F ρ .
The consistency-conservativity can be written as (2.16) F ρ (U, U, 0) = ρu, F ρu l (U, U, 0) = F ρu r (U, U, 0) = ρu 2 + p(ρ),

(2.17) F ρu r (U l , U r , ΔZ) -F ρu l (U l , U r , ΔZ) = -ρΔZ + o(ΔZ) as U l , U r → U, ΔZ → 0.

The well-balancing property can be stated as the property of having, for the considered steady states,

(2.18) F l (U l , U r , ΔZ) = F (U l ), F r (U l , U r , ΔZ) = F (U r ).
The property of satisfying a semidiscrete entropy inequality (i.e., related to the limit Δt → 0) is characterized by the existence of a numerical entropy flux G(U l , U r , Z l , Z r ) consistent with the exact flux G(U, Z), such that

(2.19) G(U r , Z r ) + η (U r , Z r )(F r (U l , U r , ΔZ) -F (U r )) ≤ G(U l , U r , Z l , Z r ), (2.20) G(U l , U r , Z l , Z r ) ≤ G(U l , Z l ) + η (U l , Z l )(F l (U l , U r , ΔZ) -F (U l )),
where η (U, Z) is the derivative of η(U, Z) with respect to U .

The hydrostatic reconstruction scheme satisfies all the above and is defined as

(2.21) F l (U l , U r , ΔZ) = F(U * l , U * r ) + 0 p(ρ l ) -p(ρ * l ) , F r (U l , U r , ΔZ) = F(U * l , U * r ) + 0 p(ρ r ) -p(ρ * r )
, where F (U l , U r ) is a numerical flux for the shallow water problem without source (Z = cst), and the reconstructed states U * l , U * r are defined by (2.22)

U * l = (ρ * l , ρ * l u l ), U * r = (ρ * r , ρ * r u r ), (2.23) e + p/ρ (ρ * l ) = (e + p/ρ)(ρ l ) -(ΔZ) + + , e + p/ρ (ρ * r ) = (e + p/ρ)(ρ r ) -(-ΔZ) + + ,
where we use the notation X + = max(0, X), and we assume that e(ρ) + p(ρ)/ρ → 0 as ρ → 0.

3. Well-balanced scheme with subsonic reconstruction. We would now like to explain how it is possible to extend the hydrostatic reconstruction scheme (2.21)-(2.23) in order to obtain a scheme that satisfies the above requirements and preserves some more general steady states than the steady states at rest. We shall obtain in particular a scheme that preserves all subsonic steady states, that is, the steady states that verify u 2 < p (ρ). This property will be called subsonic-wellbalanced. Note in particular that the steady states at rest (with u = 0) are subsonic.

Parametrization of numerical fluxes.

Following [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF], [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF], we propose and analyze finite volume schemes defined by (2.12), (2.13) with numerical fluxes (3.1)

F l (U l , U r , ΔZ) = F (U * l , U * r ) + 0 p(ρ l ) -p(ρ * l ) + T l (U l , U r , ΔZ) , F r (U l , U r , ΔZ) = F(U * l , U * r ) + 0 p(ρ r ) -p(ρ * r ) + T r (U l , U r , ΔZ)
, where F stands for a numerical flux for the homogeneous problem (Z = cst), and the interface values U * l , U * r are derived from a local reconstruction procedure. They should satisfy at least that U * l = U l , U * r = U r when ΔZ = 0. The extra terms T l , T r appear here in order to balance the advection term ∂ x (ρu 2 ) in (2.1), which was not considered in the hydrostatic scheme. Taking into account the constant discharge condition in (2.10), this balancing requirement suggests the relations

(3.2) T l = ρ l u l (u l -u * l ), T r = ρ r u r (u r -u * r )
. It is obvious from the characterization (2.18) that a steady state is maintained exactly by the scheme (2.12), (2.13), (3.1) if, for such a state, the reconstructed states satisfy

U * l = U * r , ρ * l u * l = ρ l u l ,
⎧ ⎪ ⎨ ⎪ ⎩ (u * l ) 2 2 + e + p ρ (ρ * l ) + Z * = u 2 l 2 + e + p ρ (ρ l ) + Z l , ρ * l u * l = ρ l u l , (3.4)
⎧ ⎪ ⎨ ⎪ ⎩ (u * r ) 2 2 + e + p ρ (ρ * r ) + Z * = u 2 r 2 + e + p ρ (ρ r ) + Z r , ρ * r u * r = ρ r u r , with (3.5) Z * = max(Z l , Z r ).
There exist solutions to the previous system if ΔZ is small enough (recall that ΔZ = Z r -Z l ), but this is not true for arbitrary ΔZ, as we shall see later. The idea is thus to consider generalized T l , T r , to be defined later. Their definition is motivated by the entropy inequality. Lemma 3.1. Let F(U l , U r ) be a given consistent numerical flux for the Saint-Venant problem without source that verifies a semidiscrete entropy inequality for the entropy pair (η, G) given by (2.8), and denote F = (F ρ , F ρu ). A sufficient condition for the scheme (2.12), (2.13), (3.1) to be semidiscrete entropy satisfying for the entropy pair

( η, G) in (2.8) is that for some Z * (3.6) G(U * l ) + η (U * l )(F (U * l , U * r ) -F (U * l )) + F ρ (U * l , U * r )Z * ≤ G(U l ) + η (U l )(F l (U l , U r , ΔZ) -F (U l )) + F ρ (U * l , U * r )Z l and (3.7) G(U r ) + η (U r )(F r (U l , U r , ΔZ) -F (U r )) + F ρ (U * l , U * r )Z r ≤ G(U * r ) + η (U * r )(F (U * l , U * r ) -F (U * r )) + F ρ (U * l , U * r )Z * .
Proof. The numerical flux F satisfies a semidiscrete entropy inequality associated with the entropy pair (η, G); thus

(3.8) G(U r ) + η (U r )(F (U l , U r ) -F (U r )) ≤ G(U l , U r ) ≤ G(U l ) + η (U l )(F (U l , U r ) -F (U l ))
for a numerical flux G consistent with G. For the scheme (2.12), (2.13), (3.1) to be semidiscrete entropy satisfying for the entropy pair ( η, G), (2.19), (2.20)

should hold. Let (3.9) G(U l , U r , Z l , Z r ) = G(U * l , U * r ) + F ρ (U * l , U * r )Z * . As G is consistent with G, G is consistent with G.
W l ≡ F ρ • e + p ρ (ρ l ) -e + p ρ (ρ * l ) + Z l -Z * + (u * l ) 2 2 - u 2 l 2 + (u l -u * l ) (F ρu -p(ρ * l )) + u l T l , (3.10) W r ≡ F ρ • e + p ρ (ρ r ) -e + p ρ (ρ * r ) + Z r -Z * + (u * r ) 2 2 - u 2 r 2 + (u r -u * r ) (F ρu -p(ρ * r )) + u r T r . (3.11)
A necessary and sufficient condition for (3.6), (3.7) to hold is that

(3.12) W l ≥ 0, W r ≤ 0.
Proof. From the explicit value of F , G and computing η (U ) = (e(ρ) + p(ρ)/ρu 2 /2, u), one gets the identity G(U )η (U )F (U ) = -u p(ρ). Plugging this into (3.6), (3.7) yields the result.

We remark that in the particular case when U * l , U * r satisfy (3.3)-(3.4), we have (3.13)

W l = (u * l -u l ) (u l + u * l )F ρ -F ρu + p(ρ * l ) + u l T l , W r = (u * r -u r ) (u r + u * r )F ρ -F ρu + p(ρ * r ) + u r T r .
Thus, the choice of T l , T r given by (3.2) makes W l = W r = 0 whenever U * l = U * r .

Subsonic reconstruction.

We intend here to define reconstructed states that satisfy (3.3), (3.4). Nevertheless, it will not be always possible to do so, as we will see in what follows. The possibility of achieving these relations is related to the following definitions.

Definition 3.3. Let ρ ≥ 0, and let u ∈ R. We say that (ρ, u) is a sonic, subsonic, or supersonic point for the Saint-Venant system (2.1) if we have, respectively,

u 2 = p (ρ), u 2 < p (ρ), or u 2 > p (ρ).
Definition 3.4. Let q ∈ R. We define ρ s (q) as the solution to

(3.14) ρ 2 s p (ρ s ) = q 2 , ρ s ≥ 0, and
(3.15) m s (q) = e + p ρ + p 2 (ρ s (q)).
Because of the assumptions (2.3), (2.4) on p, there exists a unique ρ s solution to (3.14). Moreover, m s (0) is well defined:

(3.16) m s (0) = e + p ρ (0).
In the particular case when p(ρ) = κρ γ , we have

(3.17) ρ s (q) = q 2 κγ 1 γ+1 , m s (q) = 1 2 + 1 γ -1 (κγ) 2 γ+1 |q| 2 γ-1 γ+1 .
The following proposition gives the necessary and sufficient conditions for the existence of a solution to (3.3) and (3.4). Proposition 3.5.

Let u 0 ∈ R, ρ 0 ≥ 0, δ ∈ R. Consider the system ⎧ ⎪ ⎨ ⎪ ⎩ (u * ) 2 2 + e + p ρ (ρ * ) = u 2 0 2 + e + p ρ (ρ 0 ) + δ, ρ * u * = ρ 0 u 0 , ρ * ≥ 0, u * ∈ R, (3.18)
and denote ρ s ≡ ρ s (ρ 0 u 0 ). There exists a solution (ρ * , u * ) to (3.18) if and only if

(3.19) u 2 0 2 + e + p ρ (ρ 0 ) + δ ≥ m s (ρ 0 u 0 ). ρ * sup ρ s ρ * sub ρ 0 |δ| f (q 0 , ρ) f (q 0 , ρ 0 ) f (q 0 , ρ 0 ) + δ m s Fig. 1. Function f (q 0 , •).
Moreover, the following hold: (i) If we have equality in (3.19), there is only one solution (ρ * , u * ) to (3.18), and it is given by

(3.20) ρ * = ρ s , u * = ρ 0 u 0 ρ s if ρ 0 u 0 = 0, 0 if ρ 0 u 0 = 0.
(ii) If we have a strict inequality in (3.19), then there are exactly two different solutions

(ρ * sup , u * sup ) and (ρ * sub , u * sub ) to (3.18), with ρ * sup ≤ ρ s < ρ * sub , and ρ * sup < ρ s for ρ 0 u 0 = 0.
(iii) A solution (ρ * , u * ) to (3.18), with ρ * u * = 0, is a sonic (resp., subsonic or supersonic) point if and only if ρ * = ρ s (resp., ρ * > ρ s or ρ * < ρ s ).

Proof. Let us suppose first that ρ 0 u 0 = 0, and consider the function

(3.21) f : R × (0, ∞) -→ R, (q, ρ) → f (q, ρ) = q 2 2ρ 2 + e + p ρ (ρ).
Then (ρ * , u * ) is a solution to (3.18) if and only if ρ * u * = ρ 0 u 0 and f (ρ 0 u 0 , ρ * ) = f (ρ 0 u 0 , ρ 0 ) + δ. We have ∂f ∂ρ (q, ρ) = (ρ 2 p (ρ)q 2 )/ρ 3 ; thus, according to (2.3), (2.4), f is strictly increasing in (ρ s (q), ∞) and strictly decreasing in (0, ρ s (q)). Therefore, ρ s (q) is a minimum point of f with minimum value f (q, ρ s (q)) = m s (q). Figure 1 shows a sketch of the function f (ρ 0 u 0 , •) and the solutions ρ * sub , ρ * sup in the case δ < 0. Thus, condition (3.19) follows, as well as (i) and (ii). Now, if we consider (ρ * , u * ) a solution of (3.18), we have

(3.22) (u * ) 2 > p (ρ * ) ⇔ (ρ 0 u 0 ) 2 > (ρ * ) 2 p (ρ * ). Since ρ 2 p (ρ) is a strictly increasing function with ρ 2 s p (ρ s ) = (ρ 0 u 0 ) 2 , (3.23) (u * ) 2 > p (ρ * ) ⇔ ρ * < ρ s ,
which proves (iii).

In the case ρ 0 u 0 = 0, the second condition in (3.18) simplifies to ρ * u * = 0. Thus the system can have solutions with either ρ * = 0, giving (u * ) 2 /2 = r.h.s.m s (0), or u * = 0, giving (e + p/ρ)(ρ * ) = r.h.s. One easily sees that a solution exists if and only if r.h.s. ≥ m s (0), proving condition (3.19). In case of equality, the only solution is ρ * = 0, u * = 0, which is sonic. In case of inequality, the solutions are given by

(3.24) ρ * sup = 0, (u * sup ) 2 2 = u 2 0 2 + e + p ρ (ρ 0 ) + δ -m s (0), (3.25) u * sub = 0, e + p ρ (ρ * sub ) = u 2 0 2 + e + p ρ (ρ 0 ) + δ.
This gives the result, with the convention that for (3.24) we identify the two solutions having density ρ * sup = 0. Note that these solutions are supersonic, while the ones corresponding to (3.25) are subsonic.

Corollary 3.6. Let ρ 0 > 0, and let u 0 ∈ R\{0}. Then (ρ 0 , u 0 ) is a sonic (resp., subsonic or supersonic) point if and only if

ρ 0 = ρ s (ρ 0 u 0 ) (resp., ρ 0 > ρ s (ρ 0 u 0 ) or ρ 0 < ρ s (ρ 0 u 0 )).
Proof. The point (ρ 0 , u 0 ) is a trivial solution to (3.18) with δ = 0. Thus the result follows from (iii) in the previous proposition.

Lemma 3.7. Suppose that in Proposition 3.5 we are in the case of two solutions ρ * sub , ρ * sup . Then one has the following ordering: (A) Case (ρ 0 , u 0 ) subsonic:

(A.1) if δ > 0, then ρ * sup < ρ 0 < ρ * sub ; (A.2) if δ ≤ 0, then ρ * sup < ρ * sub ≤ ρ 0 . (B) Case (ρ 0 , u 0 ) supersonic: (B.1) if δ ≥ 0, then ρ * sup ≤ ρ 0 < ρ * sub ; (B.2) if δ < 0, then ρ 0 ≤ ρ * sup < ρ * sub . (C) Case (ρ 0 , u 0 ) sonic: ρ * sup ≤ ρ 0 < ρ * sub .
The proof is left to the reader. Now, in order to define a scheme that preserves the nonnegativity of the density, the reconstructed states need to verify ρ * l ≤ ρ l and ρ * r ≤ ρ r (see Theorem 3.12(i)). We notice that (3.3), (3.4) correspond to the problem (3.18) with, successively, δ = Z l -Z * , δ = Z r -Z * . Since it is natural to try to choose a solution (ρ * , u * ) of the same sonicity as (ρ 0 , u 0 ), one sees with the previous lemma that in order to have ρ * ≤ ρ 0 one needs δ ≤ 0 in the subsonic case, and δ ≥ 0 in the supersonic case. Recalling the values δ = Z l -Z * and δ = Z r -Z * , this gives that one should have Z * ≥ max(Z l , Z r ) in subsonic regions, while Z * ≤ min(Z l , Z r ) in supersonic regions. We then observe that it is not possible to satisfy both conditions with Z * depending continuously on the data. Thus we make the choice of solving exactly the subsonic steady states and disregard supersonic steady states, which justifies Z * = max(Z l , Z r ), i.e., (3.5).

Consider now the function f given by (3.21). For q ∈ R fixed, f (q, •) is strictly increasing in [ρ s (q), ∞) (recall that ρ s (q) and m s (q) are defined as (3.14), (3.15)):

(3.26) f (q, •)| [ρ s (q),∞) : [ρ s (q), ∞) → [m s (q), ∞), ρ → q 2 2ρ 2 + e + p ρ (ρ).
We consider its inverse and denote it by f -1 r (q, •):

(3.27) f -1 r (q, •) : [m s (q), ∞) → [ρ s (q), ∞).
This inverse function corresponds to choosing the subsonic solution to (3.18). With the choice (3.5), equations (3.3), (3.4) correspond to the problem (3.18) with, successively,

δ = Z l -Z * = -(ΔZ) + , δ = Z r -Z * = -(-ΔZ) + .
Therefore, we define the reconstructed states U * l , U * r by (3.28)

ρ * l = min ρ l , f -1 r ρ l u l , max f (ρ l u l , ρ l ) -(ΔZ) + , m s (ρ l u l ) , u * l = ρ l u l /ρ * l (u * l = u l if ρ * l = 0), ρ * r = min ρ r , f -1 r ρ r u r , max f (ρ r u r , ρ r ) -(-ΔZ) + , m s (ρ r u r ) , u * r = ρ r u r /ρ * r (u * r = u r if ρ * r = 0).
Note that these definitions imply that

(3.29) ρ * l u * l = ρ l u l , ρ * r u * r = ρ r u r .
Lemma 3.8. The definitions (3.28) can be interpreted as follows:

(A) Case ΔZ ≤ 0:

We have the trivial solution to

(3.3) U * l = U l . (B)
The general case: By Proposition 3.5, we know that the system (3.3) has a solution if and only if

(3.30) u 2 l 2 + e + p ρ (ρ l ) -(ΔZ) + ≥ m s (ρ l u l ). (B.1) If (ρ l , u l ) is a supersonic point or a sonic point, then U * l = U l . (B.2) If (ρ l , u l
) is a subsonic point and we have strict inequality in (3.30), then (ρ * l , u * l ) is the subsonic solution to (3.3). (B.3) If (ρ l , u l ) is a subsonic point and we have equality in (3.30) or the inequality is not satisfied, then ρ * l = ρ s (ρ l u l ). Similar statements hold for (ρ * r , u * r ). The case ρ * l = 0 could pose some problems in the previous definition of u * l , but as we consider conservative variables, the product ρ * l u * l is well defined. The following result shows that there is indeed no problem of continuity.

Lemma 3.9. The reconstructed states (3.28) verify the following:

(i) (3.31) min ρ l , ρ s (ρ l u l ) ≤ ρ * l ≤ ρ l , min ρ r , ρ s (ρ r u r ) ≤ ρ * r ≤ ρ r ;
(ii) independently of the other arguments, one has

(3.32) lim ρ l →0 ρ * l = 0, lim ρ r →0
ρ * r = 0;

(iii) for ρ l ≥ 0, ρ r ≥ 0, ΔZ fixed, we have

(3.33) lim u l →0 ρ * l = e + p ρ -1 max e + p ρ (ρ l ) -(ΔZ) + , m s (0) , (3.34) lim u r →0 ρ * r = e + p ρ -1 max e + p ρ (ρ r ) -(-ΔZ) + , m s (0) ;
(iv) for ρ l , ρ r bounded, one has

(3.35) lim u l →0 u * l = 0, lim u r →0 u * r = 0.
Proof. Note that (iii) means the continuity of ρ * l and ρ * r in this asymptotics. We shall give only the proof for ρ * l , the proof for ρ * r being similar. According to Lemma 3.8, the only case when (3.31) is nontrivial is (B.2) with ΔZ ≥ 0. Then Corollary 3.6 allows one to conclude the proof of (i). Then, (ii) is a consequence of (i).

In order to prove (iii), for ρ l ≥ 0 and ΔZ fixed, we shall denote

(3.36) β l (u l ) = max f (ρ l u l , ρ l ) -(ΔZ) + , m s (ρ l u l ) , α l (u l ) = f -1 r (ρ l u l , β l (u l )).
We have (3.37) lim

u l →0 β l (u l ) = max e + p ρ (ρ l ) -(ΔZ) + , m s (0) . Since β l (u l ) ≥ m s (ρ l u l ), we have α l (u l ) ≥ ρ s (ρ l u l ), which yields (3.38) 0 ≤ (ρ l u l ) 2 α l (u l ) 2 ≤ (ρ l u l ) 2 ρ s (ρ l u l ) 2 = p ρ s (ρ l u l ) .
Therefore, (ρ l u l ) 2 /α l (u l ) 2 tends to 0 as u l → 0. Then, using the identity Since the right-hand side is at most ρ l , we deduce that ρ * l = min{ρ l , α l (u l )} has the same limit, which concludes the proof of (iii).

(3.39) (ρ l u l ) 2 2(α l (u l )) 2 + e + p ρ α l (u l ) = β l (u l ),
The last statement (iv) is also a consequence of (i) since either ρ * l ≥ ρ l , giving

|u * l | ≤ |u l |, or ρ * l ≥ ρ s (ρ l u l ), giving (3.41) (u * l ) 2 = (ρ l u l ) 2 (ρ * l ) 2 ≤ (ρ l u l ) 2 ρ s (ρ l u l ) 2 = p ρ s (ρ l u l ) .
Since the right-hand side tends to 0, this concludes the proof. We would like to end this subsection by giving an iterative procedure in order to solve the system (3.18). According to Lemma 3.8, we need to solve this system only in the case (B.2) with ΔZ > 0. This is done as follows.

Proposition 3.10. Let u 0 ∈ R\{0}, ρ 0 > 0, and δ < 0. We suppose that u 2 0 < p (ρ 0 ) and that (3.19) is strictly satisfied. Let V 0 = u 2 0 2 + (e + p ρ )(ρ 0 ) + δ, and let ψ(ρ) = ρ α (f (q 0 , ρ) -V 0 ), where f is the function given by (3.21), and q 0 = ρ 0 u 0 . Then, for α ≥ 3/2, the relation

(3.42) ρ n+1 = ρ n - ψ(ρ n ) ψ (ρ n )
(starting from ρ 0 ) defines a decreasing sequence that converges to ρ * sub , the subsonic solution to (3.18).

Proof. According to Proposition 3.5, there is a unique solution ρ * sub ∈ (ρ s , ρ 0 ) to the equation ψ(ρ) = 0. We have (3.43)

ψ (ρ) = αρ α-1 f (q 0 , ρ) -V 0 + ρ α ∂f ∂ρ (q 0 , ρ), ψ (ρ) = α(α -1)ρ α-2 f (q 0 , ρ) -V 0 + 2αρ α-1 ∂f ∂ρ (q 0 , ρ) + ρ α ∂ 2 f ∂ρ 2 (q 0 , ρ),
and, for ρ ≥ ρ * sub , (3.44)

f (q 0 , ρ) -V 0 ≥ 0, ∂f ∂ρ (q 0 , ρ) = ρ 2 p (ρ) -q 2 0 ρ 3 > 0, ∂ 2 f ∂ρ 2 (q 0 , ρ) = - 3 ρ ∂f ∂ρ (q 0 , ρ) + 1 ρ 3 ∂ ∂ρ ρ 2 p (ρ) ≥ - 3 ρ ∂f ∂ρ (q 0 , ρ).
Thus, for α ≥ 3/2, ψ is strictly increasing and convex. Therefore, the Newton method converges to the zero of the function, and this proves the result.

Definition of left and right fluxes.

In order to define the left and right numerical fluxes as (3.1), taking the definitions (3.28) for U * l , U * r (recall that Z * is given by (3.5)), we still need to define T l and T r . We denote (F ρ , F ρu ) ≡ F(U * l , U * r ) and define (3.45)

T l (U l , U r , ΔZ) = ρ l -ρ * l ρ * l F ρu -p(ρ * l ) -u * l F ρ -(u * l -u l ) F ρ + e + p ρ (ρ * l ) -e + p ρ (ρ l ) + (ΔZ) + + (u * l ) 2 2 - u 2 l 2 F ρ u l , (3.46) T r (U l , U r , ΔZ) = ρ r -ρ * r ρ * r F ρu -p(ρ * r ) -u * r F ρ -(u * r -u r ) F ρ + e + p ρ (ρ * r ) -e + p ρ (ρ r ) + (-ΔZ) + + (u * r ) 2 2 - u 2 r 2 F ρ u r .
The definition (3.45) (resp., (3.46)) is ambiguous in the case when u l = 0 or ρ * l = 0 (resp., when u r = 0 or ρ * r = 0). In order to overcome this difficulty, here we make the convention that "0/0 = 0." Then, one has to take into account the following remarks. They are stated for T l , but of course similar statements hold for T r .

1. If (ρ * l , u * l ) solves the system (3.3) (with (3.5)), the factor of F ρ u l in T l vanishes. 2. If (ρ l , u l ) is a supersonic point, we have U * l = U l and T l ≡ F ρ u l (ΔZ) + , which is well defined since u 2 l > p (ρ l ) ≥ 0. 3. If (ρ l , u l ) is a subsonic point and

u 2 l 2 +(e+ p ρ )(ρ l )-(ΔZ) + > m s (ρ l u l )
, according to (B.2) in Lemma 3.8 we have that ρ * l > 0, and the factor of F ρ u l in T l vanishes; thus T l is well defined. This is true in particular when (ρ l , u l ) is a subsonic point if we consider a continuous bottom z(x), which implies that ΔZ is small for a sufficiently fine grid.

4. If (ρ l , u l ) is a subsonic point and

u 2 l 2 + (e + p ρ )(ρ l ) -(ΔZ) + ≤ m s (ρ l u l
), but u l = 0, according to (B.3) in Lemma 3.8 we have ρ * l = ρ s (ρ l u l ) > 0; thus T l is well defined.

5. Some difficulties may arise for (ρ l , u l ) sonic close to (0, 0), or for (ρ l , u l ) subsonic with

u 2 l 2 + (e + p ρ )(ρ l ) -(ΔZ) + ≤ m s (ρ l u l
) and u l close to 0 (which also implies that (ρ l , u l ) is close to (0, 0)). In these cases, the flux F has to verify some conditions in order to define

F ρ u l and (ρ * l ) -1 F ρu -p(ρ * l ) -u * l F ρ .
As F is consistent with F , we expect these quantities to be unambiguously defined.

Lemma 3.11. The definitions (3.45), (3.46) of T l and T r imply that the conditions (3.12) are satisfied (with Z * given by (3.5)).

Proof. Consider the case of W l . We have

(3.47) u l ρ l -ρ * l ρ * l = u * l -u l ; thus (3.48) u l T l = u * l -u l F ρu -p(ρ * l ) -u * l F ρ -u l F ρ + e + p ρ (ρ * l ) -e + p ρ (ρ l ) + (ΔZ) + + (u * l ) 2 2 - u 2 l 2 F ρ = u * l -u l F ρu -p(ρ * l ) + e + p ρ (ρ * l ) -e + p ρ (ρ l ) + (ΔZ) + - (u * l ) 2 2 + u 2 l 2 F ρ .
Putting this value in (3.10) gives W l = 0.

3.4.

Properties of the subsonic reconstruction scheme. Theorem 3.12. Let F (U l , U r ) be a given consistent numerical flux for the Saint-Venant problem without source that preserves the nonnegativity of ρ by interface and satisfies a semidiscrete entropy inequality for the entropy pair (η, G) given by (2.8). Then the scheme (2.12), (2.13), with numerical fluxes (3.1), (3.28), (3.45), (3.46), (0) is conservative in density, (i) preserves the nonnegativity of ρ by interface, (ii) preserves the discrete subsonic steady states, (iii) is consistent with the Saint-Venant system away from sonic points, and (iv) satisfies a semidiscrete entropy inequality associated with the entropy pair ( η, G) in (2.8).

Proof. Notice first that when ΔZ = 0 we have U * l = U l , U * r = U r , T l = 0, and T r = 0, so that the scheme (3.1) reduces to the conservative scheme with numerical flux F.

Property (0) is obvious from the definition (3.1) and the characterization (2.15). For (i), the assumption that F (U l , U r ) preserves nonnegativity of ρ by interface means that (see [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF]) there exists some

σ l (U l , U r ) < 0 < σ r (U l , U r ) such that (3.49) ρ l + F ρ (U l , U r ) -ρ l u l σ l (U l , U r ) ≥ 0, ρ r + F ρ (U l , U r ) -ρ r u r σ r (U l , U r ) ≥ 0
(a) Case (ρ, u) supersonic. Then (ρ l , u l ) is also supersonic if it is close enough to (ρ, u) (and in particular u l = 0), and we have

(3.58) ρ * l = ρ l , u * l = u l , ρ * r = ρ r , u * r = u r , (3.59) F ρu r -F ρu l = -T l = -F ρ (U l , U r ) ΔZ u l = -ρΔZ + o(ΔZ).
(b) Case (ρ, u) subsonic. Then ρ > 0, and u 2 2 + (e + p ρ )(ρ) > m s (ρu). Therefore, for U l close enough to U and ΔZ small enough, we have (ρ l , u l ) subsonic, ρ l > 0,

u 2 l 2 + (e + p ρ )(ρ l ) -ΔZ > m s (ρ l u l ). From (B.
2) in Lemma 3.8 we have that ρ * l > 0, and we compute (3.60)

ρ * l = f -1 r ρ l u l , u 2 l 2 + e + p ρ (ρ l ) -ΔZ = f -1 r ρ l u l , f (ρ l u l , ρ l ) -ΔZ = ρ l + O(ΔZ), (3.61) u * l = ρ l u l ρ * l = u l + O(ΔZ).
Now, let us denote φ(ρ) = (e + p ρ )(ρ). Since φ is a strictly increasing function, we can consider its inverse φ -1 . According to (3.60),

(3.62) e + p ρ (ρ * l ) = e + p ρ (ρ l ) + u 2 l 2 - (u * l ) 2 2 -ΔZ;
thus using (3.61) we get (3.63)

ρ * l = φ -1 φ(ρ l ) + u 2 l 2 - (u * l ) 2 2 -ΔZ = ρ l + (φ -1 ) (φ(ρ l )) u 2 l 2 - (u * l ) 2 2 -ΔZ + O u 2 l 2 - (u * l ) 2 2 -ΔZ 2 = ρ l + ρ l p (ρ l ) u 2 l 2 - (u * l ) 2 2 -ΔZ + O(ΔZ) 2 .
Then, (3.64)

p(ρ * l ) = p(ρ l ) + p (ρ l )(ρ * l -ρ l ) + O(ΔZ) 2 = p(ρ l ) -ρ l ΔZ + ρ l u 2 l 2 - (u * l ) 2 2 + O(ΔZ) 2 .
We also have

(3.65) ρ l -ρ * l ρ * l F ρu -p(ρ * l ) -u * l F ρ = O(ΔZ) • o(1) = o(ΔZ).

and

(3.72)

F ρu r -F ρu l = p(ρ * l ) -p(ρ l ) -T l = p(ρ * l ) -p(ρ l ) + (u * l -u l )F ρ + O(ΔZ) + O |ΔZ| 1/2 |U l -U r | = (u l -u * l ) ρ l u l + u * l 2 -F ρ + O(ΔZ) + O |ΔZ| 1/2 |U l -U r | = O(ΔZ) + O |ΔZ| 1/2 |U l -U r | .
Thus, in any case, (3.67) is valid. This property (3.67) does not mean consistency in the sense of (2.17), but anyway, one can make the following observation. Assume that, at the point considered, one has dZ/dx = 0. Then ΔZ = o(Δx), U r -U l = O(Δx), and therefore (3.67) yields

(3.73) F ρu r -F ρu l = o(Δx),
which means consistency with the vanishing source. Since the condition dZ/dx = 0 is generically satisfied at sonic points (see [START_REF] Chinnayya | A well-balanced numerical scheme for the approximation of the shallow-water equations with topography: The resonance phenomenon[END_REF]), it justifies the global consistency of the scheme, except maybe close to the point (ρ, u) = (0, 0). We shall see in the numerical computations that even if the numerical fluxes can sometimes take large values, the scheme behaves reasonably well in the presence of data close to (ρ, u) = (0, 0).

Extension to the two-dimensional case.

The proposed scheme can be easily adapted to the two-dimensional (2D) case by using the usual approach of a finite volume method for multidimensional systems. We describe here the main ideas and refer the reader to [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF] and the references therein for further details.

Consider the system (4.1)

∂ t U + ∂ 1 F 1 (U, Z) + ∂ 2 F 2 (U, Z) + B 1 (U, Z)∂ 1 Z + B 2 (U, Z)∂ 2 Z = 0.
We discretize the space domain by means of a mesh made of cells C i , the control volumes. The interface Γ ij between two cells C i and C j is assumed to be a segment, and we denote by n ij the unit normal vector oriented from C i to C j . Then, define the numerical scheme (4.2)

U n+1 i = U n i - Δt |C i | j∈K i |Γ ij |F ij = 0,
where

U i represents the cell-average of U over C i , |C i | is the volume of C i , |Γ ij | is the length of Γ ij , and
F ij = F (U i , U j , Z i , Z j , n ij ) is a numerical flux approximating n 1 F 1 + n 2 F 2 which
gives the exchange term between C j and C i . Now, consider the particular case of 2D shallow water equations

(4.3) ⎧ ⎨ ⎩ ∂ t ρ + ∂ 1 (ρu) + ∂ 2 (ρv) = 0, ∂ t (ρu) + ∂ 1 (ρu 2 + p(ρ)) + ∂ 2 (ρuv) + ρg∂ 1 z = 0, ∂ t (ρv) + ∂ 1 (ρuv) + ∂ 2 (ρv 2 + p(ρ)) + ρg∂ 2 z = 0,
where ρ(t, x) ≥ 0 is the water height, u(t, x), v(t, x) are the two components of the velocity field, p(ρ) is the pressure, and z(x) is the topography. This system has the property of being invariant by rotation, which enables us to easily define a numerical flux for the 2D system from a given one-dimensional (1D) flux (see [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF]).

Thus, the fluxes at each interface can be computed with the numerical fluxes defined in the previous section after an appropriate rotation.

The scheme will then automatically preserve the nonnegativity of the water height and satisfy a semidiscrete entropy inequality.

But concerning well-balancing, it follows from (4.1) that steady states that can be preserved by such a method are those that can be seen as interface by interface steady states. In particular, the 2D scheme will be able to preserve any lake-at-rest steady state but not all subsonic equilibria. Only those equilibria that correspond to the rotation of a 1D subsonic steady state shall be preserved (provided that we choose a rectangular mesh in the direction of the fluid). But any more general steady state where the dependence in x 1 and x 2 is nontrivially balanced cannot be preserved.

Application to the

⎧ ⎨ ⎩ ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p(ρ)) = -ρ∂ x φ, -∂ 2 xx φ = ρ -ρ b ,
where ρ b ≡ ρ b (x) ≥ 0 is given. The system (5.1) is set for t > 0, 0 < x < l, with initial and boundary conditions (5.2)

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ρ(t = 0, •) = ρ b , u(t = 0, •) = u 0 , φ(l) -φ(0) = V, ρu(t, x = 0) = q 0 ≥ 0, ρ(t, x = l) = ρ b (l).
As usual, one has to complete (5.1) by an entropy inequality. In order to describe the steady states, we subtract u times the first equation in (5.1) from the second, divide the result by ρ, and get

(5.3) ∂ t u + ∂ x (u 2 /2 + e(ρ) + p(ρ)/ρ + φ) = 0.
Therefore, the steady states are determined by the relations (5.4) ρu = Cst,

u 2 2 + e(ρ) + p(ρ) ρ + φ = Cst.
We can observe that the Euler-Poisson system is of the type (2.1), where the bottom Z is replaced by a function φ that is time dependent. The subsonic reconstruction scheme can be applied to this system by "freezing" the potential on a time interval as follows. Given an approximation of φ at time t n , φ n = φ(t n , •), we solve the system (5.5)

∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p(ρ)) = -ρ∂ x φ n
in the interval [t n , t n+1 ) using the subsonic reconstruction scheme, where φ n stands for Z. We obtain approximations ρ n+1 , u n+1 . Finally, we solve the ODE in order to get the new potential. It is obvious that this algorithm is well-balanced since the freezing of the potential does not introduce any error in the case of a steady state.

(5.6) -∂ 2 xx φ n+1 = ρ n+1 -ρ b
The interest of the subsonic reconstruction scheme in this context is the ability to compute with high accuracy flows which are close to a subsonic steady state with constant discharge ρu = 0.

Numerical results.

6.1. Saint-Venant system. In order to evaluate our method, we compare the subsonic reconstruction scheme described here to the original hydrostatic reconstruction scheme of [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. We use first-order resolution, with the CFL 1 condition induced by the nonnegativity of density (see the proof of (i) in Theorem 3.12). A second-order extension can be used as in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF], and in that case no major differences are observed between the two reconstructions. The numerical flux F chosen here is the one obtained from the Suliciu relaxation system described in [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF]. We take p(ρ) = gρ 2 /2, g = 9.81 and use 200 points in the considered interval in each case.

We consider first a subsonic steady state in the interval (0, 100). The initial data are given by (6.1)

ρ 0 (x) = 3 i fx ≤ 50, 14.2175 if x > 50, u 0 (x) = 5 ifx ≤ 50, 1.055 if x > 50, (6.2) z(x) = 10 if x ≤ 50, 0 if x > 50.
The results are shown in Figures 2 and3. As we see, the subsonic reconstruction scheme maintains the subsonic steady state, as we have proved. This is an improvement with respect to the hydrostatic reconstruction scheme, which does not. Then, we consider a supersonic steady state, still in the interval (0, 100):

(6.3) ρ 0 (x) = 2 if x ≤ 50, 0.635 if x > 50, u 0 (x) = 5 i fx ≤ 50, 15.7474 if x > 50, (6.4) z(x) = 10 if x ≤ 50, 0 if x > 50.
We see in Figures 4,5, 6, and 7 that neither of the two schemes gives the right solution, but the subsonic reconstruction scheme is more accurate.

We consider now a classical transcritical shock test. The space domain is (0, 25), the initial data are ρ 0 (x) = 0.33, u 0 (x) = 0.18/0.33, and the topography is (6.5) z(x) = 0.2 -0.05(x -10) 2 if 8 < x < 12, 0 o t h e r w i s e .

The boundary conditions are taken as ρu(x = 0) = 0.18 and ρ(x = 25) = 0.33. The results are shown in Figures 8 and9. We see that the subsonic reconstruction scheme gives a solution which is sharper on the left part of the discontinuity than the hydrostatic reconstruction. In Figure 10 we show the result with the boundary discharge q 0 = 10. We use 100 points in space, and the final time is t = 100. A second-order reconstruction is used. As we see, we have reached a subsonic equilibrium where u 2 /2 + e + p/ρ + φ is constant and q is almost constant.

Finally the same test with boundary discharge q 0 = 40 is shown in Figure 11. Even if we have reached a steady state, q and u 2 /2 + e + p/ρ + φ are not constant. We observe a jump at the point where there is a change from a supersonic regime to a subsonic one. This is not surprising since the scheme is not exact for supersonic states.

  ) -(ΔZ) + , m s (0) .

  Euler-Poisson system. Let us consider the Euler-Poisson system (5.1)

Fig. 2 .

 2 Fig. 2. Subsonic steady state; ρ + z at t = 5.

Fig. 3 .

 3 Fig. 3. Subsonic steady state; ρu at t = 5.

6. 2 .

 2 Euler-Poisson system. We solve the Euler-Poisson system (5.1), (5.2) with p(ρ) = κρ γ , (6.6) κ = 1, γ = 1.1, for x ∈ (0, 0.6) and with initial and boundary conditions (6.7) ρ b (x) = 1 if x ∈ (0.1, 0.5), 100 otherwise, u 0 = 0, V = -1.

Fig. 4 .

 4 Fig. 4. Supersonic steady state; ρ + z at t = 1.

Fig. 5 .

 5 Fig. 5. Supersonic steady state; ρu at t = 1.

Fig. 6 .

 6 Fig. 6. Supersonic steady state; ρ + z at t = 25.

Fig. 7 .

 7 Fig. 7. Supersonic steady state; ρu at t = 25.

Fig. 8 .

 8 Fig. 8. Transcritical flow with shock; ρ + z at t = 200.

Fig. 9 .Fig. 10 .Fig. 11 .

 91011 Fig. 9. Transcritical flow with shock; ρu at t = 200.

for any U l and U r (with nonnegative densities ρ l , ρ r ). This implies in particular that (3.50)

According to (3.28), one has ρ * l u * l = ρ l u l , ρ * r u * r = ρ r u r , and

proving that the scheme preserves the nonnegativity of ρ by interface. The associated speeds involved in the CFL condition are σ l (U * l , U * r ), σ r (U * l , U * r ). In order to prove (ii), consider left and right states U l , U r such that the steady state equations (2.10) are satisfied,

and such that both are subsonic:

then ΔZ = 0 and according to the remark above one has (2.18). Assume now that U l = U r . Then ΔZ = 0; otherwise (3.52) would give two subsonic solutions to a system (3.18), which is not possible by Proposition 3.5. Consider first the case when ΔZ > 0. Then according to Lemma 3.8 we have

According to the remark after (3.2), the relations (2.18) are satisfied. The case ΔZ < 0 is similar, with U * l = U * r = U l . This proves (ii). Property (iv) follows from Lemmas 3.1, 3.2, and 3.11. It remains to prove the consistency (iii). The first property (2.16) is obvious according to the remark above on the case ΔZ = 0. Regarding (2.17), taking into account (3.1), we have to prove that

as U l , U r → U and ΔZ → 0. As stated, we consider only the case when U is not sonic. Let us assume that ΔZ ≥ 0, the complementary case being similar. Then U * r = U r and T r = 0. Therefore, (3.66)

which yields (iii).

Comments on the consistency at sonic points.

The proof of (iii) in the previous theorem when we are close to a sonic point (ρ, u) involves some problems. Consider the case when U l , U r → U , ΔZ → 0, with u 2 = p (ρ) = 0. Assume as previously that ΔZ ≥ 0. If (ρ l , u l ) is supersonic, then (3.59) is valid.

Otherwise, if (ρ l , u l ) is subsonic and

, the computations made in the proof of consistency in the subsonic case can be followed, with the result

Assume now that either (ρ l , u l ) is sonic, or (ρ l , u l ) is subsonic with

In any case we have

Thus, assuming that the numerical flux F is Lipschitz continuous, + O(ΔZ), and using calculations analogous to the ones applied in the proof of the consistency in the subsonic case, one gets


Friedel Hartmann
email: friedelhartmann@uni-kassel.de

Green's functions, boundary elements and finite elements

Keywords: Green's function, boundary elements, finite elements, stiffness matrices, adaptive refinement

Introduction

The Green's function G(y, x) of the Poisson equation

-Δu = p in Ω u = 0 on Γ (1)
allows to write the solution u(x) in terms of an influence function

u(x) = Ω G(y, x) p(y) dΩ y . ( 2
)
In boundary element analysis we replace the Green's function by a fundamental solution g(y, x) Solving the problem [START_REF] Tottenham | Basic Principles[END_REF] with finite elements means that we project the exact solution u(x) onto the solution and trial space

u(x) =
V h = {v h ∈ H 1 (Ω)|v h (x) = 0, x ∈ Γ} (5)
that is the FE-solution u h is the best approximation in the sense of the strain energy metric

a(u -u h , u -u h ) ≤ a(u -v h , u -v h ) v h ∈ V h (6)
where

a(u, v) := Ω (∇u • ∇v) dΩ . ( 7
)

Tottenham's equation

It seems that these two methods follow two totally different lines of reasoning but the surprising result is that also the FE-method is a Green's function method. This result goes back to Tottenham who in a very early paper, [START_REF] Tottenham | Basic Principles[END_REF], mentioned that the FE-solution u h has the form

u h (x) = Ω G h (y, x) p(y) dΩ y (8)
where G h (y, x) is the FE-Green's function that is the projection of the exact Green's function

-Δ y G(y, x) = δ(y -x) G(y, x) = 0 y ∈ Γ ( 9 )
onto the space V h . If the nodal shape functions ϕ i (x) form a basis of V h then G h (y, x) is the function

G h (y, x) = i u i (x) ϕ i (y) (10)
where the nodal values u i = u i (x) are the solution of the system

K u(x) = f k ij = a(ϕ i , ϕ j ) f i = ϕ i (x) . ( 11
)
The proof of Tottenham's equation ( 8) is easy and can be done on one line. Namely we have

u h (x) = Ω δ(y -x) u h (y) dΩ y = a(G h , u h ) = Ω G h (y, x) p(y) dΩ y . ( 12
)
First u h is considered to be the FE solution of the right-hand side p and G h ∈ V h assumes the role of a virtual displacement

a(G h , u h ) = Ω G h (y, x) p(y) dΩ y (13)
next G h is considered to be the FE solution of the right-hand side δ, and u h assumes the role of a virtual displacement

u h (x) = Ω δ(y -x) u h (y) dΩ y = a(G h , u h ) (14)
which explains the left-hand side. The symmetric strain energy a(G h , u h ) plays the role of a turnstile.

The inverse stiffness matrix

The approximate Green's function for the displacement u(x) at a node x k has the form

G h (y, x k ) = i u Gi (x k ) ϕ i (y) . ( 15
)
The vector u G = {u G1 , u G2 , . . . , u Gn } is the solution of the n × n system

K u G = e k (unit vector e k ) , (16)
which means that the columns c k of the inverse stiffness matrix K -1

u G = K -1 e k = c k (17)
are the nodal displacements which belong to the n Green's functions

G h (y, x k ) of the n nodes x k . G h (y, x k ) = i c ki ϕ i (y) = c T k ϕ i (y) . ( 18
)
This explains why the inverse of a tridiagonal matrix is fully populated. Even if only one node x k carries a point load P = 1 the whole structure deforms. An example provides the bar in Fig. 1 a, which consists of five linear elements. The stiffness matrix and its inverse are, see Fig. 1 b

K = EA l ⎡ ⎢ ⎢ ⎣ 2 -1 0 0 -1 2 -1 0 0 -1 2 -1 0 0 -1 2 ⎤ ⎥ ⎥ ⎦ , ⇒ K -1 = l 5 EA ⎡ ⎢ ⎢ ⎣ 4 3 2 1 3 6 4 2 2 4 6 3 1 2 3 4 ⎤ ⎥ ⎥ ⎦ . ( 19
)
And evidently are the columns of the inverse the Green's functions of the nodes.

Boundary elements

When the BE solution (4) is compared with the FE solution (8) there seems not to be much agreement between the two solutions. But the BE solution (4) is just the same formula

u h (x) = Ω G h (y, x) p(y) dΩ y (20)
in disguise, of course with a different G h . To see this note that the Green's function G can be split into a fundamental solution G and a regular part u R ,

u(x) = Ω G(y, x) p(y) dΩ y = Ω g(y, x) p(y) dΩ y + Ω u R (y, x) p(y) dΩ y . ( 21
)
The regular part is the solution of the problem

-Δu R = 0 in Ω u R (y, x) = -g(y, x) y ∈ Γ . ( 22
)
This splitting implies that the boundary integrals in the influence function (4) are just an equivalent expression for the work done by the distributed load p acting through the regular part u R :

Ω u R (y, x) p(y) dΩ y = Γ G(y, x) ∂u ∂n (y) ds y - Γ ∂G ∂ν (y, x) u(y) ds y . ( 23
)
Hence it can be assumed that the boundary integrals in the BE solution (4) play the same role,

Ω u h R (y, x) p(y) dΩ y : = Γ G(y, x) ∂u h ∂n (y) ds y - Γ ∂G ∂ν (y, x) u h (y) ds y , ( 24
)
and so we arrive at (20), which makes the two methods look alike. Though there is a difference in the Green's function: the FE Green's function is only mesh dependent, while the BE Green's function is also load-case dependent.

Note that the BE method uses an ingenious approach: it does not approximate u R (y, x)-this would be simply too laborious, because at every point x it would have to approximate a new function u R (y, x)-but instead substitutes for the work integral (u R [x], p) the work done by the Cauchy data u on Γ N and t = ∂u/∂n on Γ D via the conjugate terms of the fundamental solution G. That is, the program knows that it suffices to approximate the "static" data u and t leaving the effects caused by a change in the observation point x to the fundamental solution G(y, x). This is the essence of (24).

Accuracy

The FE method has a handicap because it must approximate both parts of the Green's function G = g + u R while in the BE method the fundamental solution is built into the code. This handicap can be overcome if the FE-method approximates only the regular part

u(x) -u h (x) = Ω (G(y, x) -G h (y, x)) p(x) dΩ y = Ω (g(y, x) + u R (x) -g(y, x) + u h R (y, x)) p(y) dΩ y = Ω (u R (x) -u h R (y, x)) p(y) dΩ y . ( 25
)
Of course if the source point x comes too close to the boundary then the regular part will also become singular.

is even questionable whether the Dirac delta is an integrable function, whether the expressions make sense. We simply assume they do.

Next we only need to find the function G which has the delta function as its right-hand side, that is which solves the equation

-Δ G(y, x) = δ Σ (y -x) G = 0 on Γ (30)
because Green's second identity then implies

B(G, u) = Ω δ Σ (y -x) u(y) dΩ y J(u) - Ω G(x, y) p(y) dΩ y = 0 (31)
or

J(u) = Ω G(x, y) p(y) dΩ y . ( 32
)
This is the 'trick'. The abstract Dirac delta has got a shape, G, and the scalar product of this G with the right hand side p renders the value J(u). Now we do finite elements. We approximate the Green's function G with the n shape functions in

V h G h (y, x) = n i=1 u i (x) ϕ i (y) (33)
and so the variational problem for the weak solution becomes on V h a system of n equations

a(G h , ϕ i ) = (δ Σ (y -x), ϕ i ) = J(ϕ i ) i = 1, 2, . . . n (34)
or in matrix notation

K u(x) = f (x) k ij = a(ϕ i , ϕ j ) f i = J(ϕ i ) . ( 35
)
So the extension of the FE method to Green's functions is simple and straightforward.

Infinite energy

Strictly speaking the extension is not so simple because most Green's functions have infinite energy. Hence approximating a Green's function with finite elements is an ill-posed problem. The surprising aspect is that the FE-method does not care. All the input in an FE-program is processed by these approximate Green's function which theoretically do not lie in the solution space. But it is remarkable how 'old-fashioned' the ansatz of the FE-method is. The FE-method is a numerical implementation of Green's function. And this is the same idea behind the boundary element method. So in this sense the two methods are identical twins.

Goal oriented refinement

The error of an FE-solution

u(x) -u h (x) = Ω (G(y, x) -G h (y, x)) p(x) dΩ y (36)
is simply proportional to the error in the approximate Green's function. So given that we are interested in the value u(x) of the solution at a particular point x we must refine the mesh in such a way that the error in the Green's function for the point value u(x) is minimized. This technique is called goal oriented adaptive refinement.

By making use of the Galerkin orthogonality and Schwarz' inequality we can formulate the following estimate for the error of an FE-solution

|u(x) -u h (x)| ≤ || G[x] -G h [x]|| E ||u -u h || E (37)
where ||.|| E is the standard energy norm, ||u|| 2 E := a(u, u) and G[x] := G(y, x) is the Green's function. Again we face the same problem as before: theoretically such an estimate requires the Green's function to have bounded energy what normally is not the case. But if we read u(x) as the average value of the potential over a small disk centered at x the corresponding Green's function G(y, x) would have finite energy and then (37) would make sense. This result (37) is motivation to minimize the error in the Green's function and the error in the solution u h simultaneously. That is, at each refinement step an error indicator η e G for the Green's function and an error indicator η e p for the original problem is calculated on each element and the combined error indicator on each element is then η e = η e G • η e p so that the sum of the local errors η e = η e G • η e p provides an upper bound for the error

|u(x) -u h (x)| ≤ e η e G • η e p . ( 38
)
The energy norms of the errors e G = G-G h and e u = u-u h are calculated by measuring the eigenwork done by the element residual forces and jump terms on the element edges, so that for example

η p = a(e u , e u ) 1/2 Ω e = r(e u ) 1/2 . ( 39
)
This technique was applied to the plate in Fig. 2. The first mesh, the mesh in Fig. 2 a, is the result of a standard adaptive refinement

η e p ≤ ε T OL . ( 40
)
To push the error below the preset error margin the program has to refine all those elements-in practice only the first, say, 30%-where the error η e p of the original (or primal) problem exceeds this margin. The mesh in Fig. 2 b is based on weighting the primal error η p with the error η G of the Green's function elementwise

η e G • η e p ≤ ε T OL . ( 41
)
In most parts of the mesh the error η e G of the numerical Green's function is so low that this inequality is automatically satisfied. That is many of the refinements in Fig. 2 a are not necessary if we are only interested in the stress σ yy at x.

  x) u(y)] ds y + Ω g(y, x) p(y) dΩ y (3) and we determine the unknown boundary values by solving an integral equation approximately and so construct the BE-solution u h (x) = Γ [g(y, x) ∂u h ∂n (y) -∂g ∂n (y, x) u h (y)] ds y + Ω g(y, x) p(y) dΩ y . (4)

Figure 1 :

 1 Figure 1: a) Elastic bar subdivided into five linear elements b-e) the displacements are the columns of the inverse stiffness matrix (all values times l/(5 EA)).

Figure 2 :

 2 Figure 2: Adaptive refinement a) standard refinement η p ≤ ε T OL b) goal oriented refinement η p ×η G ≤ ε T OL

Everything is a functional

The importance of Tottenham's insight lies in the fact that in linear analysis every value, every displacement, every stress, every potential value can be considered a functional J(u). For example the following expressions

are all functionals. And so as we associate with the classical functional, J(u) = u(x), a Dirac delta

so we can associate with every other linear functional J(u) a certain Dirac delta. Let the functional J(u) be for example the average value of the gradient of u over a small circular region Ω c (x) with radius c centered at a point x ∈ Ω

We postulate that this functional can be expressed as

We do not say what the Dirac delta δ Σ (y -x) looks like, what it is. We only say that the scalar product (= integral) of δ Σ with the function u provides the value J(u). In more complicated cases it


Evgeniya G Gerasimova-Chechkina

Brian Toner

Zach Marin

Benjamin Audit

Stéphane G Roux

Françoise Argoul

André Khalil

Olga Gileva

Oleg Naimark

Alain Arneodo

Combining Multifractal Analyses of Digital Mammograms and Infrared Thermograms to Assist in Early Breast Cancer Diagnosis

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

According to World Health Organization, breast cancer is one of the most common causes of cancer deaths among women worldwide. Clinical studies have demonstrated that patients' survivability increases if breast anomalies are detected at early stages [1]. Although X-ray mammography is currently the most effective way of early breast cancer detection, it remains difficult to interpret due to high tissue variability and to two-dimensional projection effects [2,3]. The miss rate in mammography is increased in dense breasts of young women [4,5]. Criticism of the use of screening mammography due to over-diagnosis made clinicians use additional modalities such as ultrasound, magnetic resonance imaging, IR imaging, etc. Since the original observation by R. Lawson [6] that skin temperature over a malignant tumor is higher than its neighborhood, possibly resulting from abnormal increase in metabolic activity and vascular circulation in the tissues beneath [7][START_REF] Keyserlingk | Functional infrared imaging of the breast: Historical perspectives, current application and future considerations[END_REF][START_REF] Sterns | [END_REF][10], IR thermography has been considered as a promising non-invasive screening method of breast cancer detection [11][START_REF] Amalu | Infrared Imaging of the Breast-An Over View[END_REF][START_REF] Ng | [END_REF]. However, the suitability of IR imaging for routine screening has been highly questioned [14], because of insufficient sensitivity for detection of deep lesions and limited knowledge of the relationship between skin surface temperature distribution and thermal diseases. This growing wave of criticism of breast cancer screening techniques has been strengthened by a similar movement criticizing the increasing use of computer-aided detection (CAD) methods that were shown to decrease specificity [15], and to increase recall rates of healthy women [16][17][18][19] resulting in an increase in unnecessary stress in women. Note that most of these CAD methods [20][21][22][23][24] were designed for texture analysis and feature extraction with the prerequisite that background roughness fluctuations of normal breast mammograms are homogeneous and uncorrelated. In this paper, we review the results of recent and ongoing analyses of dynamic infrared thermograms and X-ray mammograms based on wavelet-based multifractal methodologies [25][26][27].

Application of dynamic IR imaging as a diagnostic tool is based on the detection of variations in temperature rhythms generated by the cardiogenic and vasomotor frequencies [28,29]. Human cardiogenic frequency levels lie between 1 and 1.5 Hz and vasomotor frequencies between 0.1 and 0.2 Hz, with higher harmonic frequency rhythms. But, beyond intensity differences in these rhythms between normal and tumor tissues, there is much more to learn from dynamical IR imaging of breast cancer. Using a 1D wavelet-based multiscale analysis of the temperature intensity fluctuations about these physiologic perfusion oscillations, we propose to characterize the multifractal properties of these temperature fluctuations as a novel discriminating method that can be implemented in early screening procedures to identify women with a high risk of breast cancer [25,26].

The main premise in our investigation of mammograms [27] was related to the suggestion that healthy dense and fatty tissues, indeed display persistent long-range correlations and anti-correlations in mammogram spatial fluctuations, respectively [30,31]. Using a 2D wavelet-based multifractal methodology, we recently observed some clear loss of correlations in mammogram area surrounding the malignant tumor and which strongly correlate with regions of anomalous monofractal (as compared to healthy multifractal) temperature fluctuations previously detected in dynamic IR thermograms [27]. This loss of complexity in temperature time-series and mammographic subimages is likely to be the signature of physiological and architectural alterations of the environment of malignant breast tumors. Besides their potential clinical impact, these results open new perspectives in the investigation and understanding of the cellular mechanisms that underlie the loss of tissue homeostasis, cancer initiation and progression.

THE WAVELET TRANSFORM MODULUS MAXIMA METHOD

The wavelet transform (WT) is a mathematical microscope [32][33][34][35][36][START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] that is well suited for the analysis of complex non-stationary time-series, as those found in genomics [START_REF] Nicolay | [END_REF][39][40] and physiological systems [41][42][START_REF] Goldberger | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Richard | [END_REF], thanks to its ability to filter out low-frequency trends in the analyzed signal. Its two-dimensional (2D) generalization has been exploited for image processing [31,[START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Antoine | Two-Dimensional Wavelets and their Relatives[END_REF] including applications in cellular biology [START_REF] Khalil | [END_REF][47][48][49][50][51] and medicine [30,31,52,53].

The 1D Wavelet Transform (Time-Frequency Analysis)

The WT is a space-scale analysis which consists in expanding a signal in terms of wavelets that are constructed from a single function, the "analyzing wavelet" , by means of translations and dilations. The WT of a function ( ) t is defined as [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]:

0 0 1 [ ]( , ) ( ) d

The 2D Wavelet Transform (Space-Scale Analysis)

With an appropriate choice of the analyzing wavelet, one can reformulate Canny's multiscale edge detection [58] in terms of a 2D wavelet transform [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. The general idea is to start smoothing the discrete image data by convolving it with a filter and then to compute the gradient of the smoothed image. Let us consider two wavelets that are, respectively, the partial derivatives with respect to x and y of a 2D-smoothing function ( , )

x y :

1 ( , ) ( , ) x y

x y x and 2 ( , ) ( , ) . x y

x y y (2)

The WT of I (x, y) with respect to 1 and 2 has two components and therefore can be expressed in a vectorial form [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]:

1 2 2 0 1 2 , 2 0 2 2 1 [ ] ( ) d [ ]( , ) { }, 1 [ ] ( ) d a W I f x a a I a I W I f x a a 0 0 x x x x W x x x x (3)
from which we can extract the modulus and argument of WT:

1 2 2 2 1 / 2 [ ]( , ) [( [ ]( , )) ( [ ]( , )) ] M I a W I a W I a 0 0 0 x x x (4)
and

1 2 [ ]( , ) Arg[ [ ]( , ) [ ]( , )]. A I a W I a iW I a 0 0 0 x x x (5)

The Wavelet Transform Modulus Maxima Method

The wavelet transform modulus maxima (WTMM) method consists in computing the WT skeleton defined, at each fixed scale a, by the local maxima L(a) of the WT modulus M (•, a). In 1D, these WTMM are disposed on curves connected across scales called maxima lines which defines the WT skeleton [54-56, 59, 60]. In 2D, these WTMM lie, at a given scale, on connected chains called maxima chains. When considering the points along these maxima chains where M (•, a) is locally maximum, we define the so-called WTMMM L(a). The WTMMM are then linked through scales to form the WT skeleton [31,[61][62][63]. Along these maxima lines the WTMM (resp. WTMMM) behave as a h(t) (resp. a h(x) ), where h(t) (resp. h(x)) is the Hölder exponent characterizing the singularity of (resp. I) at time t (resp. position x). The multifractal formalism amounts to characterize the relative contributions of each Hölder exponent value via the estimate of the D(h) singularity spectrum defined as the fractal dimension of the set of points t (resp. x) where h(t) (resp. h(x)) = h. This spectrum can be obtained by investigating the scaling behavior of partition functions defined in terms of WT coefficients [31,[54][55][56][57]:

( ) ( ) ( , ) ( , ) ~, q q l L a Z q a M a a ( 6
)
where q . Then from the scaling function (q), D(h) is obtained by a Legendre transform:

D(h) min q [qh (q)]. (7)
In practice, to avoid instabilities in the estimation of the singularity spectrum D(h) through the Legendre transform, we instead compute the following expectation values [31,[54][55][56][57]:

( ) ( , ) ln ( , ) ln( ( , ) ( , , )), l L a h q a Z q a M a W q l a q (8) ( ) ˆ( , ) ln( ( , )) ( , ) ( , , )ln( ( , , )) l L a D q a q Z q a Z q a W q l a W q l a q . ( 9
)
where ˆ( , , ) ( , ) ( , ) W q l a M a Z q a are the equivalent of the Bolzmann weight in the analogy that links the multifractal formalism to thermodynamics [56]. Then from the slopes of h(q, a) and D(q, a) vs ln a, we get h(q) and D(q) and therefore the D(h) and, consequently, singularity spectrum as a curve parametrized by q.

Monofractal versus Multifractal Functions

The (q) data are fitted by the quadratic approximation 2 0 1 2 ( ) 2, q c cq c q where the coefficients 0.

n c

The corresponding singularity spectrum has a characteristic single humped shape

2 0 1 2 ( ) ( ) 2 , D h c h c c where 0 (0) c
is the fractal dimension of the support of singularities of (resp. I), c 1 is the value of h that maximizes D(h) and 2 , c the so-called intermittency coefficient [64], characterizes the width of the D(h) spectrum. Homogeneous (mono) fractal functions, i.e. functions with singularities of unique Hölder exponent H are characterized by a linear (q) curve 1 (

) . h q c H The D(h) singularity then reduces to a single point

2 ( 0) c ( ) 1 D h H (resp.
2) as the signature of no change in WTMM statistics across scales [31, 39, 54-57, 61-63, 65, 66]. A nonlinear (q) is the signature of nonhomogeneous functions exhibiting multifractal properties, meaning that the Hölder exponent h(t) [54][55][56] (resp. h(x) [31,[61][62][63]) is a fluctuating quantity that depends on t (resp. x). Then D(h) has a single humped shape 2 ( 0) c that traduces some change in WTMM statistics across scales [31,55,63,[START_REF] Arneodo | A wavelet based multifractal formalism: Application to DNA sequences, satellite images of the cloud structure and stock market data[END_REF].

RESULTS

Experimental Protocols

In this work [25][26][27], we analyzed the thermograms of 33 patients (from 37 to 83 years old) of Perm Regional Oncological Dispensary with breast cancer (invasive ductal and/or lobular cancer) and 14 healthy volunteers (from 23 to 79 years old) [67]. Tumors (1.2 to 6.5 cm in size) were found at 1 to 12 cm depths. Both breasts of each patient and volunteer were imaged with an IR camera (FLIR SC5000) at 50 Hz frame rate during the 10 min immobile patient phase. Each image set contained 30 000 image frames. The mammographic procedure involved two X-ray images of both compressed breasts of each patient, referred to as the medio-lateral oblique (MLO) and the craniocaudal (CC) views. Spatial resolution of images was 200 pixels per 1 cm. For mammogram analysis, our database includes mammograms collected from 30 of the mentioned 33 patients with breast cancer.

WTMM Multifractal Analysis of Breast Infrared Thermograms

We grouped single-pixel temperature time-series into 8×8 pixel 2 squares covering the entire breast [25,26]. The reported results are averaged over 64 temperature time-series in these 8×8 subareas. We analyzed 1-pixel temperature time-series taken from 8×8 pixels 2 -squares covering patient entire breasts. As expected, these timeseries generally fluctuate at a higher temperature when recorded in the tumor region of a cancerous breast (Fig. 1a) than in a symmetrically localized square on the contralateral unaffected breast (Fig. 1b), as well as on a healthy breast (Fig. 1c). As commonly done for noise signals [START_REF] Nicolay | [END_REF]55] and previously experienced when analyzing rainfall timeseries [68], we applied the WTMM method to the cumulative of these temperature time-series using the secondorder compactly supported version (2) (3) of the Mexican hat [26,69] (hence the singularities with a possible negative Hölder exponent 1 0, h become singularities with cum 0 1 1 h h in the cumulative). We found that the partition functions Z(q, a) (Eq. ( 6)) display nice scaling properties for 1 to 5, q over a range of time-scales limited to [0.43 s, 2.30 s] for linear regression fit estimates in a logarithmic representation. The (q) so-obtained are well approximated by quadratic spectra (Fig. 1d). For the cancerous breast of patient 20, (q) is nearly linear as quantified by a very small value of the intermittency coefficient This signature of monofractality is confirmed, when respectively plotting h(q, a) (Eq. ( 8)) and D(q, a) (Eq. ( 9)) vs 2 log a (see Fig. 2 in Ref. [26]) where the slopes 1 ( ) 0.81 0.01 h q c and ( ) 0.99 0.03, D h do not significantly depend on q, meaning that the D(h) singularity spectrum nearly reduces to a single point 1 ( 0 . 8 1 ) 1 D h c (Fig. 1e). In contrast, the (q) spectra obtained for the contralateral unaffected breast of patient 20 and the right healthy breast of volunteer 14 (Fig. 1d) are definitely nonlinear with a no longer negligible quadratic term 2 0.080 0.001 c and 2 0.069 0.002, c respectively. This hallmark of multifractal scaling is confirmed by the computation of the corresponding D(h) spectra in Fig. 1e that both display a single-humped shape [26].

The results of our comparative wavelet-based multifractal analysis of (cumulative) temperature fluctuations over the 33 cancerous breasts, the 32 contralateral unaffected breasts (no right breast for patient 6) and the 28 volunteer healthy breasts are shown in Fig. 2 [25,26]. In Fig. 2a, the corresponding histograms of c 1 values extend over a rather wide range (Fig. 2b). The small values of c 2 dominate in cancerous breasts compared to normal breasts, confirming that the cancerous breasts are enriched in squares where temperature fluctuations display significantly reduced multifractal properties. This justifies that we considered 2 0.03 c as the threshold below (resp. above) which a square was qualified as monofractal (resp. multifractal) (Figs. 3d-3f).

According to the proposed classification criteria, 8×8 pixel 2 squares covering breast thermograms were color coded in the following way: squares with 2 0.03 c were shown in black, squares with 2 0.03 c in white and squares without scaling in gray. A high percentage of squares in the cancerous breast displays monofractal temperature fluctuations with small intermittency coefficient values (49.7% in Fig. 3a), while only few of those squares are found in the contralateral unaffected breast (7.7% in Fig. 3b) and in the healthy volunteer breast (11% in Fig. 3c). Consistently, in the thermograms of the contralateral unaffected and volunteer's normal breasts, a majority of squares (89.4% and 65%, respectively) manifests multifractal scaling with an intermittency coefficient 2 0.03 c (Figs. 3b and3c). Note that 43.1% of the squares in the cancerous breast also display multifractal temperature fluctuations, as observed for normal breasts (Fig. 3a). These squares indeed cover regions of the breast that are far from the tumor area (left upper quadrant) mostly covered by monofractal square. Let us point out that for each breast, a small percentage (<20%) of (gray) squares was removed from our analysis because of lack of scaling.

A comparative statistical analysis of monofractal and multifractal square distributions in breast thermograms confirmed that the percentage of monofractal (black) squares covering cancerous breasts (26.8 ± 3.5%) is about twice higher than the ones covering contralateral unaffected breasts (13.1 ± 2.3%) and breasts of healthy volunteers (11.3 ± 2.2%). This excess is compensated by a smaller percentage of multifractal (white) squares in cancerous breasts (56.9 ± 4.4%) than in contralateral unaffected breasts (68.5 ± 3.8%) and normal breasts (69.4 ± 4.3%) in healthy volunteers (Table 1).

We have compared the percentages of monofractal squares in cancerous and contralateral breasts [26] and found that 25 patients have higher quantity of monofractal squares in the cancerous breast. 18 cancerous breasts have a percentage of monofractal squares greater than the mean value (26%). Among the other 15 cancerous breasts, 7 have a smaller percentage of monofractal squares but these squares are well localized in the tumor region. The remaining 8 cancerous breasts were detected as false negatives due to small percentage of monofractal squares and their distant location with regard to the tumor region. Four of them correspond to rather deep tumors (6 cm, 7 cm, 8 cm and 12 cm depth) in the breasts with massive adipose tissue that can explain the absence of significant changes in skin temperature dynamics. Importantly, in five of the contralateral unaffected breasts a high (>26%) percentage of monofractal squares was detected which could be an early sign of pathology and require thus further check-ups. As a control, we reproduced this comparative analysis on both normal breasts of the 14 healthy volunteers. Only 4 volunteers manifested high (>26%) percentage of monofractal squares, while in most of them this value did not exceed 10% [26].

FIGURE 3.

Breast-wide multifractal analysis of skin temperature temporal fluctuations. As estimated from the (q) spectra computed with the WTMM method (Fig. 1d), the 8×8 pixels 2 squares covering the entire breast are coded (a-c) and represented as dots in the To evaluate the efficiency of the test based on the multifractal analysis of breast infrared thermograms, we computed sensitivity (Se) and specificity (Sp). The results reported in Table 2 for our rather small set of patients, namely 76% sensitivity and 86% specificity, are quite encouraging for future applications to larger cohorts of females with breast cancer.

WTMM Multifractal Analysis of X-ray Mammograms

We divided each mammogram view into 360×360 pixel 2 squares spanning 18×18 mm 2 [27]. From the central 256×256 pixels 2 part of the WT skeleton, we computed the partition functions Z(q, a) (Eq. ( 6)) that were found to display scaling properties over a range of scales 2 1 log 3, a corresponding to [0.7 mm, 2.8 mm], for linear regression fit estimates in a logarithmic representation. Within the range 1 3, q statistical convergence was achieved and the so-obtained (q)-spectrum was found remarkably linear (Fig. 4e) for a great majority of squares in both the cancerous and contralateral unaffected breasts, as the signature of monofractal roughness fluctuations [31,61]. We systematically run our 2D WTMM algorithm for all mammographic views in our database. We confirmed the result of our preliminary study [30] that normal tissue regions display monofractal scaling properties

1 2 ( , 0) c H c
as characterized by the Hurst exponent values H < 0.5 in fatty areas which look like anti-persistent self-similar random surfaces, while H > 0.5 in dense areas which exhibit persistent long-range correlations in mammographic roughness fluctuations [31,61]. Importantly, our multiscale method provides an objective way to identify areas with some clear loss of correlations in roughness fluctuations ( 0.5).

H

We extended our waveletbased multifractal analysis of CC and MLO mammographic images to the entire sets of 256×256 pixels 2 squares that cover each breasts of the 30 patients. As illustrated in Figs. 4a-4d for patient 20, we color-coded each of these squares according to the monofractal diagnosis obtained with the 2D WTMM method. Besides confirming the existence of areas of correlated spatial roughness fluctuations, our methodology reveals that peculiar areas of uncorrelated roughness fluctuations are mainly found in the neighborhood of malignant tumors and strongly correlate with the regions of anomalous monofractal (as compared to healthy multifractal) temperature fluctuations previously detected in dynamic IR thermograms (Fig. 3a).

We computed the percentages of monofractal uncorrelated H = 0.5 squares in the mammograms of the two breasts of 30 patients with breast cancer. Good agreement was observed between the percentages of uncorrelated H = 0.5 squares in the CC and the corresponding MLO mammograms of both breasts of each patient (Fig. 5a). Obviously, the percentage of uncorrelated H = 0.5 squares in the cancerous mammograms is much larger than in the opposite breast mammograms. When comparing the two breasts of patients with breast cancer (Fig. 5b), we definitely evidenced some important asymmetry with a significant excess of H = 0.5 squares in the mammograms of the cancerous breast compared to the contralateral unaffected breast. This is a very promising indication that the loss of correlations in the mammographic spatial density fluctuations can be used as an indicator of malignancy. We correlated our findings in IR thermograms and X-ray mammograms of the 30 patients (Fig. 5c). Z q a vs. It is clear that the majority of patients have consistent large numbers of "risky" uncorrelated H = 0.5 squares in mammograms and of "risky" monofractal squares in thermograms of their cancerous breast as compared to their contralateral unaffected breast. We further investigated the spatial distribution of uncorrelated H = 0.5 squares over the entire breast defining a H = 0.5 cluster as an aggregation of H = 0.5 squares sharing a common edge [27]. Isolated H = 0.5 squares were considered as cluster-singlet. As reported in Table 3, the mean size of a cluster in a cancerous breast is larger than in the contralateral unaffected breast. This clustering effect is even more evident for the size of the largest H = 0.5 square that is twice bigger for the cancerous breasts than for the opposite ones. These findings suggest that uncorrelated H = 0.5 squares conglomerate in rather big clusters in cancerous breasts likely defining "risky" area that can be used as a new indicator of malignancy.

CONCLUSIONS

By combining a 1D WTMM multifractal analysis of IR temperature time-series and a 2D WTMM multifractal analysis of X-ray mammograms, we have revealed some change in both temporal and spatial scaling properties in the region of the malignant tumor in cancerous breasts. This is a clear indication of some physiological alterations and architectural disorganization of the "niche" surrounding breast tumor that may be the signature of a transition from an antitumorigenic to a tumorigenic microenvironment underlying the process of oncogenic transformation, tissue invasion and metastasis invasion during cancer progression [70][71][START_REF] Maguer-Satta | The stem cell niche: the black master of cancer[END_REF][START_REF] Lu | [END_REF]. Despite the exploratory nature of this preliminary study involving a rather small cohort of females who all went through surgery to remove the histologically confirmed malignant tumor, the obtained results are quite encouraging in the perspective of generalizing our approach to the analysis of longitudinal studies. This will allow us to quantify the statistical performance (sensitivity, specificity) of our computational approach and its potentiality to assist in early breast cancer diagnosis, progression and treatment.

FIGURE 1 .FIGURE 2 .

 12 FIGURE 1. Multifractal analysis of cumulative IR temperature time-series in a 8×8 pixels 2 square in the cancerous right breast (•) and opposite left breast ( ) of patient 20 and in the healthy right breast ( ) of volunteer 14. (a-c) 1 min portions of pixel temperature time-series in the cancerous breast (top), contralateral unaffected breast (median) and healthy volunteer breast (bottom), (d) (q) vs q estimated by linear regression fit of log 2 Z(q, a) vs 2 log a over a range of time-scales [0.3, 3] s, (e) D(h) vs h. Adapted from Gerasimova et al.[26]

  FIGURE 3.Breast-wide multifractal analysis of skin temperature temporal fluctuations. As estimated from the (q) spectra computed with the WTMM method (Fig.1d), the 8×8 pixels 2 squares covering the entire breast are coded (a-c) and represented as dots in the 1 2 ( , ) c c plane (d-f). The coding has the following meaning: black: 2 0.03, c white: 2 0.03 c and gray: no scaling (see text). (a, d) Cancerous right breast of patient 20: the tumor is located in the upper outer quadrant, (b, e) contralateral unaffected left breast of patient 20, (c, f) healthy right breast of volunteer 14. Adapted from Gerasimova et al.[26]

FIGURE 4 .

 4 FIGURE 4. Wavelet-based multifractal segmentation of X-ray mammograms. Patient 20: cancerous right breast (a, b) and contralateral unaffected left breast (c, d). (a, c) As estimated from the (q) spectra of CC mammographic view computed with the 2D WTMM method, 360×360 pixel 2 squares spanning 18×18 mm 2 were colored coded according to monofractal H < 0.45 (blue), 0.45 H 0.55 (yellow), H > 0.55 (red) and no scaling (pink). (b, d) Same as (a, c) for MLO mammographic view. (e) (q) vs. q estimated by linear regression fit of 2 log ( , )Z q a vs.

2 logFIGURE 5 .

 25 FIGURE 5. Analysis of the percentage of "risky" monofractal uncorrelated H = 0.5 squares in mammograms and "risky" monofractal squares in thermograms of the breasts of 30 patients with breast cancer. (a) Percentage of H = 0.5 squares in cancerous (•) and contralateral unaffected ( ) breasts: MLO view vs. CC view. (b) Mean percentage (MLO + CC)/2 of H = 0.5 squares in MLO and CC views: contralateral unaffected breast vs. cancerous breast. (c) Percentage of monofractal squares in thermograms vs. percentage of H = 0.5 squares in mammograms of the cancerous ( ) and contralateral unaffected ( ) breasts of the 30 patients with breast cancer. Adapted from Gerasimova et al. [27]

6 Cancerous breast Contralateral unaffected breast Normal breast 31 32 33 34 35 36 28 29 30 31 32 33 0.12 D E F 0.09

	C	0.06
		0.03
		0
		0.5	0.9	1.3	1.7	0.5	0.9	1.3	1.7	0.5	0.9	1.3	1.7
			C 1				C 1				C 1

TABLE 1 .

 1 Percentage

		of monofractal 2 ( c	0.03),	multifractal 2 ( c	0.03)	and no-scaling 8×8 pixels 2
	squares in cancerous, contralateral unaffected and healthy breasts [26]
		Cancerous breasts	Opposite breasts	Healthy breasts
	Monofractal	26.8 ± 3.5		13.1 ± 2.3		11.3 ± 2.2
	Multifractal	56.9 ± 4.4		68.5 ± 3.8		69.4 ± 4.3
	No-scaling	16.3 ± 2.9		18.4 ± 2.8		20.0 ± 5.0

TABLE 2 .

 2 Sensitivity Se = TP/(TP+FN) and specificity Sp = TN/(TN+FP) of the test based on multifractal analysis of breast infrared thermograms

		Has the disease	Does not have the disease
	Positive test score	True Positives (TP)	False Positives (FP)
		25	4
	Negative test score False Negatives (FN)	True Negatives (TN)
		8	24

TABLE 3 .

 3 Mean values of H = 0.5 square clusters characteristics

		Cancerous	Contralateral
		breasts	unaffected breasts
	Size of a cluster	3.7	2.3
	Size of the largest cluster	8.9	4.0

ACKNOWLEDGMENTS

This work was supported by INSERM, ITMO Cancer for its financial support under contract PC201201-084862 "Physiques, mathématiques ou sciences de l'ingénieur appliqués au Cancer", the Russian Foundation for Basic Research (grant No. 16-41-590235) and the Maine Cancer Foundation.

The study reported in this article was conducted according to accepted ethical guidelines involving research in humans and/or animals and was approved by an appropriate institution or national research organization.

The study is compliant with the ethical standards as currently outlined in the Declaration of Helsinki. All individual participants discussed in this study, or for whom any identifying information or image has been presented, have freely given their informed written consent for such information and/or image to be included in the published article.