A nonlinear weakly dispersive method for recovering the elevation of irrotational surface waves from pressure measurements

Keywords: surface wave, measurements, nonlinear, weakly dispersive, nonhydrostatic, shallow water, LiDAR, Acoustic Surface Tracking

We present the derivation of a nonlinear weakly dispersive formula to reconstruct, from pressure measurements, the surface elevation of nonlinear waves propagating in shallow water. The formula is simple and easy to use as it is local in time and only involves first and second order time derivatives of the measured pressure. This novel approach is evaluated on laboratory and field data of shoaling waves near the breaking point. Unlike linear methods, the nonlinear formula is able to reproduce at the individual wave scale the peaked and skewed shape of nonlinear waves close to the breaking point. Improvements in the frequency domain are also observed as the new method is able to accurately predict surface wave elevation spectra over four harmonics. The nonlinear weakly dispersive formula derived in this paper represents an economic and easy

Introduction

Near-bottom-mounted pressure sensors have long been used for measuring surface wave in the nearshore. However, the relationship between bottom pressure and sea surface elevation is not straightforward. This relationship is commonly assumed to be given by linear wave theory, the so-called transfer function method (e.g. Bishop and Donelan [START_REF] Bishop | Measuring waves with pressure transducers[END_REF] and Tsai et al. [START_REF] Tsai | On the recovery of surface wave by pressure transfer function[END_REF]). The validity of this linear reconstruction has been extensively studied in field conditions for waves propagating in relatively shallow water [START_REF] Hom-Ma | Response characteristics of underwater wave gauge[END_REF][START_REF] Esteva | Comparison of pressure and staff wave gage records[END_REF][START_REF] Cavaleri | Measurement of the Pressure and Velocity Field Below Surface Waves[END_REF][START_REF] Guza | Local and shoaled comparisons of sea surface elevations, pressures, and velocities[END_REF]. Although discrepancies were greater close to the break point, Guza and Thornton [START_REF] Guza | Local and shoaled comparisons of sea surface elevations, pressures, and velocities[END_REF] found a good agreement in and outside the surf zone between sea surface elevation spectra derived from pressure data and from direct elevation measurements. Errors in both total variance and energy density in a particular frequency band were less than 20%. In a more controlled environment, Bishop and Donelan [START_REF] Bishop | Measuring waves with pressure transducers[END_REF] estimated that using linear wave theory was leading to error of about 5% of the wave height; uncertainty in the deployment of in situ instruments and the data itself was thought to be responsible for the varying error estimates found in the literature. Following these seminal studies, the linear reconstruction method has become the main approach for characterizing shallow-water surface-wave elevation in field conditions. This approach is commonly used for determining bulk wave parameters such as the significant wave height H s , but it has also served as a basis for studying nonlinear wave interactions in the field (e.g. Elgar and Guza [START_REF] Elgar | Observations of bispectra of shoaling surface gravity waves[END_REF], Elgar et al. [START_REF] Elgar | Spectral evolution of shoaling and breaking waves on a barred beach[END_REF], Senechal et al. [START_REF] Senechal | Field experiment on secondary wave generation on a barred beach and the consequent evolution of energy dissipation on the beach face[END_REF], Henderson et al. [START_REF] Henderson | Nonlinear generation and loss of infragravity wave energy[END_REF]). However, we know that wave nonlinearities can be strong in the shoaling zone, especially in the region close to the onset of breaking, and thus the use of a linear theory to reconstruct wave elevation can be questioned. For instance, Bonneton and Lannes [START_REF] Bonneton | Recovering water wave elevation from pressure measurements[END_REF] and Martins et al. [START_REF] Martins | On the influence of swash-based reflection on surf zone hydrodynamics: a wave-by-wave approach[END_REF] have shown that the linear reconstruction fails to describe the peaky and skewed shape of nonlinear waves, and lead to an underestimation of the individual wave height by up to 30% just prior the breaking point [START_REF] Martins | On the influence of swash-based reflection on surf zone hydrodynamics: a wave-by-wave approach[END_REF].

Such measurement errors are problematic for many coastal applications, such as studies on wave overtopping and submersion which require accurate measurements of the highest wave crests. Furthermore, a correct description of wave asymmetry and skewness is of paramount importance for understanding sediment dynamics (e.g. Dubarbier et al. [START_REF] Dubarbier | Processbased modeling of cross-shore sandbar behavior[END_REF]). Finally, an accurate description of the wave elevation field is also crucial for the validation of the new generation of fully-nonlinear phase-resolving wave models (e.g. Zijlema et al. [START_REF] Zijlema | SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters[END_REF], Bonneton et al. [START_REF] Bonneton | Recent advances in Serre-Green Naghdi modelling for wave transformation, breaking and runup processes[END_REF] or Shi et al. [START_REF] Shi | A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF]).

Even if some methods are now available for a direct measurement of the surface elevation, such as acoustic surface tracking [START_REF] Birch | The capabilities of Doppler current profilers for directional wave measurements in coastal and nearshore waters[END_REF] or LiDAR scanning [START_REF] Martins | Monitoring Individual Wave Characteristics in the Inner Surf with a 2-Dimensional Laser Scanner (LiDAR)[END_REF], pressures sensors remain a very useful tool for coastal wave applications. Indeed, they are cheap, robust, not sensitive to air bubbles or turbidity, and easy to deploy since they do not require the presence of nearshore infrastructure, as it can be the case for LiDAR technology [e.g., see 21]. Bonneton and Lannes [START_REF] Bonneton | Recovering water wave elevation from pressure measurements[END_REF] recently derived a method which allows a fully dispersive nonlinear reconstruction of the surface elevation from pressure measurements. Comparisons with numerical Euler solutions and laboratory data showed that this nonlinear method provides much better results than the classical linear approach. It gives an accurate prediction of the maximum elevation and, contrary to the nonlinear heuristic method proposed by Oliveras et al. [START_REF] Oliveras | Recovering the water-wave profile from pressure measurements[END_REF], it accurately reproduces the skewed shape of nonlinear dispersive wave fields. However, this method requires, like the classical linear transfer approach [START_REF] Bishop | Measuring waves with pressure transducers[END_REF][START_REF] Tsai | On the recovery of surface wave by pressure transfer function[END_REF] and the heuristic method [START_REF] Oliveras | Recovering the water-wave profile from pressure measurements[END_REF][START_REF] Vasan | A method to recover water-wave profiles from pressure measurements[END_REF], the use of a frequency cutoff which becomes a limiting factor for the reconstruction of strongly nonlinear waves. In the present paper we derive a nonlinear weakly-dispersive method which allows an accurate reconstruction of nonlinear waves in shallow water, especially just prior to breaking.

2a x z z = 0 z = -h 0 δ m z = ζ(x,t)
Figure 1: Sketch of the Cartesian coordinate system: x is the horizontal axis along which waves propagate and z points vertically upwards, with z = 0 being the mean water level and z = -h 0 the distance to the bottom. The wave amplitude is noted a and δm represents the distance to the bottom at which the pressure is measured.

Nonlinear weakly-dispersive reconstruction formula

In this section we derive a formula working in the time domain, which allows the elevation reconstruction of nonlinear shallow water waves from pressure measurements. This formula is an approximate expression, in the shallow water regime, of the fully dispersive formula derived in [START_REF] Bonneton | Recovering water wave elevation from pressure measurements[END_REF]. The derivation presented in this section is much simpler and straightforward compared to the general fully dispersive derivation.

We consider that the wave field is locally close to a two-dimensional wave field. We choose Cartesian coordinates (x, z), where x is the horizontal axis along which waves propagate and z the upward vertical coordinate. We denote z = ζ(x, t) the elevation of the free surface above the still water level z = 0, and by z = -h 0 the constant bottom elevation (see Figure 1). The water depth h can be expressed as h(x, t) = h 0 + ζ(x, t). The pressure P m is measured at a distance δ m above the bottom, P m = P | z=-h 0 +δm , where P (x, z, t) is the pressure field.

The fluid motion is governed by the free-surface incompressible irrotational Euler equations; if the flow is irrotational, as it is assumed here, it is conve-nient to work with a velocity potential instead of the velocity field, and with Bernoulli's equation instead of Euler equations. If φ denotes the velocity potential, these equations can be recast in the form:

∂ 2 x φ + ∂ 2 z φ = 0 ∂ t φ + gz + 1 2 |∂ x φ| 2 + 1 2 |∂ z φ| 2 = - 1 ρ (P -P atm ),
where ρ is the water density, g the gravity and P atm the (constant) atmospheric pressure. These equations are complemented by boundary conditions. At the bottom we have

∂ z φ = 0 on z = -h 0 ; ( 1 )
at the surface, we have the classical kinematic equation on ζ,

∂ t ζ = ∂ z φ -∂ x ζ • ∂ x φ on z = ζ, (2)
and the pressure continuity,

P = P atm on z = ζ. ( 3
)
Three main length scales are involved in this problem: the characteristic horizontal length L (L = 1/k, where k is the typical wave number), the amplitude a of the wave, and the depth at rest h 0 . The problem is then controlled by two dimensionless parameters:

ε = a h 0 , μ= h 2 0 L 2 = (h 0 k) 2
where ε is a nonlinearity parameter while μ is the shallowness parameter. The different variables and functions involved in this problem can be put in dimensionless form using the relations

x = x L , z = z h 0 , δ m = δ m h 0 , t = √ gh 0 L t, ζ = ζ a , h = h h 0 = 1 + ζ , φ = h 0 aL √ gh 0 φ, P = P ρgh 0 ,
where the primes are used to denote dimensionless quantities.

Omitting the primes for the sake of clarity, the vertical momentum equation, in dimensionless form, writes

εΓ = -1 -∂ z P,
where Γ = ∂ t w + εu∂ x w + ε μ w∂ z w is the vertical acceleration, u = ∂ x φ the horizontal velocity and w = ∂ z φ the vertical velocity. Integrating this equation over z we get

ζ = ζ H - εζ -1+δm Γdz ( 4
)
where ζ H is the dimensionless hydrostatic reconstruction

ζ H = 1 ε (P m -P atm -1 + δ m ). (5)
The formula ( 4) is exact but involves quantities that cannot be expressed in terms of the measured pressure P m . Our goal is to derive approximate formulas 75 that can be expressed as a function of P m , or equivalently ζ H . Following [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF],

we shall perform an asymptotic expansion of (4) in terms of the shallowness parameter μ.

The velocity potential φ is given at second order by

φ = ψ - μ 2 (z + 1) 2 -h 2 ∂ 2 x ψ + O(μ 2 ), (6)
with ψ = φ | z=εζ . From this equation we can deduce that

u = U + O(μ) w = -μ(z + 1)∂ x U + O(μ 2 ),
where U is the depth-averaged horizontal velocity. Using these relations Eq. ( 4)

becomes

ζ = ζ H - μ 2 (h 2 -δ 2 m ) ∂ 2 x ζ + 2 ∂ x U 2 + O(μ 2 ). ( 7
)
From the linearized water waves equations in shallow water, we know that

∂ 2 t ζ - ∂ 2 x ζ = O(ε, μ
), and the above formula therefore yields

∂ 2 x ζ H = ∂ 2 t ζ H + O(ε, μ). ( 8
)
For weakly nonlinear waves (ε = O(μ)) we deduce a linear shallow water reconstruction formula from ( 7) and ( 8),

ζ SL = ζ H - μ 2 (1 -δ 2 m )∂ 2 t ζ H , (9)
which is valid up to terms of order O(εμ, μ 2 ).

80

For moderately nonlinear waves (ε = O(μ) 1/2 ) the O(εμ) terms cannot be neglected, and we therefore seek a higher order correction of ( 9) under the form

ζ SNL = ζ SL + εμC,
where the corrector term C is chosen in such a way that Eq. ( 7) is satisfied up

to O(μ 2 ) terms if ζ is replaced by ζ SNL . This yields the condition - μ 2 (1 -δ 2 m )(∂ 2 t ζ H -∂ 2 x ζ H ) + εμC = -εμ ζ∂ 2 x ζ + (1 -δ 2 m ) ∂ x U 2 + O(μ 2 )
From the mass conservation equation

∂ t ζ + ∂ x (hU ) = 0, we deduce that ∂ x U = -∂ t ζ + O(ε) which, together with (8) yields ∂ 2 t ζ H -∂ 2 x ζ H = 2ε 1 -δ 2 m C + ∂ t ζ∂ t ζ -δ 2 m ∂ t ζ 2 + O(μ).
The O(ε) terms gathers the quadratic interactions of the wave field. Since at order O(ε, μ), the waves are governed by a linear wave equation with speed 1, these quadratic terms are either the product of co-propagating or counterpropagating waves and therefore of the form

∂ 2 t ζ H -∂ 2 x ζ H = ε F (x -t) + G(x + t) + F (x -t)G(x + t) + O(μ),
for some functions F and G that depend quadratically on ζ and some functions F and G that depend linearly on ζ.

The only possibility that does not create secular growth for ζ H is if the righthand-side is equal to zero (see for instance [START_REF] Lannes | Secular growth estimates for hyperbolic systems[END_REF] or Lemma 7.20 in [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF]). This corresponds to

C = -∂ t ζ∂ t ζ + δ 2 m ∂ t ζ 2 . Since ζ = ζ SL + O(μ),
we therefore obtain the following nonlinear shallow water reconstruction formula

ζ SNL = ζ SL -εμ ∂ t ζ SL ∂ t ζ SL -δ 2 m ∂ t ζ SL 2 , ( 10
)
which is valid up to O(μ 2 ) terms for moderately nonlinear waves. A generalization of this formula in presence of a background current is given in [START_REF] Bonneton | Recovering water wave elevation from pressure measurements[END_REF].

The reconstruction formula [START_REF] Elgar | Spectral evolution of shoaling and breaking waves on a barred beach[END_REF] has been derived under the assumption of a locally flat bottom. We can easily extend this approach by taking into account a slowly varying bathymetry. The bottom is given by z = -h 0 + b(x), where b is a slowly varying function of x. The nonlinear reconstruction becomes in dimensionless form

ζ = ζ SNL + √ μσ b ∂ x b∂ x ζ + O(μ 2 , μσ b ),
where σ b is the characteristic slope of the bathymetry at the measurement location. For unidirectional traveling waves the nonlinear reconstruction simplifies to the form which follows

ζ = ζ SNL - √ μσ b ∂ x b∂ t ζ SNL + O(μ 2 , μσ b ).
The parameter √ μσ b being very small for many coastal applications, as those discussed in section 3, we neglect throughout the paper the bottom contribution.

In variables with dimension the reconstruction formulas ( 5), ( 9) and [START_REF] Elgar | Spectral evolution of shoaling and breaking waves on a barred beach[END_REF] become

ζ H = P m -P atm ρg -h 0 + δ m ( 11
)
ζ SL = ζ H - h 0 2g (1 -(δ m /h 0 ) 2 )∂ 2 t ζ H ( 12
)
ζ SNL = ζ SL - 1 g ∂ t ζ SL ∂ t ζ SL -(δ m /h 0 ) 2 ∂ t ζ SL 2 . ( 13
)
Contrary to fully dispersive reconstructions (see [START_REF] Bonneton | Recovering water wave elevation from pressure measurements[END_REF]) these formulas can be applied locally in time and do not necessarily require Fourier transforms (see Appendix A). The ability of these formulas to reconstruct shallow water waves was first assessed by [START_REF] Bonneton | Recovering water wave elevation from pressure measurements[END_REF] from comparison with solitary wave solutions computed from the Euler equations. In the next section the validation is extended to nonlinear shallow water waves propagating on gently sloping bottoms.

Applications to laboratory and field data

To assess the ability of the formulas derived in the preceding section we use laboratory and field data corresponding to weakly dispersive waves (μ ≤ 0.28)

propagating in the shoaling zone close to the breaking point. Such test cases, which are associated with strong nonlinearities (i.e. large ), represent highly demanding cases for surface wave reconstruction.

The reconstructed elevation is obtained from pressure signals measured under waves and then compared to direct elevation measurements. The application of formulas ( 12) and ( 13) requires to compute first and second order time derivatives of the measured pressure P m , or equivalently ζ H . To compute these derivatives, it is recommended to filter the measurement noise, either using a local filter (e.g. moving average filter) or a spectral filter. In the present study we use a low-pass spectral filter with a cut-off frequency f c . Linear and nonlinear fully dispersive reconstructions (e.g. [START_REF] Bishop | Measuring waves with pressure transducers[END_REF][START_REF] Guza | Local and shoaled comparisons of sea surface elevations, pressures, and velocities[END_REF][START_REF] Bonneton | Recovering water wave elevation from pressure measurements[END_REF][START_REF] Vasan | A method to recover water-wave profiles from pressure measurements[END_REF]) also require introducing a cut-off frequency. However, it is worth noting that this fully-dispersive cut-off is much more restrictive (i.e. much lower f c used) than that required for the present weakly dispersive reconstructions. Indeed, the frequency cut-off used in fully dispersive reconstructions is not fixed by the level of noise in the measurements but by the level of wave nonlinearity and the error originating from the use of the linear dispersion relation for the estimation of the nonhydrostatic correction factor cosh(kh 0 ).

In order to accurately compute the time derivatives involved in equations ( 12) and ( 13) we have chosen to use a Fourier transform method. The use of a local time discretization approach for solving the time derivatives is also discussed in appendix A. For the shallow water applications presented in this paper, the shallow water linear reconstruction [START_REF] Guza | Local and shoaled comparisons of sea surface elevations, pressures, and velocities[END_REF] gives performance at least as good as those given by the fully dispersive linear reconstruction (the classical so-called transfer function formula, see [START_REF] Bishop | Measuring waves with pressure transducers[END_REF]). For that reason, we do not present linear transfer function reconstruction in the following validations.

Monochromatic waves

We first use experimental data from the A7-mono test obtained during the BARDEXII set of experiments [START_REF] Masselink | Large-scale Barrier Dynamics Experiment II (BARDEX II): Experimental design, instrumentation, test program, and data set[END_REF], performed at prototype-scale in the Delta flume, Vollenhove, The Netherlands. Martins et al. [START_REF] Martins | On the influence of swash-based reflection on surf zone hydrodynamics: a wave-by-wave approach[END_REF] used pressure and sur-face elevation data obtained from a LiDAR scanner from this test to validate the simulated wave transformation across the wave flume. Although the datasets are described in Martins et al. [START_REF] Martins | On the influence of swash-based reflection on surf zone hydrodynamics: a wave-by-wave approach[END_REF], some basic information are reminded here.

Second-order monochromatic Stokes waves (T = 12.1 s, H = 0.68 m) were generated in the experimental flume. The pressure was measured just before the break point, with a sampling frequency of 20 Hz. The pressure sensor was situated at a distance δ m = 0.33 m above the bottom. For this wave test, the still water depth was h 0 = 1.17 m, the characteristic bottom slope was σ b = 0.041 and the values of the dimensionless parameters were μ = 0.032 and ε = 0.65. It is worth noting that the pressure data was obtained close to the breaking point [START_REF] Martins | On the influence of swash-based reflection on surf zone hydrodynamics: a wave-by-wave approach[END_REF], where nonlinearities are the strongest and close to a crest antinode, due to the strong reflection experienced during the test.

The LiDAR scanner consisted of a LMS511 SiCK eye-safe (λ = 905 nm)

2D scanner deployed 3.9 m above the still water level, a meter landward of the pressure transducer. Unlike the pressure transducer, the scanner directly detects the surface by capturing back the light beam scattered at the surface by the presence of ripples/roughness or air bubbles in case of breaking. The LiDAR data were collected at 37.5 Hz, and only the data at the pressure transducer location is used. For this specific experiment the LiDAR measurement accuracy was about 0.02 m. The reader is referred to Martins et al. [START_REF] Martins | Monitoring Individual Wave Characteristics in the Inner Surf with a 2-Dimensional Laser Scanner (LiDAR)[END_REF][START_REF] Martins | On the influence of swash-based reflection on surf zone hydrodynamics: a wave-by-wave approach[END_REF] for more details on the scanner data processing.

Figure 2 presents a comparison, over four wave periods, between reconstructed wave elevation and direct LiDAR measurements. We can see that in this shallow water regime the linear reconstruction ( 12) brings little improvement in comparison with the hydrostatic formula [START_REF] Esteva | Comparison of pressure and staff wave gage records[END_REF] and strongly underestimates the crest elevation. A wave by wave analysis shows that the root mean square (RMS) error for the wave height estimated from the linear reconstruction reaches 28% (see also [START_REF] Martins | On the influence of swash-based reflection on surf zone hydrodynamics: a wave-by-wave approach[END_REF]). By contrast, the nonlinear formula (13) gives much better results. In particular, we can see in Figure 3 that the peaked and skewed shape of the wave is much better described by the nonlinear formulation Eq. ( 13) than by the linear one Eq. [START_REF] Guza | Local and shoaled comparisons of sea surface elevations, pressures, and velocities[END_REF]. In order to quantify the wave asymmetry ζH, eq. ( 11) ζSL, eq. ( 12) ζSNL, eq. ( 13) ζLiDAR with respect to the horizontal axis we have computed the skewness parameter:

S k 0.
S k = (ζ -ζ ) 3 (ζ -ζ ) 2 3/2 ,
where . is the time-averaging operator. Table 1 shows that the linear reconstruction strongly underestimates, by 22.5%, the wave skewness. By contrast the skewness is very well reproduced by the nonlinear reconstruction with an error, for this specific test case, smaller than 1%. The nonlinear formula also improves the evaluation of the wave crest elevation (RMS error of 17%), but less significantly compared with the skewness.

The skewness parameter is related to the velocity skewness at the bed (see e.g. Rocha et al. [START_REF] Rocha | Improving the parameterization of wave nonlinearities -The importance of wave steepness, spectral bandwidth and beach slope[END_REF]) which, along with the wave asymmetry, can be used to compute the sediment transport rate (e.g. Abreu et al. [START_REF] Abreu | Bed shear stress under skewed and asymmetric oscillatory flows[END_REF]). The wave high order harmonics also contribute to the acceleration effects that can promote bed motion (e.g. Berni et al. [START_REF] Berni | Effects of horizontal pressure gradients on bed destabilization under waves[END_REF]). This underlines the need to accurately characterize wave non-linearities for predicting sediment transport and morphodynamics (e.g. Dubarbier et al. [START_REF] Dubarbier | Processbased modeling of cross-shore sandbar behavior[END_REF]).

Bichromatic waves

In this section, the ability of our shallow water formulas to reconstruct nonlinear waves is assessed with respect to a bichromatic wave field propagating over a gently sloping (1/20) movable bed. The small-scale experimental set-up is described in Michallet et al. [START_REF] Michallet | Bed motion under waves: plug and sheet flow observations[END_REF] 128 Hz) in the shoaling zone, just prior to wave breaking. The surface elevation was monitored with a capacitive wave gauge. Its accuracy was estimated to be of approximately 1 mm from comparison with high-speed video recordings. The pressure sensor was located at the bed. The still water depth at the measurement location was h 0 = 0.185 m, the characteristic bottom slope was σ b = 0.08 and the values of the dimensionless parameters were μ = 0.28 and, for the highest wave of the group, ε = 0.37.

Figure 4 presents a comparison between shallow water reconstructions and direct elevation measurements. The nonlinear formula (13) significantly improves the elevation reconstruction compared with the linear formula [START_REF] Guza | Local and shoaled comparisons of sea surface elevations, pressures, and velocities[END_REF], especially for the highest waves. The prediction is slightly poorer at the wave crests in the second half of the wave packet. It was observed that the erosion depth and the mobile sediment layer thickness were enhanced after the highest wave. Further investigations would be needed to confirm how the sediment dynamics might contribute to damp high frequencies at the sensor location.

Figure 5 shows a comparison between the measured surface elevation energy density spectrum and the spectra obtained from reconstruction formulas. We can see that the linear [START_REF] Guza | Local and shoaled comparisons of sea surface elevations, pressures, and velocities[END_REF] and nonlinear (13) reconstructions properly describe the elevation energy around the first (i.e. fundamental) harmonic. For higher harmonics, the agreement between the linear reconstruction [START_REF] Guza | Local and shoaled comparisons of sea surface elevations, pressures, and velocities[END_REF] and direct measurements significantly decreases. Figure 5 shows that the linear theory underestimates the energy around the 4 th harmonic by more than one order of magnitude. By contrast, our nonlinear formula (13) accurately predicts the wave spectrum over the four harmonics. This result demonstrates that a nonlinear approach is essential for nonlinear wave reconstruction.

Field data

The data was collected during a field campaign performed on April 13-14 2017, at La Salie beach, situated on the southern part of the French Atlantic coast. Several instruments were deployed at low tide to characterize the shallow water wave field: a stereophotogrammetry system, a Nortek Signature 1000

T m f E (m 2 s) The wave conditions at the water depth h = 2.25 m were characterized by a significant wave height of 0.70 m, with a peak period of 11.1 s (μ = 0.075); large-amplitude wave groups were also observed (see Figure 6). The maximum wave height in these wave groups is 1.4 m, which corresponds to a nonlinearity parameter ε of 0.31. Figure 7 shows that the linear reconstruction [START_REF] Guza | Local and shoaled comparisons of sea surface elevations, pressures, and velocities[END_REF] gives good results for the lowest waves of the wave group but strongly underestimates the elevation at the crest of the highest waves. On the other hand, our nonlinear formula [START_REF] Henderson | Nonlinear generation and loss of infragravity wave energy[END_REF] gives excellent results even for the highest waves. A zoom around the highest wave of the wave group observed in Figure 7 is presented in Figure 8. We can see in this figure that the linear reconstruction fails to predict the maximum elevation and above all the skewed shape of this nonlinear wave.

By contrast, these wave properties are very well reproduced by the nonlinear reconstruction [START_REF] Henderson | Nonlinear generation and loss of infragravity wave energy[END_REF]. The measured surface elevation energy density spectrum and the spectra obtained from reconstruction formulas are compared in Figure 9. In agreement with observations for laboratory bichromatic waves (see Figure 5) we can see that a nonlinear method is required to properly reconstruct the surface elevation spectrum, especially for the highest harmonics.

Conclusion

We have derived a nonlinear weakly dispersive formula (Eq. ( 13)) to recon- the breaking point, and even in the presence of strong wave reflection (Section 3.1), the novel method was found to perform very well.

We have shown that the nonlinear formula (13) provides much better reconstructed surface elevation than linear methods. Our nonlinear reconstruction is able to accurately reproduce the peaked and skewed shape of nonlinear waves prior to breaking. Although the linear reconstruction method properly describes the surface elevation energy density around the fundamental harmonic, it strongly underestimates the energy density for the higher harmonics. In contrast, the nonlinear weakly dispersive reconstruction (13) accurately predicts the surface wave elevation spectrum over four harmonics. This demonstrates that a nonlinear approach is essential for nonlinear wave reconstruction in the shoaling zone and especially near the breaking point. However, the derivation of our nonlinear formula (Eq. ( 13)) is based on a shallowness assumption. For nonlinear waves propagating in intermediate water depth a fully dispersive nonlinear reconstruction is required (see Bonneton and Lannes [START_REF] Bonneton | Recovering water wave elevation from pressure measurements[END_REF]).

An accurate description of surface wave elevation is crucial for many coastal applications, such as wave-induced sediment transport or overtopping and submersion associated with extreme waves. In this context, the application of nonlinear reconstruction methods on pressure data represent an economic and easy to use alternative to direct wave elevation measurement methods, such as acoustic surface tracking or LiDAR scanning that can be hard to deploy in some circumstances. 12) and ( 13) can be computed, at order O(Δt 2 ), as follows

Acknowledgments

∂ t α| n = α n+1 -α n-1 2Δt ∂ 2 t α| n = α n+1 -2α n + α n-1 Δt 2 (A.1) ∂ t α∂ t α | n = α n+1 (α n+2 -α n ) -α n-1 (α n -α n-2 ) 4Δt 2 ,
where α n denotes the discrete value of a flow variable α(t) at t = nΔt and Δt is the measurement time step. This method properly reconstructs the wave 285 field (see figure A.10), except that the measurement noise induces some spurious oscillations (see green line). This drawback is easily overcome (see blue line in figure A.10) by applying a moving average filter to the pressure measurements.

The local method gives similar results than the Fourier one, except at the wave crest where the Fourier method is more accurate.

290

Figure 2 :Figure 3 :

 23 Figure 2: Surface elevation reconstruction of monochromatic waves. A7-mono test obtained during BARDEXII, h 0 = 1.17 m, Tp = 12.1 s and δm = 0.33 m. Dimensionless cut-off frequency Tpfc = 20. black line: direct LiDAR measurement of ζ ; dashed black line: hydrostatic reconstruction ζ H , Eq. (11); green line: ζ SL , Eq. (12); magenta line: ζ SNL , Eq. (13).

Figure 4 : 2 - 1 )

 421 Figure 4: Surface elevation reconstruction of bichromatic waves, f 1 = 0.5515 Hz, f 2 = 0.6250 Hz (Tm = f 1 +f 2 2

Figure 5 : 2 - 1 ), h 0

 5210 Figure 5: Surface elevation energy density spectra, E(f ), as a function of the dimensionless frequency Tmf , for bichromatic waves, f 1 = 0.5515 Hz, f 2 = 0.6250 Hz (Tm = f 1 +f 2 2

Figure 8 :Figure 9 :

 89 Figure 8: Reconstruction of the highest wave observed in a wave group. Cut-off frequency fc = 1 Hz, h = 2.25 m, δm = 0.69 m. dot: direct acoustic measurement of h ; dashed black line: hydrostatic reconstruction ζ H , Eq. (11); green line: h SL = h + ζ SL , Eq. (12); magenta line: h SNL = h + ζ SNL , Eq. (13).

Philippe

  Bonneton and David Lannes have been partially funded by the ANR-13-BS01-0009-01 BOND. David Lannes also acknowledges support from the ANR-13-BS01-0003-01 DYFICOLTI, Philippe Bonneton from the ANR-14-ASTR-0019 COASTVAR and Kévin Martins from the University Research Services scholarship from University of Bath. The BARDEXII data set used in this study is available from http://www.hydralab.info, and it was obtained with the support of the European Community's 7th Framework Programme through the grant to the budget of the Integrating Activity HYDRALAB IV, contract no. 261520. The authors are grateful to all those who contributed to the field experiment at La Salie Beach: Bonneton N., Castelle B., Detandt G., Dubarbier B., Marieu V. and Poncet P.A..Appendix A. Local time discretization method for the surface reconstructionIn contrast to fully dispersive reconstructions, the weakly dispersive methods derived in the present paper are local in time and thus do not necessarily require the use of Fourier transforms. However, to accurately compute time derivatives and to easily filter out the measurement noise we have used Fourier transforms in this paper (see section 3). When the measured pressure time series is too short for a Fourier analysis a local time discretization method is required. The time derivatives involved in equations (

Table 1 :

 1 Sea surface skewness. Comparison between reconstructed elevation and direct LiDAR elevation measurements.

		98	1.12	1.46	1.45
	S k error	32.2 %	22.5 %	0.7 %

  SL = h + ζ SL , Eq. (12); magenta line: h SNL = h + ζ SNL , Eq. (13).

		3.2 3.2 3.2
		3 3 3
		2.8 2.8 2.8
	h (m) h (m) h (m)	2.4 2.6 2.4 2.6 2.4 2.6
		2.2 2.2 2.2
		2 2 2
		1.8 1.8 1.8
		0 200	100 Figure 6: Reconstruction of water depth time series of waves observed in the field. Cut-off 200 300 400 500 t(s) frequency fc = 1 Hz, h = 2.25 m, δm = 0.69 m. dot: direct acoustic measurement of h; magenta line: h SNL = h + ζ SNL , Eq. (13). 210 220 230 240 250 260 t(s) Figure 7: Reconstruction of water depth time series of a group of waves observed in the field. Cut-off frequency fc = 1 Hz, h = 2.25 m, δm = 0.69 m. dot: direct acoustic measurement of 250 252 254 256 h ; green line: h 248 t(s)	600 270 258

  struct, from pressure measurements, the surface elevation of nonlinear waves propagating in shallow water. This simple and easy to use formula is local in time, only involving first and second order time derivatives of the measured pressure P m , or equivalently ζ H . Contrary to fully dispersive reconstruction formulas it does not necessarily require Fourier transform. The ability of this formula to reconstruct surface elevation has been assessed from laboratory and field experiments of weakly dispersive waves (μ ≤ 0.28) propagating in the shoaling zone just prior to breaking. Despite the strong nonlinearities naturally found close to


Yves Coudière
email: yves.coudiere@inria.fr

Charlie Douanla-Lontsi
email: charlie.douanla-lontsi@inria.fr

Charles Pierre
email: charles.pierre@univ-pau.fr

Exponential Adams-Bashforth integrators for stiff ODEs, application to cardiac electrophysiology

Keywords: stiff equations, explicit high-order multistep methods, exponential integrators, stability and convergence, Dahlquist stability

published or not. The documents may come     Exponential Adams Bashforth integrators for stiff ODEs, application to cardiac electrophysiology

Introduction

Computations in cardiac electrophysiology have to face two constraints. Firstly the stiffness due to heterogeneous time and space scales. This is usually dealt with by considering very fine grids. That strategy is associated with large computational costs, still challenging in dimension 3. Secondly, the resolution of the reaction terms that appears in the ionic models has an important cost. That resolution occur at each grid node. The total amount of evaluation of the reaction terms has to be maintained as low as possible by the considered numerical method. Implicit solvers therefore are usually avoided. Exponential integrators are well adapted to cope with these two constraints. Actually they both allow an explicit resolution of the reaction term and display strong stability properties. In this article is studied and analyzed exponential timestepping methods dedicated to the resolution of reaction equations.

Models for the propagation of the cardiac action potential are evolution reaction diffusion equations coupled with ODE systems. The widely used monodomain model [START_REF] Clements | Activation dynamics in anisotropic cardiac tissue via decoupling[END_REF][START_REF] Colli-Franzone | Monodomain simulations of excitation and recovery in cardiac blocks with intramural heterogeneity[END_REF][START_REF] Colli-Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models[END_REF] formulates as,

∂v ∂t = Av + f 1 (v, w, x, t), ∂w ∂t = f 2 (v, w, x, t),
with space and time variables x ∈ Ω ⊂ R d and t ∈ R. The unknowns are the functions v(t, x) ∈ R (the transmembrane voltage) and w(t, x) ∈ R N (a vector that gathers variables describing pointwise the electrophysiological state of the heart cells). In the monodomain model, the diffusion operator is A(:= div(g(x)∇•)), and the reaction terms are the nonlinear functions f 1 , f 2 . These functions model the cellular electrophysiology. They are called ionic models. Ionic models are of complex nature, see e.g. [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF][START_REF] Luo | A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes[END_REF][START_REF] Tusscher | A model for human ventricular tissue[END_REF][START_REF] Iyer | A computational model of the human left-ventricular epicardial myocyte[END_REF]. A special attention has to be paid to the number of evaluations of the functions f 1 and f 2 , and implicit solvers are usually avoided. Though we ultimately use an implicit/explicit method to solve the PDE, we need an efficient, fast and robust method to integrate the reaction terms. Therefore, this article focuses on the time integration of the stiff ODE system dy dt = f (t, y), y(0) = y 0 ,

in the special cases where f (t, y) is an ionic model from cellular electrophysiology. For that case, stiffness is due to the co existence of fast and slow variables. Fast variables are given in [START_REF] Auzinger | Convergence of rational multistep methods of Adams-Padé type[END_REF] by equations of the form,

dy i dt = f i (t, y) = a i (t, y)y i + b i (t, y). (2)
Here a i (t, y) ∈ R is provided by the model. This scalar rate of variation will be inserted in the numerical method to stabilize its resolution. Exponential integrators are a class of explicit methods meanwhile exhibiting strong stability properties. They have motivated many studies along the  Exponential Adams-Bashforth integrators for stiff ODEs past 15 years, among which we quote e.g. [START_REF] Hochbruck | Exponential integrators for large systems of differential equations[END_REF][START_REF] Cox | Exponential time differencing for stiff systems[END_REF][START_REF] Hochbruck | Explicit exponential Runge-Kutta methods for semilinear parabolic problems[END_REF][START_REF] Hochbruck | Exponential Rosenbrock-type methods[END_REF][START_REF] Tokman | New adaptive exponential propagation iterative methods of Runge-Kutta type[END_REF][START_REF] Luan | Explicit exponential Runge-Kutta methods of high order for parabolic problems[END_REF] and refer to [START_REF] Minchev | A review of exponential integrators for first order semi-linear problems[END_REF][START_REF] Hochbruck | Exponential integrators[END_REF][START_REF] Hochbruck | A short course on exponential integrators[END_REF] for general reviews. They already have been used in cardiac electrophysiology, as e.g. in [START_REF] Perego | An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations[END_REF][START_REF] Börgers | Exponential time differencing for Hodgkin-Huxley-like ODEs[END_REF]. Exponential integrators are based on a reformulation of (1) as, dy dt = a(t, y)y + b(t, y), y(0) = y 0 ,

(with f = ay + b) where the linear part a(t, y) is used to stabilize the resolution. Basically a(t, y) is assumed to capture the stiffest modes of the Jacobian matrix of system [START_REF] Auzinger | Convergence of rational multistep methods of Adams-Padé type[END_REF]. Stabilization is brought by performing an exact integration of these modes. This exact integration involves the computation of the exponential exp(a(t n , y n )h) at the considered point. This computation is the supplementary cost for exponential integrators as compared to other time stepping methods. Exponential integrators of Adams type are explicit multistep exponential integrators. They were first introduced by Certaine [START_REF] Certaine | The solution of ordinary differential equations with large time constants[END_REF] in 1960 and Nørsettin [START_REF] Norsett | An A-stable modification of the Adams-Bashforth methods[END_REF] in 1969 for a constant linear part A = a(t, y) in [START_REF] Börgers | Exponential time differencing for Hodgkin-Huxley-like ODEs[END_REF]. The schemes are derived using a polynomial interpolation of the non linear term b(t, y). It recently received an increasing interest [START_REF] Tokman | Efficient integration of large stiff systems of ODEs with exponential propagation iterative (EPI) methods[END_REF][START_REF] Calvo | A class of explicit multistep exponential integrators for semilinear problems[END_REF][START_REF] Ostermann | A class of explicit exponential general linear methods[END_REF] and various convergence analysis have been completed in this particular case [START_REF] Hochbruck | Exponential multistep methods of Adams-type[END_REF][START_REF] Auzinger | Convergence of rational multistep methods of Adams-Padé type[END_REF][START_REF] Koskela | Exponential Taylor methods: analysis and implementation[END_REF]. Non constant linear parts have been less studied. Lee and Preiser [START_REF] Lee | A class of non linear multistep A-stable numerical methods for solving stiff differential equations[END_REF] in 1978 and by Chu [START_REF] Chu | An automatic multistep method for solving stiff initial value problems[END_REF] in 1983 first suggested to rewrite the equation [START_REF] Auzinger | Convergence of rational multistep methods of Adams-Padé type[END_REF] at each time instant t n as, rewritten as, dy dt = a n y + g n (t, y), y(t n ) = y n ,

with a n = a(t n , y n ) and g n (t, y) = b(t, y) + (a(t, y) -a n )y. In the sequel, a n is referred to as the stabilizer. It is updated after each time step. Recently, Ostermann et al [START_REF] Hochbruck | Exponential multistep methods of Adams-type[END_REF][START_REF] Koskela | Exponential Taylor methods: analysis and implementation[END_REF] analyzed the linearized exponential Adams method, where the stabilizer a n is set to the Jacobian matrix of f (t, y) in (1). This choice requires the computation of a matrix exponential at every time step. Moreover, when the fast variables of the system are known, stabilization can be performed only on these variables. Considering the full Jacobian as the stabilizer implies unnecessary computational efforts. To avoid these problems, an alternative is to set the stabilizer as a part or as an approximation of the Jacobian. This has been analyzed in [START_REF] Tranquilli | Rosenbrock-Krylov methods for large systems of differential equations[END_REF] and in [START_REF] Rainwater | A new class of split exponential propagation iterative methods of Runge-Kutta type (sEPIRK) for semilinear systems of ODEs[END_REF] for exponential Rosenbrock and exponential Runge Kutta type methods respectively. That strategy is well adapted to cardiac electrophysiology, where a diagonal stabilizer associated with the fast variable is directly provided by the model with equation [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF]. The present contribution is to analyze general varying a(t, y) in (3) for exponential integrators of Adams type, referred to as exponential Adams-Bashforth, and shortly denoted EAB. Together with the EAB scheme, we introduce a new variant, that we called integral exponential Adams-Bashforth, denoted I-EAB. The convergence analysis held in [START_REF] Hochbruck | Exponential multistep methods of Adams-type[END_REF] extends to the case of general varying stabilizers. However there is a lack of results concerning the stability in that  Yves Coudière, Charlie Douanla-Lontsi and Charles Pierre case: for instance, consider the simpler exponential Euler method, defined by

y n+1 = s(t n , y n , h) := e anh y n + hϕ 1 (a n h)b n , ϕ 1 (z) = (e z -1)/z.
Stability under perturbation (also called 0-stability) can be easily proven provided that the scheme generator s(t, y, h) is globally Lipschitz in y with a constant bounded by 1 + Ch. Therefore stability under perturbation classically is studied by analyzing the partial derivative ∂ y s. This can be done in the case where a(t, y) is either a constant operator or a diagonal varying matrix. In the general case however things turn out to be more complicated. The reason is that we do not have the following general series expansion, as → 0,

e M + N = e M + e M N + O( 2 ),
unless the two matrices M and N are commuting. As a consequence differentiating e a(t,y)h in y cannot be done without very restrictive assumptions on a(t, y). We present here a stability analysis for general varying stabilizers. This will be done by introducing relaxed stability conditions on the scheme generator s(t, y, h). Together with a consistency analysis, it provides a new proof for the convergence of the EAB schemes, in the spirit of [START_REF] Hochbruck | Exponential multistep methods of Adams-type[END_REF].

Stability under perturbation provides results of qualitative nature. In addition, the Dahlquist stability analysis strengthens these results. It is a practical tool that allows to dimension the time step h with respect to the variations of f (t, y) in equation [START_REF] Auzinger | Convergence of rational multistep methods of Adams-Padé type[END_REF]. The analysis is made by setting f (t, y) = λy in [START_REF] Auzinger | Convergence of rational multistep methods of Adams-Padé type[END_REF]. For exponential integrators with general varying stabilizer, the analysis must incorporate the decomposition of f (t, y) = λy used in [START_REF] Börgers | Exponential time differencing for Hodgkin-Huxley-like ODEs[END_REF]. The stability domain of the considered method will depend on the relationship between λ and a(t, y), following a concept first introduced in [START_REF] Perego | An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations[END_REF]. We numerically establish that EAB methods are A(α) stable provided that the stabilizer is sufficiently close to the system Jacobian matrix (precise definitions are in section 4). Moreover the angle α approaches π/2 when the stabilizer goes to the system Jacobian matrix. In contrast, there exists no A(0) stable explicit linear multistep method (see [15, chapter V.2]. This property is remarkable for explicit methods.

Numerical experiments for the EAB and I-EAB scheme are provided in section 6, in the context of cardiac electrophysiology. Robustness to stiffness is studied with this choice. It is numerically shown to be comparable to implicit methods both in terms of accuracy and of stability condition on the time step. We conclude that EAB methods are well suited for solving stiff differential problems. In particular they allow computations at large time step with good accuracy properties and cheap cost.

The article is organized as follows. The EAB and I-EAB methods are introduced in section 1. The general stability and convergence results are stated and proved in section 2. The EAB and I-EAB stability under perturbation and convergence are proved in section 3. The Dahlquist stability is investigated in section 4, and the numerical experiments end the article, in section 6.

 Exponential Adams-Bashforth integrators for stiff ODEs

In all this paper, h > 0 is a constant time-step and t n = nh are the time instants associated with the numerical approximate y n of the solution of the ODE (1).

1 Scheme definitions 1.1 The EAB k method

The exact solution at time t n+1 to the equation (4) (with a n = a(t n , y n )) is given by the variation of the constants formula:

y(t n+1 ) = e anh y(t n ) + h 0 e -anτ g n (t n + τ, y(t n + τ )) dτ . ( 5
)
Using the k approximations y n-j y(t n-j ) for j = 0 . . . k -1, we build the Lagrange polynomial gn of degree at most k -1 that satisfies, gn (t n-j ) = g nj := g(t n-j , y n-j ), 0

≤ j ≤ k -1. ( 6
)
It provides the numerical approximation y n+1 y(t n+1 ) as

y n+1 = e anh y n + h 0 e -anτ gn (t n + τ ) dτ . (7)
The Taylor expansion of the polynomial gn is gn (t n +τ ) = k j=1 γ nj (j-1)! (τ /h) j-1 , where the coefficients γ nj are uniquely determined by [START_REF] Chu | An automatic multistep method for solving stiff initial value problems[END_REF], and actually given in table 1 for k = 1, 2, 3, 4. An exact integration of the integral in equation ( 7) may be performed:

y n+1 = e anh y n + h k j=1 ϕ j (a n h)γ nj , (8)
where the functions ϕ j , originally introduced in [START_REF] Norsett | An A-stable modification of the Adams-Bashforth methods[END_REF], are recursively defined (for j ≥ 0) by

ϕ 0 (z) = e z , ϕ j+1 (z) = ϕ j (z) -ϕ j (0) z and ϕ j (0) = 1 j! . (9)
The equation ( 8) defines the Exponential Adams-Bashforth method of order k, denoted by EAB k .

Remark 1. When a(t, y) = diag(d i ) is a diagonal matrix, ϕ k (a n h) = diag(ϕ k (d i ))
can be computed component-wise. Its computation is straightforward.

Remark 2. With the definition (9), the functions ϕ k are analytic on the whole complex plane. Therefore the EAB k scheme definition (8) makes sense for a matrix term a(t, y) in equation (3) without particular assumption.

Remark 3. The computation of y n+1 in the formula (8) requires the computation of ϕ j (a n h) for j = 0, . . . , k. This computational effort can be reduced with the recursive definition [START_REF] Colli-Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models[END_REF]. In practice only ϕ 0 (a n h) needs to be computed. This is detailed in section 6.1.



Yves Coudière, Charlie Douanla-Lontsi and Charles Pierre g n -g n-1

3 2 g n -2g n-1 + 1 2 g n-2 11 6 g n -3g n-1 + 3 2 g n-2 -1 3 g n-3 γ n3 g n -2g n-1 + g n-2 2g n -5g n-1 + 4g n-2 -g n-3 γ n4 g n -3g n-1 + 3g n-2 -g n-3
1.2 A variant: the I-EAB k method

If the matrix a(t, y) is diagonal, we can take advantage of the following version for the variation of the constants formula:

y(t n+1 ) = e An(h) y(t n ) + h 0 e -An(τ ) b(y(t n + τ ), t n + τ ) dτ ,
where A n (τ ) = τ 0 a(t n + σ, y(t n + σ))dσ. An attempt to improve the EAB k formula ( 8) is to replace a(t, y) and b(t, y) in the integral above by their Lagrange interpolation polynomials. At time t n , we define the two polynomials ãn and bn of degree at most k -1 so that: ãn (t n-j ) = a(t n-j , y n-j ), bn (t n-j ) = b(t n-j , y n-j ), j = 0 . . . k -1, and the primitive Ãn (τ ) = τ 0 ã(t n +σ, y(t n +σ))dσ. The resulting approximate solution at time t n+1 is finally given by the formula

y n+1 = e Ãn(h) y n + h 0 e -Ãn(τ ) bn (t n + τ ) dτ . ( 10
)
The method is denoted I-EAB k (for Integral EAB k ). Dislike for the formula [START_REF] Clements | Activation dynamics in anisotropic cardiac tissue via decoupling[END_REF], no exact integration formula is available, because of the term e -Ãn(τ ) . A quadrature rule is required for the actual numerical computation of the integral in formula [START_REF] Coudiere | Rush Larsen time stepping methods of high order for stiff problems in cardiac electrophysiology[END_REF]. Implementation details are given in section 6.1.

Stability conditions and convergence

The equation ( 1) is considered on a finite dimensional vector space E with norm | • | E . We fix a final time T > 0 and assume that equation (1) has a solution y on [0, T ]. We adopt the general settings for the analysis of k-multistep methods following [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF]. The space E k is equipped with the maximum norm

|Y | ∞ = max 1≤i≤k |y i | E with Y = (y 1 , . . . , y k ) ∈ E k . A k-multistep scheme is defined by a mapping, s : (t, Y, h) ∈ R × E k × R + → s(t, Y, h) ∈ E.
 Exponential Adams-Bashforth integrators for stiff ODEs For instance, the EAB k scheme rewrites as y n+1 = s(t n , Y, h) with Y = (y n-k+1 , . . . , y n ) in the formula [START_REF] Colli-Franzone | Monodomain simulations of excitation and recovery in cardiac blocks with intramural heterogeneity[END_REF]. The scheme generator is the mapping S given by

S : (t, Y, h) ∈ R × E k × R + -→ (y 2 , . . . , y k , s(t, Y, h)) ∈ E k . A numerical solution is a sequence (Y n ) in E k for n ≥ k -1 so that Y n+1 = S(t n , Y n , h) for n ≥ k -1, (11)
Y k-1 = (y 0 , . . . , y k-1 ) being its initial condition. A perturbed numerical solu- tion is a sequence (Z n ) in E k for n ≥ k -1 such that Z k-1 = Y k-1 + ξ k-1 , Z n+1 = S(t n , Z n , h) + ξ n+1 for n ≥ k -1, (12)
with

(ξ n ) ∈ E k for n ≥ k -1.
The scheme is said to be stable under perturbation (or 0-stable) if: being given a numerical solution (Y n ) as in [START_REF] Cox | Exponential time differencing for stiff systems[END_REF] there exists a (stability) constant L s > 0 so that for any perturbation (Z n ) as defined in [START_REF] Douanla-Lontsi | Efficient high order schemes for stiff ODEs in cardiac electrophysiology[END_REF] we have,

max k-1≤n≤T /h |Y n -Z n | ∞ ≤ L s k-1≤n≤T /h |ξ n | ∞ . (13)
Proposition 1. Assume that there exists constants C 1 > 0 and C 2 > 0 such that

1 + |S(t, Y, h)| ∞ ≤ (1 + |Y | ∞ ) (1 + C 1 h), ( 14
) |S(t, Y, h) -S(t, Z, h)| ∞ ≤ |Y -Z| ∞ (1 + C 2 h(1 + |Y | ∞ )) , (15)
for 0 ≤ t ≤ T , and for Y, Z ∈ E k . Then, the numerical scheme is stable under perturbation with the constant L s in (13) given by

L s = e C T , C := C 2 e C 1 T (1 + |Y k-1 | ∞ ) . (16)
Proof. Consider a numerical solution (Y n ) in [START_REF] Cox | Exponential time differencing for stiff systems[END_REF]. A recursion on condition [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF] gives,

1 + |Y n | ∞ ≤ (1 + |Y k-1 | ∞ ) (1 + C 1 h) n-k+1 ≤ e C 1 T (1 + |Y k-1 | ∞ ) , since (1 + x) p ≤ e px and (n -k + 1)h ≤ nh ≤ T . Now, consider a perturbation (Z n ) of (Y n )
given by [START_REF] Douanla-Lontsi | Efficient high order schemes for stiff ODEs in cardiac electrophysiology[END_REF]. Using the condition [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF] together with the previous inequality,

|Y n+1 -Z n+1 | ∞ ≤ |S(t n , Y n , h) -S(t n , Z n , h)| ∞ + |ξ n+1 | ∞ ≤ |Y n -Z n | ∞ (1 + C 2 h(1 + |Y n | ∞ )) + |ξ n+1 | ∞ ≤ |Y n -Z n | ∞ 1 + C 2 e C 1 T (1 + |Y k-1 | ∞ )h + |ξ n+1 | ∞ ≤ |Y n -Z n | ∞ (1 + C h) + |ξ n+1 | ∞ ,  where C := C 2 e C 1 T (1 + |Y k-1 | ∞ )
. By recursion we get,

|Y n -Z n | ∞ ≤ (1 + C h) n-k+1 |Y k-1 -Z k-1 | ∞ + n-k i=0 (1 + C h) i |ξ n-i | ∞ ≤ (1 + C h) n n i=k-1 |ξ i | ∞ ≤ e c T n i=k-1 |ξ i | ∞ ,
which ends the proof.

Like in the classical cases, stability under perturbation together with consistency ensures convergence. Let us specify this point. For the considered solution y(t) of problem ( 1) on [0, T ], we define,

Y (t) = (y (t -(k -1)h) , . . . y(t)) ∈ E k for (k -1)h ≤ t ≤ T. (17)
The local error at time t n is,

(t n , h) = Y (t n+1 ) -S(t n , Y (t n ), h). ( 18
)
The scheme is said to be consistent of order p if there exists a (consistency) constant L c > 0 only depending on y(t) such that max

k-1≤n≤T /h | (t n , h)| ∞ ≤ L c h p+1 .
Corollary 1. If the scheme satisfies the stability conditions ( 14) and [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF], and is consistent of order p, then a numerical solution (Y n ) given by ( 11) satisfies,

max k-1≤n≤T /h |Y (t n ) -Y n | ∞ ≤ L s L c T h p + L s |ξ 0 | ∞ , (19)
where ξ 0 = Y (t k-1 ) -Y k-1 denotes the error on the initial data, and the constant L s is as in equation ( 16).

Remark 4. Note that the stability constant L s in ( 16) depends on |Y k-1 | ∞ , and then on h. This is not a problem since L s can be bounded uniformly as

h → 0 for Y k-1 in a neighbourhood of y 0 . Proof. We have Y (t k-1 ) = Y k-1 + ξ 0 and Y (t n+1 ) = S(t n , Y (t n ), h) + (t n , h).
Therefore the sequence (Y (t n )) is a perturbation of the numerical solution (Y n ) in the sense of [START_REF] Douanla-Lontsi | Efficient high order schemes for stiff ODEs in cardiac electrophysiology[END_REF]. As a consequence, proposition 1 shows that max

k-1≤n≤T /h |Y n -Y (t n )| E ≤ L s   |ξ 0 | + k≤n≤T /h | (t n , h)|   ≤ L s |ξ 0 | + L s L c   k≤n≤T /h h   h p ,
and the convergence result follows.

 Exponential Adams-Bashforth integrators for stiff ODEs

EAB k and I-EAB k schemes analysis

The space E is assumed to be E = R N with its canonical basis and with | • | E the maximum norm. The space of operators on E is equipped with the associated operator norm, and associated to N × N matrices. Thus a(t, y) is a N × N matrix and its norm |a(t, y)| is the matrix norm associated to the maximum norm on R N . It is commonly assumed for the numerical analysis of ODE solvers that f in the equation ( 1) is uniformly Lipschitz in its second component y. With the formulation (3), the following supplementary hypothesis will be needed: on

R × E, |a(t, y)| ≤ M a , a(t, y), b(t, y) and f (t, y) uniformly Lipschitz in y. ( 20
)
We denote by K f , K a and K b the Lipschitz constant for f , a and b respectively.

Theorem 1. With the assumptions (20), the EAB k and I-EAB k schemes are stable under perturbations. Moreover, if a and b are C k regular on R × E, then the EAB k and I-EAB k schemes are consistent of order k. Therefore they converge with order k in the sense of inequality [START_REF] Hochbruck | Exponential integrators[END_REF], by applying corollary 1.

The stability and consistency are proved in sections 3.3 and 3.4, respectively. Preliminary tools and definitions are provided in the sections 3.1 and 3.2.

Interpolation results

Consider a function x : R × E -→ R and a triplet (t, Y, h) ∈ R × E k × R + with Y = (y 1 , . . . , y k ). We set to x[t,Y,h] the polynomial with degree less than k -1 so that

x[t,Y,h] (t -ih) = x(t -ih, y k-i ), 0 ≤ i ≤ k -1.
We then extend component-wise this definition to vector valued or matrix valued functions x (e.g. the functions a, b or f ). Lemma 1. There exists an (interpolation) constant L i > 0 such that, for any function

x : R × E → R, sup t≤τ ≤t+h |x [t,Y,h] (τ )| ≤ L i max 0≤i≤k-1 |x(t -ih, y k-i )| , (21)
sup

t≤τ ≤t+h |x [t,Y 1 ,h] -x[t,Y 2 ,h] | ≤ L i max 0≤i≤k-1 |x(t -ih, y 1,k-i ) -x(t -ih, y 2,k-i )| , (22)
Consider a function y : [0, T ] → E and assume that x and y have a

C k regularity. Then, when [t -(k -1)h, t + h] ⊂ [0, T ], sup t≤τ ≤t+h x(t, y(t)) -x[t,Y (t),h] E ≤ sup [0,T ] d k f dt k (f (t, y(t))) h k , (23)
 with Y (t) defined in [START_REF] Hochbruck | Exponential integrators for large systems of differential equations[END_REF]. For a vector valued function in R d the previous inequalities hold when considering the max norm on R d . For a matrix valued function in R d × R d this is also true for the operator norm on R d × R d when multiplying the constants in the inequalities ( 21), ( 22) and ( 23) by d.

Proof. The space P k-1 of the polynomials p with degree less than k -1 is equipped with the norm sup

[0,1] |p(τ )|. On R k is considered the max norm. To R = (r 1 , . . . , r k ) ∈ R k is associated LR ∈ P k-1 uniquely determined by LR(-i) = r k-1 for i = 0 . . . k -1.
The mapping L is linear. Let C L be its continuity constant (it only depends on k).

We fix

x : R × E → R and (t, h) ∈ R × R + . Consider the vector Y = (y 1 , . . . , y k ) ∈ E k and set the vector R ∈ R k by R = (x(t -(k -1)h, y 1 ), . . . , x(t, y k )). We have x[t,Y,h] (t + τ ) = LR(τ /h).
The relation [START_REF] Hochbruck | Exponential Rosenbrock-type methods[END_REF] exactly is the continuity of L and

L i = C L . Consider Y 1 , Y 2 ∈ E k and the associated vectors R 1 , R 2 as above. We have (x [t,Y 1 ,h] -x [t,Y 2 ,h] )(t + τ ) = L (R 1 -R 2 ) (τ /h).
Again, relation [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] comes from the continuity of L.

Let ϕ : R → R be a C k function, its interpolation polynomial φ at the points t -(k -1)h, . . . , t is considered. A classical result on Lagrange interpolation applied to ϕ states that, for all τ ∈ (t, t+h), there exists ξ [START_REF] Iyer | A computational model of the human left-ventricular epicardial myocyte[END_REF] by setting ϕ(t) = x(t, y(t)). For a vector valued function x : R × E → R d , these three inequalities holds by processing component-wise and when considering the max norm on R d . For a matrix valued function

∈ (t-(k-1)h, t+h), such that (ϕ -φ) (τ ) = 1 k! ϕ (k) (ξ)π(τ ), where π(τ ) = k i=1 (τ -t i ). For τ ∈ (t, t + h), we have |π(τ )| ≤ k! h k . This proves
x : R × E → R d × R d , the extension is direct when considering the max norm | • | ∞ on R d × R d (i.e.

the max norm on the matrix entries). The operator norm

| • | is retrieved with the inequality | • | ∞ ≤ d| • | 3.2 Scheme generators Let us consider (t, Y, h) ∈ R × E k × R + with Y = (y 1 , . . . , y k ).
With the notations used in the previous subsection, we introduce the interpolations ã[t,Y,h] and b[t,Y,h] for the functions a and b in [START_REF] Börgers | Exponential time differencing for Hodgkin-Huxley-like ODEs[END_REF]. Thanks to the definition [START_REF] Coudiere | Rush Larsen time stepping methods of high order for stiff problems in cardiac electrophysiology[END_REF], the I-EAB k scheme generator is defined by

s(t, Y, h) = z(t + h) with dz dτ = ã[t,Y,h] (τ )z(τ ) + b[t,Y,h] (τ ), z(t) = y k , (24)
We introduce the polynomial ḡ[t,Y,h] with degree less than k -1 that satisfies,

ḡ[t,Y,h] (t -ih) = f (t -ih, y k-i ) -a(t, y k )y k-i , i = 0 . . . k -1.
The function gn in ( 6) is given by gn = ḡ[tn,Yn,h] with Y n = (y n-k+1 . . . , y n ).

With the definition [START_REF] Clements | Activation dynamics in anisotropic cardiac tissue via decoupling[END_REF], the EAB k scheme generator is defined by

s(t, Y, h) = z(t + h) with dz dτ = a(t, y k )z(τ ) + ḡ[t,Y,h] (τ ), z(t) = y k , (25)
 Exponential Adams-Bashforth integrators for stiff ODEs

We will use the fact that ḡ[t,Y,h] is the interpolator of the function g t,y k : (τ, ξ) → f (τ, ξ) -a(t, y k )ξ, precisely:

ḡ[t,Y,h] = g t,y k [t,Y,h] .
These scheme generator definitions will allow us to use the following Gronwall's inequality (see [START_REF] Dragomir | Some Gronwall type inequalities and applications[END_REF]Lemma 196,p.150]).

Lemma 2. Suppose that z(t) is a C 1 function on E. If there exist α > 0 and β > 0 such that |z (t)| E ≤ α|t| + β for all t ∈ [t 0 , t 0 + h], then:

|z(t)| E ≤ |z(t 0 )| E e αh +βh e αh for t ∈ [t 0 , t 0 + h]. ( 26
)

Stability

With proposition 1 we have to prove the stability conditions ( 14) and [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF]. Note that it is sufficient to prove these relations for h ≤ h 0 for some constant h 0 > 0 since the limit h → 0 is of interest here.

Case of the I-EAB k scheme

Consider (t, h) ∈ R × R + and a vector Y = (y 1 , . . . , y k ) ∈ E k . We simply denote ã = ã[t,Y,h] and b = b[t,Y,h] .
The scheme generator is given by [START_REF] Koskela | Exponential Taylor methods: analysis and implementation[END_REF]. We first have to bound z(t + h) where z is given by z = ãz + b, z(t) = y k .

On the first hand, with the interpolation bound [START_REF] Hochbruck | Exponential Rosenbrock-type methods[END_REF],

sup t≤τ ≤t+h |ã(τ )| ≤ L i max 0≤i≤k-1 |a(t -ih, y k-1 )| ≤ L i M a .
On the second hand, the function b(t, y) is globally Lipschitz in y and thus can be bounded as

|b(t, y)| E ≤ |b(t, 0)| E + K b |y| E ≤ R b (|y| E + 1
), for 0 ≤ t ≤ T and for some constant R b only depending on K b and on T . Then with the bound [START_REF] Hochbruck | Exponential Rosenbrock-type methods[END_REF],

sup t≤τ ≤t+h | b(τ )| E ≤ L i max 0≤i≤k-1 R b (|y k-i | E + 1) ≤ L i R b (|Y | ∞ + 1).
By applying the Gronwall inequality [START_REF] Luan | Explicit exponential Runge-Kutta methods of high order for parabolic problems[END_REF], for 0 ≤ τ ≤ h,

|z(t + τ )| E ≤ e L i Mah (|y k | E + hL i R b (|Y | ∞ + 1)) .
Thus, there exists a constant C 1 only depending on L i , M a , R b and h 0 so that, for 0 ≤ τ ≤ h and for 0 ≤ h ≤ h 0 ,

|z(t + τ )| E ≤ C 1 h + |Y | ∞ (1 + C 1 h). (27)
This gives the condition ( 14) taking τ = h.

For j=1, 2 We consider Y j = (y j,1 , . . . y j,k ) ∈ E k and denote ãj = ã[t,Y j ,h] and bj = b[t,Y j ,h] the interpolations of the functions a and b. With the definition [START_REF] Koskela | Exponential Taylor methods: analysis and implementation[END_REF] of the I-EAB k scheme we have:

|s(t, Y 1 , h) -s(t, Y 2 , h)| E = |δ(t + h)| with δ = z 1 -z 2 and
with z j given by z j = ãj z j + bj , z j (t) = y j,k .

We then have,

δ = ã1 δ + r, r := (ã 1 -ã2 )z 2 + ( b1 -b2 ).
Using that a and b are Lipschitz in y and with the interpolation bound [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF],

sup t≤τ ≤t+h | b1 (τ ) -b2 (τ )| E ≤ L i K b |Y 1 -Y 2 | ∞ , sup t≤τ ≤t+h |ã 1 (τ ) -ã2 (τ )| ≤ L i K a |Y 1 -Y 2 | ∞ .
With the upper bound [START_REF] Luo | A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes[END_REF], for t ≤ τ ≤ t + h ≤ T and for h ≤ h 0 ,

|r(τ )| E ≤ L i |Y 1 -Y 2 | ∞ (K b + K a (C 1 h + |Y 2 | ∞ (1 + C 1 h))) ≤ C|Y 1 -Y 2 | ∞ (1 + |Y 2 | ∞ ) , (28)
For a constant C only depending on h 0 , K a , K b , L i and C 1 . We finally apply the Gronwall inequality [START_REF] Luan | Explicit exponential Runge-Kutta methods of high order for parabolic problems[END_REF]. It yields,

|δ(t + h)| ≤ e L i Mah (|y 1,k -y 2,k | E + Ch|Y 1 -Y 2 | ∞ (1 + |Y 2 | ∞ )) ≤ |Y 1 -Y 2 | ∞ e L i Mah (1 + Ch (1 + |Y 2 | ∞ )) ,
This implies the second stability condition [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF] for h ≤ h 0 .

Case of the EAB k scheme

Consider (t, h) ∈ R × R + and a vector Y = (y 1 , . . . , y k ) ∈ E k . Following the definition of the EAB k scheme given in section 1.1, we denote ā = a(t, y k ), g the function g(τ, ξ) = f (τ, ξ) -āξ and ḡ = ḡ[t,Y,h] . We have that ḡ is the Lagrange interpolation polynomial of g, specifically ḡ = g[t,Y,h] . The scheme generator is then given by equation [START_REF] Lee | A class of non linear multistep A-stable numerical methods for solving stiff differential equations[END_REF]:

s(t, Y, h) = z(t + h) with, z = āz + ḡ, z(t) = y k .
We first have the bound |ā| ≤ M a . As in the previous subsection, f being globally Lipschitz in y, one can find a constant R f so that for 0

≤ t ≤ T , |f (t, y)| E ≤ R f (1 + |y| E ). It follows that |g(τ, ξ)| E ≤ R f (|y| E + 1) + M a |y| E ≤ C 0 /L i (|y| E + 1), with C 0 /L i = R f + M a .
Therefore, with the interpolation bound [START_REF] Hochbruck | Exponential Rosenbrock-type methods[END_REF]:

sup t≤τ ≤t+h |ḡ(τ )| E ≤ L i max 0≤i≤k-1 |g(t -ih, y k-i )| E ≤ C 0 (|y| E + 1).

 Exponential Adams-Bashforth integrators for stiff ODEs

By applying the Gronwall inequality [START_REF] Luan | Explicit exponential Runge-Kutta methods of high order for parabolic problems[END_REF], for 0 ≤ τ ≤ h,

|z(t + τ )| E ≤ e Mah (|y k | E + hC(|Y | ∞ + 1)) .
Thus, there exists a constant C 1 only depending on M a and C 0 so that, for 0 ≤ τ ≤ h and for 0 ≤ h ≤ h 0 , [START_REF] Luo | A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes[END_REF] holds. This gives the condition [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF].

We now consider Y 1 , Y 2 ∈ E k for j =1, 2 and denote as previously, āj = a(t, y j,k ), g j the function g j (τ, ξ) = f (τ, ξ) -āj ξ and ḡj = ḡ[t,Y j ,h] . With [START_REF] Lee | A class of non linear multistep A-stable numerical methods for solving stiff differential equations[END_REF],

|s(t, Y 1 , h) -s(t, Y 2 , h)| E = |δ(t + h)| with δ = z 1 -z 2 and
with z j given by z j = āj z j + ḡj , z j (t) = y j,k .

The function δ satisfies the ODE

δ = ā1 δ + r(t), r(t) := (ā 1 -ā2 )z 2 + (ḡ 1 -ḡ2 ).
We have

|g 1 (τ, y 1,i ) -g 2 (τ, y 2,i )| E ≤ |f (τ, y 1,i ) -f (τ, y 2,i )| E + |ā 1 ||y 1,i -y 2,i | E + |ā 1 -ā2 ||y 2,i | E ≤ |Y 1 -Y 2 | ∞ (K f + M a + K a |Y 2 | ∞ ) .
Thus, with equation ( 22), for some C > 0,

sup t≤τ ≤t+h |ḡ 1 (τ ) -ḡ2 (τ )| E ≤ L i max 0≤i≤k-1 |g 1 (t -ih, y 1,k-i ) -g 2 (t -ih, y 2,k-i )| E ≤ C|Y 1 -Y 2 | ∞ (1 + |Y 2 | ∞ ) .
Meanwhile we have the upper bound (27) that gives, for t ≤ τ ≤ t + h ≤ T and h ≤ h 0 ,

|(ā 1 -ā2 )z 2 | E ≤ M a |Y 1 -Y 2 | ∞ |z 2 (τ )| E ≤ M a |Y 1 -Y 2 | ∞ (C 1 h + |Y 2 | ∞ (1 + C 1 h)) .
Altogether, we retrieve the upper bound (28) on r(t). We can end the proof as for the I-EAB k case and conclude that the stability condition (15) holds for the EAB k scheme.

Consistency

A solution y(t) to problem (1) on [0, T ] is fixed. The functions a and b in (3) are assumed to be k regular so that y is C k+1 regular.

Case of the EAB k scheme

The local error [START_REF] Hochbruck | Explicit exponential Runge-Kutta methods for semilinear parabolic problems[END_REF] for the EAB k scheme has been analyzed in [START_REF] Hochbruck | Exponential multistep methods of Adams-type[END_REF]. That analysis remains valid for the case presented here and we only briefly recall it. The local error is obtained by subtracting (7) to (5):

| (t n , h)| E ≤ h 0 e Ma(h-τ ) |g n (t + τ, y(t + τ )) -gn (t + τ )| E dτ ≤ hϕ 1 (M a h)h k sup [0,T ] d k dt k (g n (t, y(t))) ,
thanks to the interpolation error estimate [START_REF] Iyer | A computational model of the human left-ventricular epicardial myocyte[END_REF]. Finally, with the upper bound M a on a n , the last term can be bounded independently on n for h ≤ h 0 . The following constants only depend on the considered exact solution y, on the functions a and b in problem (3) and on T ,

C y = sup [0,T ] |y| E , C a,y = sup [0,T ] d k dt k a(t, y(t)) , C b,y = sup [0,T ] d k dt k b(t, y(t)) .
With the interpolation bound [START_REF] Iyer | A computational model of the human left-ventricular epicardial myocyte[END_REF] Thus the EAB k scheme is consistent of order k.

4 Dahlquist stability

Background

The classical framework for the Dahlquist stability analysis is to set f (t, y) in problem (1) to f = λy. For linear multistep methods, see e.g. [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF], the numerical solutions satisfy |y n+1 /y n | ≤ ρ(λh), where ρ : C → R + . The function ρ is the stability function. Its definition is detailed below. The stability domain is defined by D = {z ∈ C, ρ(z) < 1}. The scheme is said to be:

• A stable if C -⊂ D,
• A(α) stable if D contains the cone with axis R -and with half angle α,

• A(0) stable if R -⊂ D,
• stiff stable if D contains a half plane Re z < x ∈ R -. For exponential integrators, when setting a(t, y) = λ in the reformulation (3) of problem (1), the scheme is exact, and therefore also A stable. Such an equality does not hold in general. Then for exponential integrators the Dahlquist stability analysis has to incorporate the relationship between the stabilization term a(t, y) in (3) and the test function f = λy. This is done here by considering the splitting,



f = λy = ay + b, a = θλ and b = λ(1 -θ)y, (29)
The parameter θ > 0 controls with what accuracy the exact linear part of f in equation 1 is captured by a in equation 3. In practice θ = 1, though we may hope that θ -1 is small. In that framework, the stability function and the stability domain depend on θ, following the idea of Perego and Venezziani in [START_REF] Perego | An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations[END_REF] (see the remark below). For a fixed θ, the stability function is ρ θ so that

y n+1 y n ≤ ρ θ (λh),
and the stability domain is

D θ = {z ∈ C, ρ θ (z) < 1}.
Remark 5. In [START_REF] Perego | An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations[END_REF] an alternative splitting was introduced:

f = λy = ay + b, a = r 1 + r λ and b = 1 1 + r λy,
involving the real parameter r. We modified that splitting in order to adress more easily the limit where a = λ that corresponds to r → ±∞ here whereas it corresponds to θ = 1 in our case. Nevertheless the results can be translated from one splitting to the other with the change of variable r = θ/(1 -θ).

The stability function ρ θ (z) is defined as in the classical case for linear multistep methods, see [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF]Ch. V.1]. For a fixed value of θ, the EAB k scheme applied to the right hand side in equation ( 29) reads

y n+1 + c θ,1 (λh)y n + . . . + c θ,k (λh)y n-k+1 = 0.
The coefficients c θ,j are explicit (for the EAB 2 scheme they are given by c θ,1 (z) = -1 -ϕ 1 (θz)z -ϕ 2 (θz)(1 -θ)z and c θ,2 (z) = ϕ 2 (θz)(1 -θ)z). We fix z ∈ C and consider the polynomial ξ k + c θ,1 (z)ξ k-1 + . . . + c θ,k . It has k complex roots ξ 1 θ (z), . . . , ξ k θ (z). The stability function is defined as

ρ θ (z) = max 1≤j≤k |ξ j θ (z)|.
The stability function ρ θ is numerically analyzed. For a given vakue of θ. The software Maple is used to compute the polynomial roots on a given grid. That grid either is a 1D grid of [-30, 0] ⊂ R -with size ∆x = 0.01 to study the A(0)-stability. Or it is a 2D cartesian grid in C to study the stability domain.

In that case the grid size is ∆x = 0.05 and the domain is [-40, 2] × [0, 60] (roughly 900 000 nodes). The isoline ρ θ (z) = 1 then is constructed with the software gnuplot.

A(0) stability

The stability functions ρ θ (z) are numerically studied for z ∈ R -. These functions have been plotted for different values of the parameter θ. The results are depicted on figure 1. A limit lim -∞ |ρ θ | is always observed. The scheme is A(0) stable when this limit is lower than 1. From Figure 1,

• EAB 2 scheme is A(0) stable if θ ≥ 0.75, • EAB 3 scheme is A(0) stable if 0.88 ≤ θ ≤ 1.9,



Exponential Adams-Bashforth integrators for stiff ODEs

• EAB 4 scheme is A(0) stable if 0.94 ≤ θ ≤ 1.2.
Roughly speaking, A(0) stability holds for the EAB k scheme if the exact linear part of f (t, y) in problem ( 1) is approximated with an accuracy of 75 %, 85 % or 95% for k = 2, 3 or 4 respectively. • Figure 2 shows that the EAB 2 scheme is A(α) stable when θ = 0.75, 0.8 and 0.9 with α 50, 60 and 80 angle degrees respectively.

• Figure 3 displays A(α) stability with α 60, 70 and 60 angle degrees for θ = 0.88, 0.9 and 1.9 respectively for the EAB 3 scheme.

• For the EAB 4 scheme eventually, A(α) stability holds with an angle α approximately of 65, 70 and 60 degrees for θ = 0.94, 0.95 and 1.2 respectively, as shown on figure 4.

In all cases, when A(α) stability is observed, the unstable region inside C - is made of a discrete collection of uniformly bounded sets located along the imaginary axes. Hence, the stability domain D θ also contains half planes of the form Re(z) ≤ a < 0. We conjecture that, when θ is so that the EAB k scheme is A(α)-stable, then it is also stiff stable.

Conclusion

For explicit linear multistep methods, A(0) stability cannot occur, see [ It is interesting to compare the results on the EAB 2 scheme with the results in [START_REF] Perego | An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations[END_REF] concerning another explicit exponential scheme of order 2 named AB2*. The AB2* scheme displays stronger absolute stability properties than the EAB 2 scheme. The AB2* scheme is A(0)-stable for θ ≥ 2/3 whereas the EAB 2 scheme requires θ ≥ 0.75. The AB2* scheme is A-stable for θ ≥ 1 whuich property is not achieved by the EAB 2 scheme when θ = 1. Extensions of that scheme to higher orders have also been developped in [START_REF] Coudiere | Rush Larsen time stepping methods of high order for stiff problems in cardiac electrophysiology[END_REF] that however loses the A(0)-stability property.

Positivity properties

Consider the scalar equation y = f (t, y), with f (t, y) chosen so that there exists a function a(t, y) satisfying a(t, y) ≤ 0 and a(t, y)(y

-K 1 ) ≤ f (t, y) ≤ a(t, y)(y -K 2 ) (30)
for two constants

K 1 ≤ K 2 . Then the solution y(t) satisfies K 1 ≤ y(t) ≤ K 2 for t > 0 provided that y(0) ∈ [K 1 , K 2 ].
A proof for that property is given in [START_REF] Perego | An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations[END_REF].

It can be important for the time-stepping method to preserve that positivity property. This is in particular the case for the application in section 6 where the gating variables w(t) satisfy the property [START_REF] Ostermann | A class of explicit exponential general linear methods[END_REF] with K 1 = 0, K 2 = 1 and with a(t, y) defining the stabilizer of the EAB k scheme.

For the AB * 2 scheme introduced in [START_REF] Perego | An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations[END_REF] (that is explicit exponential and with order 2), that property has been shown to be preserved for a constant linear part a(t, y) = a provided that h ≤ log(3)/|a| and adding a condition on the initial data. These results can be extended to the EAB 2 and to the EAB 3 schemes. It provides sufficient condition for the positivity. However the technique of the proof fails to apply to the EAB 4 scheme. For simplicity in the sequel, we simply denote ϕ k = ϕ k (ah). Proposition 2. Consider the equation y = ay + b(t, y) with a ≤ 0 and

-aK 1 ≤ b(t, y) ≤ -aK 2 ,
so that the condition (30) holds for f (t, y) = ay + b(t, y). Denoting y n the numerical solution for the EAB k scheme, we have that y n ∈ [K 1 , K 2 ] for n ≥ k provided that the following conditions are satisfied.

• EAB 2 scheme: h ≤ 1/|a| and

C p K 1 ≤ e ah y 1 -hϕ 2 b 0 ≤ C p K 2 . (31)
with C p = e ah + ahϕ 2 .

• EAB 3 scheme: h ≤ β 3 /|a| for a constant β 3 0.331 and

C p K 1 ≤ e ah y 2 -2h(ϕ 2 + ϕ 3 )b 1 ≤ C p K 2 . (32)
with C p = e ah + 2ah ϕ 2 + ϕ 3 .

The proof of proposition 2 will be based on the following result.

Lemma 3. Consider the equation y = ay +b(t, y) with a ≤ 0 as in proposition 2 and moreover assume that b(t, y) ≥ 0. Its numerical solution with the EAB k scheme is denoted y n . For the EAB 2 scheme, y n ≥ 0 for n ≥ 2 if:

h ≤ 1/|a| and 0 ≤ e ah y 1 -hϕ 2 b 0 , (33)
For the EAB 3 sceme:

y n ≥ 0 for n ≥ 3 if h ≤ β 3 /|a| and 0 ≤ e ah y 2 -2h(ϕ 2 + ϕ 3 1 , (34)
for a constant β 3 0.331.

Proof of lemma 3. In the present case we have g n (t, y) = b(t, y). For the EAB 2 scheme we have

y n+1 = e ah y n -hϕ 2 b n-1 :=sn + h ϕ 1 b n + ϕ 2 b n ≥0 .
It suffices to prove that s n ≥ 0, which is done by induction. Assume that s n ≥ 0. We have s n+1 = e ah y n+1 -hϕ 2 b n and we get s n+1 = e ah s n + hb n e ah ϕ 1 + e ah ϕ 2 -ϕ 2 = e ah s n + hb n ϕ 1 (e ah +ahϕ 2 ).

By assumption e ah s n ≥ 0, b n ≥ 0 and so hb n ϕ 1 ≥ 0. We then must find the condition so that e ah +ahϕ 2 ≥ 0. A simple function analysis gives ah ∈ [-1, 0] which is the stability condition in lemma 3. The condition [START_REF] Spiteri | Stiffness analysis of cardiac electrophysiological models[END_REF] states that s 1 ≥ 0.

For the EAB 3 scheme we define

s n = e ah y n -2hb n-1 (ϕ 2 + ϕ 3 )
so that we have

y n+1 = s n + hb n ϕ 1 + 3/2ϕ 2 + ϕ 3 + hb n-2 ϕ 2 + ϕ 3
The two terms on the right are non-negative. We will prove by induction that s n ≥ 0. Assume that s n ≥ 0, then:

s n+1 = e ah s n + hb n ϕ 1 e ah +ϕ 2 (3/2 e ah -2) + ϕ 3 (e ah -2) + hb n-1 ϕ 2 /2 + ϕ 3 .
In the right hand side, the first and last terms are non-negative. Thus the condition for s n+1 to be non-negative is that ϕ 1 e ah +ϕ 2 (3/2 e ah -2) + ϕ 3 (e ah -2) ≥ 0.

The study of that inequality gives h ≤ β 3 /|a| and a numerical evaluation of β 3 is β 3 0.331. With the initial condition s 2 ≥ 0 we recorev [START_REF] Tusscher | A model for human ventricular tissue[END_REF].

Proof of proposition 2. Consider the EAB 2 scheme. We define u n = y n -K 1 .

We have

u n+1 =(y n -K 1 ) + h ϕ 1 a(y n -K 1 ) + (b n + aK 1 ) + ϕ 2 (b n + aK 1 ) -(b n-1 + aK 1 )
=u n + h ϕ 1 au n + bn + ϕ 2 bnbn-1 , with bj = b j + aK 1 . Thus u n is the numerical solution with the EAB 2 scheme of the equation u = au + b(t, u) with b(t, u) = b(t, u) + aK 1 . By assumption b(t, u) ≥ 0 and the lemma 3 applies to u n . It gives the lower bound in condition [START_REF] Perego | An efficient generalization of the Rush-Larsen method for solving electro-physiology membrane equations[END_REF]. The upper bound is similarly proven by considering v n = K 2 -y n . The proof is the same for the EAB 3 scheme.

Numerical results

We present in this section numerical experiments that investigate the convergence, accuracy and stability properties of the I-EAB k and EAB k schemes. The membrane equation in cardiac electrophysiology is considered for two ionic models, the Beeler-Reuter (BR) and to the ten Tusscher et al. (TNNP) models. We refer to [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF] and [START_REF] Tusscher | A model for human ventricular tissue[END_REF] for the definition of the models. The stiffness of these two models is due to the presence of different time scales ranging from 1 ms to 1 s, as depicted on figure 5. The stabilizer a n always is a diagonal matrix in this section. The membrane equation has the general form, see [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF][START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF][START_REF] Luo | A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes[END_REF][START_REF] Tusscher | A model for human ventricular tissue[END_REF]:

dw i dt = w ∞,i (v) -w i τ i (v) , dc dt = q(w, c, v), , dv dt = -I ion (w, c, v)+I st (t), ( 35
)
where w = (w 1 , . . . , w p ) ∈ R p is the vector of the gating variables, c ∈ R q is a vector of ionic concentrations or other state variables, and v ∈ R is the cell membrane potential. These equations model the evolution of the transmembrane potential of a single cardiac cell. The four functions w ∞,i (v), τ i (v), q(w, c, v) and I ion (w, c, v) are given reaction terms. They characterize the cell model. The function I st (t) is a source term. It represents a stimulation current.

The formulation (3) is recovered with,

a(t, y) =   -1/τ (v) 0 0 0 0 0 0 0 0   , b(t, y) =   w ∞ (v)/τ (v) q(y) -I ion (y) + I st (t)   , for y = (w, c, v) ∈ R N (N = p+q+1) and where -1/τ (v) = diag (-1/τ i (v)) i=1...p and w ∞ (v)/τ (v) = w ∞,i (v)/τ i (v) i=1...p .

Implementation and computational cost

The computation of y n+1 with the I-EAB k and EAB k schemes requires the data y n-i , a n-i := a(t n-i , y n-i ) and b n-i := b(t n-i , y n-i ) for i = 0 . . . k.

EAB k practical implementation

Firstly, the g ni = b n-i + (a n-i -a n )y n-i are updated at each time step. Then the coefficients γ nj in table 1 are computed. Secondly, the computation of y n+1 by formula (8) also requires the computation of the ϕ j (a n h)γ nj . This is a matrix-vector product in general. In the present case of a diagonal stabilizer, it becomes a scalar-scalar product per row. The ϕ j (a n h) are computed on all diagonal entries of a n h. This computation simply necessitates to compute ϕ 0 (a n h) = e anh (one exponential per non zero diagonal entry) thanks to the recursion rule [START_REF] Colli-Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models[END_REF]. In general, the relation ( 9) can be used to replace the computation of the ϕ j (a n h)γ nj for j = 0 . . . k by the computation of a single product ϕ k (a n h)w k . Denoting by w 1 = a n y n + b n and w j = γ nj + a n hw j-1 :

EAB 2 : y n+1 = y n + h [w 1 + ϕ 2 (a n h)w 2 ] , EAB 3 : y n+1 = y n + h [w 1 + w 2 /2 + ϕ 3 (a n h)w 3 ] , EAB 4 : y n+1 = y n + h [w 1 + w 2 /2 + w 3 /6 + ϕ 4 (a n h)w 4 ] .

I-EAB k practical implementation

In addition, the I-EAB k method (10) requires a quadrature rule of sufficient order to preserve the scheme accuracy and convergence order. We used the Simpson quadrature rule for the cases k = 2, 3 and the three point Gaussian quadrature rule for k=4. We point out that a n is assumed diagonal here so that the matrix exponentials below actually are scalar exponential. The I-EAB k method with Simpson quadrature rule reads, y n+1 = e g1 (y n + b n h/6) + b1 + 4 e δ b1/2 h/6, these operations, especially for large time steps when convergence is harder to reach.

In addition, the I-EAB k and EAB k schemes need several specific operations. In the case of a diagonal function a(t, y) they have been previously described: the EAB k require one scalar exponential computation per non zero row of a(t, y), the I-EAB k with Simpson rule needs twice more and the I-EAB 3 with 3 point Gaussian quadrature rule four times more. Such a cost is not negligible, but is at worst of same order than computing (t, y) → f (t, y) for complex models. For the TNNP model considered here, computing (t, y) → f (t, y) costs 50 scalar exponentials whereas the EAB k implementation adds 7 supplementary scalar exponentials per time step. In terms of cost per time step, the EAB k method is rather optimal. The relationship between accuracy and cost of the EAB k method has been investigated in [START_REF] Douanla-Lontsi | Efficient high order schemes for stiff ODEs in cardiac electrophysiology[END_REF]: more details are available in section 6.4. On every interval (t 3n , t 3n+3 ) the polynomial y of degree at most 3 so that y(t 3n+i ) = y 3n+i , i = 0 . . . 4 is constructed. On (0, T ), y is a piecewise continuous polynomial of degree 3. Its values at the reference time instants nh ref are computed. This provides a projection P (y) of the numerical solution y on the reference grid. Then P (y) can be compared with the reference solution y ref . The numerical error is defined by

Convergence

e(h) = max |v ref -P (v)| max |v ref | , (36)
where the potential v is the last and stiffest component of y in equation 35.

The convergence graphs for the BR model are plotted on figure 6. All the schemes display the expected asymptotic behavior e(h) = O(h k ) as h → 0, as proved in theorem 1.

Stability

The stiffness of the BR and TNNP models along one cellular electrical cycle (as depicted on figure 5 has been evaluated in [START_REF] Spiteri | Stiffness analysis of cardiac electrophysiological models[END_REF]. The largest negative real part of the eigenvalues of the Jacobian matrix during this cycle is of -1170 and -82 for the TNNP and BR models respectively. This means that the TNNP model is 15 times stiffer than the BR model ( 15 1170/82). We want to evaluate the impact of this increase of stiffness in terms of stability for the EAB k and I-EAB k schemes and to provide a comparison with some other classical time stepping methods. To this aim we consider the critical time step ∆t 0 . It is defined as the largest time step such that the numerical simulation runs without overflow nor non linear solver failure for h < ∆t 0 . The numerical evaluation of ∆t 0 is easy for explicit methods. For implicit methods, the choice of the non linear solver certainly impacts ∆t 0 . Without considering more deeply this problem, we just carefully set up the non linear solver, so as to provide the largest ∆t 0 . In practice, we have been using a Jacobian free Krylov Newton method. Results are on table 2. The Adams-Bashforth (AB k ) and the backward differentiation (BDF k ) methods of order k have been considered, together with the RK 4 scheme. The AB k and the RK 4 schemes have bounded stability domain (see [15, p. 243]). Then it is expected for the critical time step to be divided by a factor close to 15 between the BR and TNNP models. Results presented in table 2 show this behavior. The BDF 2 scheme is A-stable whereas the BDF 3 and BDF 4 are A(α)-stable with large angle α (see [15, p. 246]). Hence the critical time step is expected to remain unchanged between the two models. 2. The critical time steps for the I-EAB k schemes remain almost unchanged from the BR to the TNNP model. For the EAB k , they are divided by approximately 2, which behavior is similar as for the BDF k method. As a conclusion, for the present application, the EAB k and I-EAB k methods are as robust to stiffness than the implicit BDF k schemes, though being explicit. As a matter of fact, section 4 shows that the stability domains for the I-EAB k and EAB k schemes depend on the discrepancy between the complete Jacobian matrix and a(t, y). In the present case, a(t, y) only contains a part of the Jacobian matrix diagonal. It is very interesting to notice that robustness to stiffness is actually achieved with this choice. It is finally also interesting to see that the critical time steps of implicit and exponential methods are of the same order.

Accuracy

In terms of accuracy, the schemes can be compared using the relative error e(h) in equation 36. The EAB k and I-EAB k schemes can be compared with the AB k methods only at very small time steps, because of the lack of stability of AB k schemes (see table 2). In table 3 are given the accuracies of these methods for a given time step h = 10 -3 and for the BR model. Ii is observed that the same level of accuracy is obtained with AB k and EAB k at fixed k. These figures illustrate that inside the asymptotic convergence region, the EAB k , I-EAB k and AB k schemes are equivalent in terms of accuracy.

Comparison at large time steps between the EAB k and BDF k for the TNNP model is shown in table 4. These figures show that for large time steps BDF k is more accurate than EAB k . A gain in accuracy of factor 2.5, 5 and 6 is observed for h = 0.05 and for k=2, 3 and 4 respectively. However, compare row 3 for EAB k (h = 0.025) with row 2 for BDF k (h = 0.05). It shows that the numerical solutions with an accuracy close to 0.01 are obtained when dividing the time step by (less than) 2 between BDF k and EAB k . Meanwhile, EAB k with h = 0.025 costs less than BDF k with h = 0.05, as developed in section We conclude that EAB k schemes provide a cheaper way to compute numerical solutions at large time step for a given accuracy. The same conclusion also holds for the BR model, see table 5. A deeper analysis of the relationship between accuracy and computational cost for the EAB k scheme as compared to other methods is available in [START_REF] Douanla-Lontsi | Efficient high order schemes for stiff ODEs in cardiac electrophysiology[END_REF] with the same conclusion. In table 5 are given the accuracies at large time step now considering the BR model. Comparison with table 4 shows that accuracy is preserved by dividing h by 2 when switching from the BR to the TNNP model. As already said, the TNNP model is 15 times stiffer than the BR model. We conclude that the EAB k schemes also exhibit a large robustness to stiffness in terms of accuracy. This robustness is equivalent as for the implicit BDF k schemes. This is remarkable for an explicit scheme, as for the robustness to stiffness in terms of critical time step discussed in the previous subsection.

3. 4 . 2

 42 Case of the I-EAB k scheme We denote ã = ã[tn,Y (tn),h] and b = b[tn,Y (tn),h] . The local error (18) for the I-EAB k scheme satisfies (t n , h) = |δ(t n+1 )| E with δ = y -z and where z is defined by z = ãz + b, z(t n ) = y(t n ), so that with (24) we have s(t n , Y (t n ), h) = z(t n+1 ). The function δ is defined with δ(t n ) = 0 and, δ = ãδ + r, r(τ ) := (a(τ, y(τ )) -ã(τ ))y(τ ) + (b(τ, y(τ )) -b(τ )).

Figure 1 :

 1 Figure 1: Stability function ρ θ (z) for z ∈ R -, for various values of θ and for the three schemes EAB 2 , EAB 3 and EAB 4 .

Figure 2 :θFigure 3 : 20 Figure 4 :

 23204 Figure 2: EAB 2 : isolines ρ θ (z) = 1 for two different ranges. The stability domain D θ is on the left of the isoline.

Figure 5 :

 5 Figure 5: Two consecutive action potentials for the TNNP model: transmembrane potential v (left) and the fast sodium current I Na (right), that is the main component of I ion during the fast upstroke of the action potential.

Figure 6 :

 6 Figure 6: Relative L ∞ error e(h) for the I-EAB k and EAB k schemes, k = 2, 3 and 4, and for the BR model.

Table 3 : 5 . 32 ×

 3532 Accuracy e(h) for the AB k , I-EAB k and EAB k schemes: using the BR model and fixed time step h = 10 -3 10 -6 4.33 × 10 -8 8.69 × 10 -10 I-EAB k 8.55 × 10 -6 4.44 × 10 -8 7.30 × 10 -10 EAB k 7.90 × 10 -6 7.00 × 10 -8 1.16 × 10 -9

3 (

 3 9.01 × 10 -2 5.59 × 10 -2 8.93 × 10 -2 0.025 2.14 × 10 -2 7.34 × 10 -3 8.34 × 10 -3.57 × 10 -2 1.15 × 10 -2 1.44 × 10 -2 0.025 1.10 × 10 -2 2.58 × 10 -3 2.38 × 10 -3 6.1.

3 (

 3 26 × 10 -2 9.17 × 10 -2 0.119 0.05 8.20 × 10 -2 1.09 × 10 -2 8.96 × 10 -74 × 10 -2 4.09 × 10 -2 4.98 × 10 -2 0.1 3.44 × 10 -2 1.04 × 10 -2 1.27 × 10 -2 0.05 9.74 × 10 -3 2.29 × 10 -3 2.02 × 10 -3

Table 1 :

 1 Coefficients γ nj for the EAB k schemes

	k	1 2	3	4
	γ n1 g n g n	g n	g n
	γ n2

  Bashforth integrators for stiff ODEs an absolute discrepancy lower than 25 %, 10 % and 5 % for k =2, 3 and 4 respectively.

	15,
	chapter V.2]. In contrast, EAB

k and I-EAB k methods exhibit much better stability properties. When θ is close enough to 1, they are A(α) stable and stiff stable. Such stability properties are comparable with those of implicit linear multistep methods. Basically speaking, these properties will hold if the stabilization term a(t, y) in (3) approximates the Jacobian of f (t, y) in (1) with  Exponential Adams-

Table 2

 2 shows that the ∆t 0 actually are divided by approximately 2. The critical time steps for the I-EAB k and EAB k models are presented in

table  Table 2 :

 2 Critical time step ∆t 0

		(a) Classical methods	(b) I-EAB k and EAB k exponential methods
		BR	TNNP		BR	TNNP
	AB 2	0.124 × 10 -1 0.850 × 10 -3	I-EAB 2 0.121 0.103
	BDF 2 0.306	0.158	EAB 2	0.424 0.233
	AB 3	0.679 × 10 -2 0.464 × 10 -3	I-EAB 3 0.103 0.123
	BDF 3 0.362	0.181	EAB 3	0.203 0.108
	AB 4	0.372 × 10 -2 0.255 × 10 -3	I-EAB 4 0.133 0.106
	RK 4	0.338 × 10 -1 0.255 × 10 -2	EAB 4	0.122 0.756 × 10 -1
	BDF 4 0.423	0.201

Table 4 :

 4 Accuracy for the TNNP model

Table 5 :

 5 Accuracy for the BR model

Acknowledgments. This study received financial support from the French Government as part of the "Investissement d'avenir" program managed by the National Research Agency (ANR), Grant reference ANR-10-IAHU-04. It also received fundings of the project ANR-13-MONU-0004-04.

where (with the notations of section 1.2 g1 = gn (t n+1 ), δ = g1 -gn (t n + h/2), b1 = bn (t n+1 ) and b1/2 = bn (t n + h/2). These coefficients are given for k = 2 by

and for k = 3 by

The I-EAB k method with the three point Gaussian quadrature rule reads,

5 bl e -g l +8 b0 e -g 0 +5 br e -gr , with bs = bn (t s ), gs = gn (t s ) for s ∈ {l, 0, r} where and bl is the radical conjugate of br (the radical conjugate of x + √ y is x -√ y). Finally, the parameters gs definition is:

and gl is the radical conjugate of gr .

Computational cost

Consider an ODE system (1) whose numerical resolution cost is dominated by the computation of (t, y) → f (t, y). This might be the case in general for "large and complex models". For such problems explicit multistep methods are relevant since they will require one such operation per time step. In contrast, implicit methods, associated to a non linear solver, may necessitate a lot of