Benoît Chappet De Vangel

Jérémy Fix
email: jeremy.fix@centralesupelec.fr

In the quest of efficient hardware implementations of dynamic neural fields: an experimental study on the influence of the kernel shape

Dynamic neural field (DNF) is a popular mesoscopic model for cortical column interactions. It is widely studied analytically and successfully applied to physiological modelling, bioinspired computation and robotics. DNF behavior emerges from distributed and decentralized interactions between computing units which makes it an interesting candidate as a cellular building-block for unconventional computations. That is why we are studying the hardware implementation of DNF on digital substratum (eg. FPGA). As shown in previous papers, this implementation requires several modifications to the equations in order to obtain decent hardware surface utilisation and clock speed.

Here we show that the modification of the lateral weights kernel function is possible as long as certain conditions, enumerated in Amari's seminal work are respected. Thank to metaheuristic optimisation it is possible to find the right parameters for two behavioral scenarii of bio-inspired computation interest. We show that the two most hardware-friendly kernels (difference of linear functions and piece-wise function) are as easy to tune as the traditional Mexican hat kernel. However the difference of exponential kernel is more difficult to tune.

I. INTRODUCTION

Dynamic neural fields (DNF) is a mean field model for cortical neural population behavior. In DNF, both excitatory and inhibitory neural population share the same equation and the propagation delays are neglected for the afferent and lateral contributions. It has interesting applications for physiological modeling [START_REF] Erlhagen | Dynamic field theory of movement preparation[END_REF], autonomous robotics [START_REF] Engels | Dynamic fields endow behavior-based robots with representations[END_REF]- [START_REF] Vitay | A distributed model of spatial visual attention[END_REF] and bio-inspired computing like multi-level clustering [START_REF] Jin | A new clustering approach on the basis of dynamical neural field[END_REF], classification [START_REF] Cerda | Bio-inspired visual sequences classification[END_REF] or tracking [START_REF] Rougier | Emergence of attention within a neural population[END_REF]. This wide spectrum of application makes it an interesting tool for unconventional computing substratum. Unconventional computing is the field which study new computing paradigms to escape from von-Neumann bottleneck. New hardware have to be developed together with new algorithms both often based on a biological inspiration. Unconventional computing include quantum computing, chemical computer for reaction diffusion computations, cellular computing, membrane computing etc. Neuro inspiration can also be included in unconventional computing with neuro-morphic hardware [START_REF] Mead | Neuromorphic electronic systems[END_REF] together with neuron-based computations. For most of these new computing paradigms, the common idea is to have massively parallel computation in a decentralized and unsupervised fashion. The resulting behavior is an emergent property of the hardware substratum. It is attractive because it is generally very robust, scalable and fast on a certain category of problems [START_REF] Sipper | The emergence of cellular computing[END_REF]. We are studying the dynamic neural fields as a candidate for emergent computation with efficient hardware implementation [START_REF] Chappet De | Randomly spiking dynamic neural fields[END_REF] and by studying the spectrum of DNF computation abilities.

After the seminal work of Amari, a lot of analytical results were shown on the dynamical properties of the DNF. Bump stability, oscillations, travelling waves etc [START_REF] Coombes | Waves, bumps, and patterns in neural field theories[END_REF], [START_REF] Bressloff | Spatiotemporal dynamics of continuum neural fields[END_REF]. On the contrary, in [START_REF] Schöner | Dynamics of behavior: Theory and applications for autonomous robot architectures[END_REF], [START_REF] Rougier | Emergence of attention within a neural population[END_REF], [START_REF] Fix | A dynamic neural field approach to the covert and overt deployment of spatial attention[END_REF], the concern was more on the utilisation of the DNF as a building block for bio-inspired cognition behaviors with visual attention, working memory and tracking working together to obtain the right behavior.

Hardware implementation of DNF is not straightforward. The heaviside activation function has to be use to diminish the bandwidth (some time even with a change in the potential equation to reduce the activation time with a spiking behavior like in [START_REF] Vazquez | Visual attention using spiking neural maps[END_REF]). But the main issue persists: the global connectivity between the computing units. As it is not imaginable to connect so many units together we are looking for a von-Neumann neighbourhood communications between our units thus facilitating the hardware implementation. In previous work we proposed an approach to bypass the global connectivity in a cellular way by using a stochastic bit stream approach [START_REF] Chappet De | Stochastic and asynchronous spiking dynamic neural fields[END_REF]. A direct consequence was a modification of the lateral-weights function of the DNF from a difference of Gaussian to a difference of exponential.

More generally, for any hardware implementation we often want to modify the equations and mainly the lateral weights kernel equation. Are the main properties of the DNF still present with unconventional lateral weights kernel shapes?

To answer this problematic we will use meta-optimisation to show that it is possible to find parameters which preserve the DNF behaviors for the most interesting abilities of the DNF: selection and memory.

II. METHODS

A. Dynamic neural fields

Dynamic neural fields were introduced in [16]- [START_REF] Beurle | Properties of a mass of cells capable of regenerating pulses[END_REF] as a mean-firing rate model of coupled populations of inhibitory and excitatory neurons. The evolution of the potential u(x, t) of a neural population reads:

τ ∂u ∂t (x, t) = -u(x, t) + y w(x, y)f (u(y, t)) + I(x, t) + h (1)
where f is an activation function (generally a sigmoid or a step function), I is the afferent input and h is the resting potential. In this paper we define the transfer function f , which links the firing rate to the membrane potential, as a Heaviside function (f (x) = H(x)) as its implementation is much easier in hardware. The lateral weights kernel w(x, y) is responsible for the lateral interactions within the field.

In [START_REF] Amari | Dynamics of pattern formation in lateral-inhibition type neural fields[END_REF], [START_REF] Guo | Existence and stability of standing pulses in neural networks: I. existence[END_REF] the lateral weights kernel general shape is described to ensure convergence and stability: 1) w(x,y) depends only on x-y (homogeneous); 2) w(x,y) is positive around its origin (excitatory part); 3) w(x,y) is negative and tends towards 0 outside this neighbourhood (inhibitory part); 4) w(x) = w(-x) (symmetry).

A traditional choice can be a Mexican hat function (see fig. 1a) but we will study in this paper different other shapes which are convenient for hardware and software implementations. These weight functions keep the above properties and roughly the excitatory and inhibitory influences of the Mexican hat.

(x) = A e exp -x 2 2σ 2 e -A i exp -x 2 2σ 2 i (2)
b) Wizard hat or difference of exponential (DoE): this is a kernel that we already studied experimentally and implemented in hardware in [START_REF] Chappet De | Randomly spiking dynamic neural fields[END_REF]. It has also been studied analytically by several authors because it is easily tractable [START_REF] Coombes | Waves, bumps, and patterns in neural field theories[END_REF], [START_REF] Guo | Existence and stability of standing pulses in neural networks: I. existence[END_REF].

w(x) = A e exp -4|x| σ 2 e -A i exp -4|x| σ 2 i (3)
The terms to the power of x have been defined in order to keep a shape of the kernel similar to the difference of Gaussian when using the same parameters. c) Difference of linear function (DoL): this is a new kernel that we will use in a new architecture. The linearity of the function greatly increases the computation speed in a hardware substratum.

w(x) = A e 1 - |x| 2σ e + -A i 1 - |x| 2σ i + (4)
with A e > A i , σ i > σ e and [.] + denotes the positive part of the inner term.

d) Rectangular Mexican Hat or step function (

Step): already used in [START_REF] Giese | Technical applications of neural fields[END_REF] to generate a working memory behavior, the step function is a piecewise constant kernel whose equation reads:

w(x) = A e H(|x| -σ e ) -A i H(|x| -σ i ) (5)
with A e > A i , σ i > σ e and H denotes the Heaviside function. This is the cheapest implementation possible as its computation is based on simple comparisons.

B. Optimization problem

Dynamic neural fields have been shown to exhibit various properties such as the ability to create a localized decision (a so-called bump) where the input is locally maximal as well as to retain an information and therefore to act as a memory of some gated inputs. In order to determine if a dynamic neural field with the previously introduced kernel functions are still able to exhibit these properties, we follow the optimisation procedure of [START_REF] Fix | Template based black-box optimization of dynamic neural fields[END_REF]. It consists in introducing some scenarios with a set of inputs and particular cost functions to be minimized. The cost function or fitness captures the expected behavior when the neural field is fed with the proposed input. The optimization of the scenarios is performed with the particle swarm optimizer (PSO) 1 , a metaheuristic which only requires to evaluate the cost function to be minimized and not any of its derivatives which is handful given the non linearities in the cost function and neural field equations. In this paper we consider two scenarios: a competition scenario and a working memory scenario. In the competition scenario, the neural field has five excited regions among which one region is slightly more excited than the others, as depicted on fig. 2. The input consists in four Gaussian stimuli of variance σ = 4 and amplitude A ∈ {0.4, 0.6, 0.8} and one Gaussian stimulus with the same variance but amplitude A + 0.2.

The cost function to be minimized is designed to ensure that a bump appears in the neural field at the location where the input is maximal while vanishing when the input is switched off and, formally, reads :

J = x,|x-x0|≤σ (1 -f (u(x, T 2 )) + x,|x-x0|>σ f (u(x, T 2
))

+ x f (u(x, T ))
1 but any other black box optimizer such as CMA-ES would have been well suited where T = 40 denotes the maximal time of the scenario and x 0 the mean of the input stimulus with the highest amplitude. The cost function simply captures the requirement to get a bump centered on x 0 at half of the scenario, i.e.

f (u(x, T 2 )) = 1 for |x -x 0 | ≤ σ and f (u(x, T 2
)) when |xx 0 | > σ, and no firing rates at the end of the scenario and hence the term

x f (u(x, T )) to be minimized. In order to get a neural field that is more robust than one we would obtain by optimizing on a single input condition, we simultaneously optimize three conditions with various input amplitudes. Indeed, if the neural field is optimized on a scenario with strong inputs it might not be able to perform a competition when the inputs are weak. Therefore, three conditions are optimized simultaneously with the amplitude of the weak stimuli set to A ∈ {0.4, 0.6, 0.8} and the amplitude of the strong stimulus set to A + 0.2. Intuitively, such a behavior is obtained with lateral weights that are essentially globally inhibitory and locally weak but not necessarily excitatory as a positive firing rate is sufficient to get a bump with a Heaviside transfer function. It might be noted that if a single scenario was considered, it would have been straightforward to find optimal parameters. Indeed, setting the baseline to h = -(A + 0.1) and discarding the lateral weights should produce parameters pretty close to an optimal solution. Optimizing simultaneously three scenarios discards such a trivial solution.

For the working memory scenario, two regions are excited in the input with an initial weak amplitude. At t = 30, the leftmost stimulus get transiently excited to trigger the entrance in memory. At t = 50, the second stimulus is triggered. At t = 75, the right most stimulus starts moving to the right by δx = 20, ending at t = 155. The neural field is expected to track the motion of the input and to update its content respectively. The scenario ends at t = 195 after a period during which the inputs are switched off and the stimuli are expected to disappear from the memory. The scenario is challenging for a neural field because the memory requires stronger lateral excitation than the competition and this strong lateral excitation might actually induce such a strong local positive feedback that the field might get blind to motions of input stimuli. The inhibitory component of the lateral weights must be local and rather strong to keep the positive firing rates local to the memorized stimulus and to avoid interference with the other locations in the field. In order to capture this behavior, an appropriate cost function is designed. Its definition is similar to the one for the competition for capturing the presence or absence of a bump. However, it is here built from the contribution of the neural field activities at 5 stages :

• before the first stimulus gets focused, we minimize

x f (u(x))
• after the first stimulus was focused, we minimize

x,|x-x1|≤σ (1 -f (u(x)) + x,|x-x1|>σ f (u(x))
• after the second stimulus was focused, we minimize

x∈|x-x1|≤σ∪|x-x2|≤σ (1 -f (u(x))) + x / ∈|x-x1|≤σ∪|x-x2|≤σ f (u(x)) • after the second stimulus moved to x 2 = x 2 + dx,
while the first stimulus is at

x 1 = x 1 , we minimize x∈|x-x 1 |≤σ∪|x-x 2 |≤σ (1 -f (u(x))) + x / ∈|x-x 1 |≤σ∪|x-x 2 |≤σ f (u(x))
• at the end of the scenario, we minimize

x f (u(x))
For all the simulations, a time discrete neural field equation is used :

u(x, t + 1) = u(x, t) + ∆u(x, t) τ ∆t ∆u(x, t) = -u(x, t) + y w(x, y)f (u(y, t)) + I(x, t) + h
The cost functions introduced previously are optimized with respect to 6 parameters: the time constant ∆t τ , baseline h as well as a specific parametrization of the weight functions to ensure the locally-excitatory/globally-inhibitory shape of the kernel. The kernels have been introduced with A e , σ e for the excitatory component and A i , σ i for the inhibitory component. Actually, to ensure the appropriate shape of the kernels, we set

A i = k a A e , k a ∈ [0, 1] and σ e = k σ σ i , k σ ∈]0, 1].
Therefore, the cost functions are optimized with respect to ∆t τ , h, A e , k σ , k a , σ i with the bounds given in the table below.

Parameter Bounds

∆t τ [0, 0.3] h [-1, 1] A e [0, 5] k σ [1e -3 , 1] k a [0, 1] σ i [1, 100]
C. Optimization algorithm : Particle Swarm Optimization

As in [START_REF] Fix | Template based black-box optimization of dynamic neural fields[END_REF], the cost functions previously introduced are minimized using the Particle Swarm Optimization (PSO) algorithm. The PSO algorithm is an optimization heuristic inspired by the flocking behaviors of birds. Formally, it consists in considering a collection of N so called particles which have a position p i (t) in the parameter space which evolves according to the following equations :

p i (t + 1) = p i (t) + v i (t + 1) v i (t + 1) = wv i (t) + c 1 r 1 (p l i (t) -p i (t)) + c 2 r 2 (p g i (t) -p i (t))
where w, c 1 , c 2 are parameters of the algorithm, r 1 , r 2 are vectors in [0, 1] n randomly generated for every particle at every iteration, p l i (t) is the best position the particle i ever had and p g i (t) is the best position ever occupied by a particle in the neighborhood of the i-ith particle. The particles are not necessarily interacting directly with all the other particles and a neighborhood is introduced. Several neighborhoods have been considered in the literature, for example all-to-all, ring, von Neumann (grid with cyclic boundary conditions), random, etc. Here we stick to the standard PSO algorithm 2006 as defined in [START_REF] Clerc | Standard particle swarm optimization[END_REF] and therefore use a random informants topology where each particle informs K = 3 other particles. In addition, at every iteration, if no improvement in the best fitness is observed, the topology is regenerated. The parameters are set to w = 1 2 log 2 and c 1 = c 2 = 1 2 + log 2 as suggested in [START_REF] Clerc | Standard particle swarm optimization[END_REF].

III. RESULTS

A. Competition

The competition scenario was optimized with n = 20 particles during 100 epochs (i.e. 2000 evaluations of the cost functions) and repeated 1000 times. In order to get a robust analysis, we keep only the trials which ended with a cost function at 0. This is clearly a conservative choice but actually, most of the runs successfully reached a null cost. Out of the 1000 trials, the success rate was respectively 98.9 % for the DOG, 84.1 % for the DOE, 95.2 % for the DOL and 95.2% for the step weight function. On the figure 3 we plot the mean fitness and their standard deviation as a function of the epoch of PSO for all the successful trials for each weight kernel.

To better compare the evolution of the cost functions for each kernel, the mean cost function of the successful trials for all the kernels is shown on the figure 4. The optimization is slower for the step function and harder for the wizard hat (difference of exponential). More interestingly, there are common patterns that can be observed from the distribution of the optimal parameters. We shall comment first the distribution of the parameters that are not shown in the paper but the optimal parameters are available online [START_REF] Fix | The source codes and results for the simulations of the paper[END_REF]. There is not much to say about the time constant ∆t τ which is confined in [0.15, 0.3]. The optimal time constant is clearly influenced by the dynamics that the scenario imposes as the field is requested to build up its decision within the T /2 first time steps. For all the kernels, the baseline h is confined within [-0.5, 0], monomodal with a mean around h 0 ≈ -0.15. For the weights, there are also common patterns. The amplitude of the inhibitory component is always close to the amplitude of the excitatory component with k a ≈ 1. For the DOG, DOE and DOL, the excitatory amplitude tends to be uniformly distributed within the bounds while the variance of the excitatory component tends to be half the variance of the inhibitory component and the variance of the inhibitory component tends to be uniformly distributed within the bounds but lower bounded by 20 which is clearly linked to the inter-stimulus distance. The distribution of the parameters is different for the step function for which the excitatory amplitude A e is much more concentrated to small values, the variance of the excitatory component is strongly concentrated around one tenth of the inhibitory variance and the inhibitory variance is still uniformly distributed within the bounds but lower bounded by 40. As a summary, the mean weights and their standard deviation, still for the successful trials, are depicted on the figure 5. The plotted mean to which the standard deviation is added can be sometimes misleading. For example, for the step function, none of the weights have a positive component for large distances but the distribution of the weight values are actually strongly asymmetric, densely packed around small values but with strong negative components. The plots of all the considered weights are given at [START_REF] Fix | The source codes and results for the simulations of the paper[END_REF].

B. Working memory

The working memory scenario was optimized with n = 200 particles during 1000 epochs (i.e. 200000 evaluations of the cost function) and repeated 1000 times 2 . Contrary to the selection scenario, the choice of the threshold cannot be made too conservative otherwise too few runs would be considered. The 1000 runs of the 4 kernel setups (i.e. a total of 4000 neural fields) were ordered by increasing fitness and the behavior of all the runs were visually inspected in order to set a threshold on the fitness of the trials to consider. It turns out that all the runs with a fitness lower than 8 successfully performed the scenario while all the scenarios with a fitness larger failed with a systematic behavior. The neural fields of the failing runs were all able to memorize the input stimuli after the transient increase of amplitude but failed to track the change of position and lost the moving stimulus. Therefore, a fitness threshold of 8 was considered. Out of the 1000 trials, the success rate were respectively 100 % for the DOG, 46 % for the DOE, 100 % for the DOL and 100% for the step weight function. It is clear that the neural fields are much harder to optimize with the DOE. On the figure 6 we plot the mean fitness and their standard deviation as a function of the epoch of PSO for all the successful trials for each weight kernel. 2 one trial took approximately one hour on an Intel(R) Xeon(R) 3075 at 2.66GHz. As a comparison, a trial with the competition scenario took one minute. To better compare the evolution of the cost functions for each kernel, the mean cost function of the successful trials for all the kernels is shown on the figure 7. The speed of convergence of the optimization is similar whatever the kernel function. However, as mentioned previously, in the case of the DOE (wizard hard) only half of the runs successfully converged to the desired behavior. The optimization problem is harder for the difference of exponential probably because the kernel is very pointy. In order to correctly track a moving stimulus, the kernels has to spread excitation around excited regions in order to stimulate the propagation of excitation within the field, an unfavorable requirement for the difference of exponential.

As in the competition scenario, there are common patterns in the optimal parameters. From the histograms of the optimal parameters (not shown but available at [START_REF] Fix | The source codes and results for the simulations of the paper[END_REF]) the time constant ∆t τ is usually distributed within [0.15, 0.3] but tends to be more packed toward large values (compared to the competition scenario) especially for the step kernel. The baseline h is still monomodal but much densely centered around -0.5 except for the step kernel for which the mean value is -0.33. For all the kernels, there is a clear tendency to get small amplitudes for the excitatory components. The excitatory component tends to be local and usually 40% to 90% of the variance of the inhibitory component. The inhibitory component amplitude tends to be 70% of the excitatory component for the DOG and DOL and more around 90% for the step. For the DOE, the inhibitory amplitude has much more freedom and tends to be uniformly distributed between 10% and 90% of the excitatory amplitude. The variance of the inhibitory component is rather local with mean values around 7 to 14 for the DOG, DOL and DOE and more around 20 for the step. The variance of the inhibitory component has clearly much less freedom for the working memory than for the competition scenario. Indeed, the role of the inhibitory component is to keep the memorized stimulus local while not interacting with the neighboring stimulus that is memorized which imply that the inhibitory component must exist and be local. As a summary, the mean weights and their standard deviation, still for the successful trials, are depicted on the figure 8.

IV. DISCUSSION

Following optimization guidelines introduced in [START_REF] Fix | Template based black-box optimization of dynamic neural fields[END_REF] we were able to find parameters for 4 different kernel on two relevant scenarios. These results are very encouraging for efficient hardware implementation as they indicate that designers have a large choice of lateral weights function shape and precision and can therefore adapt the kernel to the targeted hardware substratum.

For instance, in our previous work we compared centralized hardware dynamic neural field implementation with a decentralized and a cellular approach [START_REF] Chappet De | Spiking dynamic neural fields architectures on fpga[END_REF]. The centralized implantation is based on address event representation (AER) communication between units [START_REF] Mahowald | Vlsi analogs of neuronal visual processing: A synthesis of form and function[END_REF]. We adapted it so as every time a neuron receives an activation information, the distance of the event (and not its address) is provided to the target unit. The modification of the neuron's potential then depends on a local look-up-table which associates a weight to the distance. Present results show that we can dramatically reduce the precision of the weights and even reduce them to a step function with two values using a simple test on the distance thus reducing the memory needs of the implementation.

We also investigated a decentralized, cellular architecture. We are interested by cellular computing as it is a massively parallel architecture with emergent behavior and a good candidate to escape the von Neumann bottleneck as memory and computations are distributed on the same substratum [START_REF] Sipper | The emergence of cellular computing[END_REF]. In this approach the challenge is to design all-to-all lateral interactions with local update rule. It is thus necessary to propagate the activation information sufficiently fast so as the propagation delay can be neglected. In [START_REF] Chappet De | Randomly spiking dynamic neural fields[END_REF] we introduced the randomly spiking dynamic neural fields (RSDNF) were random propagation of the information results in a lateral weights function in a shape of a difference of exponential (DoE) w(d) = k e p d ek i p d i were p e and p i are the probability of propagation for the excitation and inhibitions information respectively. We shown the behavior of this lateral weights function was correct for robust tracking. However the implantation performance in term of speed and area is limited by the generation of a tremendous amount of pseudo random numbers using a cellular automaton which takes half of the area on a FPGA chip. Moreover we show in the present paper that the optimization of the DoE is more difficult.

That is why we are interested by the DoL or step lateral weights kernel which are easier to optimize than the DoE and would not require random numbers generation. The gains in area and speed are not yet evaluated and will be shown in a future work. The first approach is to use a counter to count time since the activation of the units. Assuming that the data is propagated without delay, the cells know the distance of activated unit as it is the value of the counter. With this approach it is easy to derive the lateral weight using a look-uptable to associate one weight with one distance. Using a step function would again replace the look-up-table by a simple test and diminish the memory needs. This approach is cellular but not decentralized as a global synchronization signal has to trigger the start of the propagation. We are currently looking for a way to decentralize this computation as in [START_REF] Chappet De | Stochastic and asynchronous spiking dynamic neural fields[END_REF]. It might however be noticed that using a step function for the lateral weights significantly improves the performance for software implementation. Indeed, when computing the lateral contributions for all the cells with a step function, starting from the lateral contributions for the first cell (involving 3 products and at most N sums), the lateral contributions for the other cells can be deduced, step by step, by only considering the change in contribution of the four units on the splits between the positive/negative/null parts of the weights, which makes the complexity of evaluating such a neural field linear in the number of units rather than N log(N ) with a general convolutive kernel. Future work will also explore analytical tools to derive the lateral weights parameters in order to obtain an expected behavior. This problem known as the inverse problem is well defined and was studied in the context of DNF in [START_REF] Potthast | Inverse problems in neural field theory[END_REF].

To our knowledge similar study does not exist as we are among the first one to tackle the problematic of neural field implementation in digital hardware. However, a lot of work was done to implement spiking recurrent neural networks on digital or analogous substratum using locally cooperative globally competitive connectivity. For example the work of Neftci et al. [START_REF] Neftci | Synthesizing cognition in neuromorphic electronic systems[END_REF], [START_REF] Indiveri | A vlsi neuromorphic device for implementing spike-based neural networks[END_REF] which uses soft winner-take-all (sWTA) network as a building block for a VLSI neuromorphic chip. The sWTA dynamic is the result of global inhibitory connections and local cooperation with excitatory connections to close neighbours. This is similar to our step lateral weights function with infinite w i for global inhibition.

Fig. 1 .

 1 Fig. 1. The different lateral weights function that we are studying in this paper. The difference of Gaussian is our control and the three other are modified to improve hardware implementation. a) Mexican hat or difference of Gaussian (DoG): this is our control for the behavior of the DNF.

  w

Fig. 2 .

 2 Fig.2. In the competition scenario, three conditions (a),(b),(c) are simultaneously optimized with various input amplitudes in order to optimize the ability of the neural field to select the most salient input. In the working memory scenario (d), the neural field is optimized to keep track, in memory, input stimuli that get transiently over-excited.

Fig. 3 .

 3 Fig. 3. Mean and standard deviation of the cost function for every weight kernel. The optimization reached the global minimum of the cost function for 98.9 % of the trials for the DOG (a), 84.1 % for the DOE (b), 95.2 % for the DOL (c) and 95.2% for the step (d) weight function.

Fig. 4 .

 4 Fig. 4. Comparison of the evolution of the mean cost function for every the successful trials of every kernels for the competition scenario.

Fig. 5 .

 5 Fig. 5. Mean and standard deviation of the weights for each of the kernels for the successful trials of the competition scenario. More details about the shapes of the kernels are given in the text.

Fig. 6 .

 6 Fig. 6. Mean and standard deviation over 1000 trials of the cost function for every weight kernel for the working memory scenario. The optimization reached the global minimum of the cost function for 100 % of the trials for the DOG (a), 46 % for the DOE (b), 100 % for the DOL (c) and 100% for the step (d) weight function. The y-scale is identical for the DOG, DOE and DOL but different for the Step function.

Fig. 7 .

 7 Fig. 7. Comparison of the evolution of the mean cost function for every successful trials of every kernels for the working memory scenario.

Fig. 8 .

 8 Fig.8. Mean and standard deviation of the weights for every kernel for the successful trials. More details about the shapes of the kernels are given in the text. The y-scale is the same for all but the step kernel.


David Gasperini

Hans-Peter P Beise

Udo Schröeder

Xavier Antoine

U Schroeder

C Geuzaine

An analysis

Keywords: high-frequency scattering, acoustic single-layer integral operator, steepest descent method, highly oscillatory integrals

Introduction

We consider the problem of the efficient evaluation of the elementary acoustic single-layer integral operator I r,∆ (k) = ∆ e ik( r-r +θ•r ) rr dr , [START_REF] Antoine | An introduction to operator preconditioning for the fast iterative integral equation solution of time-harmonic scattering problems[END_REF] where r ∈ R 3 is an observation point, ∆ is a triangle in R 3 and θ is a given directional unitary vector in R 3 . In addition, we set i := √ -1, and x := √ x • x is the euclidian norm defined through the dot product x • y = 3 j=1 x j y j , for x and y in R 3 . The problem of the computation of (1) arises for example in the framework of the Physical Optics approximation and integral equations [START_REF] Antoine | An introduction to operator preconditioning for the fast iterative integral equation solution of time-harmonic scattering problems[END_REF][START_REF] Chew | Fast and Efficient Algorithms in Computational Electromagnetics[END_REF][START_REF] Nédélec | Acoustic and Electromagnetic Equations[END_REF][START_REF] Wu | The modern high frequency methods for solving electromagnetic scattering problems[END_REF], where the real-valued positive constant k is then the wave number. It is well-known that evaluating [START_REF] Antoine | An introduction to operator preconditioning for the fast iterative integral equation solution of time-harmonic scattering problems[END_REF] for large values of k requires a prohibitive computational cost [START_REF] Chandler-Wilde | Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering[END_REF][START_REF] Wu | The modern high frequency methods for solving electromagnetic scattering problems[END_REF], and reducing the complexity is considered as being a challenging and important numerical problem, most particularly when trying to weaken the k-dependence of the computational cost for example by getting a computation of I r,∆ (k) independent of k. This is the reason why the development of new approximation methods to compute highly oscillatory integrals [START_REF] Deaño | Computing Highly Oscillatory Integrals[END_REF] has received a lot of attention during the last two decades. Among others and without being exhaustive, important achievements in this direction are related to the Levin-type methods [START_REF] Deaño | Computing Highly Oscillatory Integrals[END_REF][START_REF] Levin | Fast integration of rapidly oscillatory functions[END_REF][START_REF] Xiang | On the Filon and Levin methods for highly oscillatory integral[END_REF], the Filon-type discretization technique [START_REF] Asheim | A combined Filon/asymptotic quadrature method for highly oscillatory problems[END_REF][START_REF] Deaño | Computing Highly Oscillatory Integrals[END_REF][START_REF] Gao | Error analysis of the extended Filon-type method for highly oscillatory integrals[END_REF][START_REF] Gao | A generalization of Filon-Clenshaw-Curtis quadrature for highly oscillatory integrals[END_REF][START_REF] Iserles | On quadrature methods for highly oscillatory integrals and their implementation[END_REF][START_REF] Iserles | Efficient quadrature of highly oscillatory integrals using derivatives[END_REF][START_REF] Iserles | Quadrature methods for multivariate highly oscillatory integrals using derivatives[END_REF][START_REF] Iserles | From high oscillation to rapid approximation III: multivariate expansions[END_REF][START_REF] Xiang | On the Filon and Levin methods for highly oscillatory integral[END_REF][START_REF] Zhao | An adaptive Filon-type method for oscillatory integrals without stationary points[END_REF] and the method of steepest descent/stationary phase [START_REF] Asheim | Applying the numerical method of steepest descent on multivariate oscillatory integrals in scattering theory[END_REF][START_REF] Asheim | A remedy for the failure of the numerical steepest descent method on a class of oscillatory integrals[END_REF][START_REF] Asheim | Asymptotic analysis of numerical steepest descent with path approximations[END_REF][START_REF] Huybrechs | On the evaluation of highly oscillatory integrals by analytic continuation[END_REF][START_REF] Wu | The modern high frequency methods for solving electromagnetic scattering problems[END_REF][START_REF] Wu | Computing highly oscillatory physical optics integral on the polygonal domain by an efficient numerical steepest descent path method[END_REF][START_REF] Wu | The numerical steepest descent path method for calculating physical optics integrals on smooth conducting quadratic surfaces[END_REF][START_REF] Wu | Frequency-independent approach to calculate physical optics radiations with the quadratic concave phase variations[END_REF][START_REF] Zhang | Efficient evaluation of the physical-optics integrals for conducting surfaces using the uniform stationary phase method[END_REF]. Other contributions include e.g. the following references [START_REF] Asheim | A Gaussian quadrature rule for oscillatory integrals on a bounded interval[END_REF][START_REF] Bruno | An o (1) integration scheme for three-dimensional surface scattering problems[END_REF][START_REF] Bruno | Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case[END_REF][START_REF] Huybrechs | Highly oscillatory quadrature[END_REF][START_REF] Huybrechs | Superinterpolation in highly oscillatory quadrature[END_REF][START_REF] Huybrechs | The construction of cubature rules for multivariate highly oscillatory integrals[END_REF][START_REF] Iserles | On the computation of highly oscillatory multivariate integrals with stationary points[END_REF][START_REF] Iserles | Quadrature methods for multivariate highly oscillatory integrals using derivatives[END_REF][START_REF] Iserles | From high oscillation to rapid approximation III: multivariate expansions[END_REF]. For completeness, we refer to [START_REF] Deaño | Computing Highly Oscillatory Integrals[END_REF][START_REF] Huybrechs | Highly oscillatory quadrature[END_REF] for a general overview and more references. We focus in this paper on the method of steepest descent [START_REF] Deaño | Computing Highly Oscillatory Integrals[END_REF] for the computation of [START_REF] Antoine | An introduction to operator preconditioning for the fast iterative integral equation solution of time-harmonic scattering problems[END_REF].

Before attacking the full three-dimensional problem, it is standard to start in dimension one. The basic idea of the steepest descent method is then to use the Cauchy Integral Theorem (CIT) to replace the univariate oscillatory integral

I(k, a, b) := b a f (x)e ikg(x) dx (2)
by integrals along paths in the complex plane where e ikg(x) is not oscillatory. Throughout this work, f and g are assumed to be real-valued on [a, b] and analytic in a sufficiently large complex domain that contains [a, b]. Moreover, f is supposed to be non-oscillatory on [a, b]. Now, suppose that an integration is carried out along a path where g has constant real part and positive increasing imaginary part. Then, the key observation is that e ikg(z) is non-oscillatory and exponentially decaying for complex values of g(z), since e ikg(z) = e -kIm(g(z)) e ikRe(g(z)) . Therefore, the steepest descent method is based on building extensions of [a, b] to a complex integration contour which fulfills the previous conditions. Finally, this leads to the following definition of the steepest descent path.

Definition 1. Let us consider the previous notations and let I = [0, T [ be an interval such that T ∈]0, ∞]. A function h : I → C which satisfies g(h(t)) = g(h(0))

+ it = g(x) + it, ( 3
)
for all t ∈ I, is called [START_REF] Huybrechs | On the evaluation of highly oscillatory integrals by analytic continuation[END_REF] a steepest descent path of g on I with trace h * := h(I).

In the following, we denote by h x a steepest descent path starting at x ∈ [a, b] (i.e. h x (0) = x). In the steepest descent method, the goal is to build steepest descent paths originating at the edges of [a, b], and parametrized by a positive real number t, that can be connected by a path ϕ a,b of negligible contribution to the integral. More precisely, we need to find steepest descent paths h a , h b , both being defined on [0, ∞[, and such that for every t ≥ 0, there exists a path ϕ a,b defined on [0, 1] and two real numbers t a , t b ∈ [t, ∞[, with h a (t a ) = ϕ a,b (0) and h b (t b ) = ϕ a,b [START_REF] Antoine | An introduction to operator preconditioning for the fast iterative integral equation solution of time-harmonic scattering problems[END_REF], such that the CIT applies on the induced contour and that the following path integral (z) dz → 0 as t tends to ∞. Then, the integral (2) can be rewritten as

ϕ a,b f (z)e ikg
I(k, a, b) = ha f (z)e ikg(z) dz - h b f (z)e ikg(z) dz = ∞ 0 h a (t)f (h a (t))e ikg(ha(t)) dt - ∞ 0 h b (t)f (h b (t))e ikg(h b (t)) dt. (4)
In [START_REF] Huybrechs | Highly oscillatory quadrature[END_REF], Huybrechs and Olver present a way to write explicitly the steepest descent paths in the two-dimensional case. However, it is not known whether there exist some extensions of such paths for higher dimensional problems, leading to further investigations to handle multivariate integrals [START_REF] Asheim | Applying the numerical method of steepest descent on multivariate oscillatory integrals in scattering theory[END_REF][START_REF] Deaño | Computing Highly Oscillatory Integrals[END_REF][START_REF] Huybrechs | The construction of cubature rules for multivariate highly oscillatory integrals[END_REF][START_REF] Iserles | On the computation of highly oscillatory multivariate integrals with stationary points[END_REF][START_REF] Iserles | Quadrature methods for multivariate highly oscillatory integrals using derivatives[END_REF][START_REF] Iserles | From high oscillation to rapid approximation III: multivariate expansions[END_REF]. Nevertheless, it is shown in [START_REF] Huybrechs | The construction of cubature rules for multivariate highly oscillatory integrals[END_REF] that the method of steepest descent can be applied to multivariate cubatures by means of nested integrals. In order to apply this method to the integrals of the form (1), the restriction to the corresponding univariate form has to be investigated. Here, it turns out that in general a straightforward application of (4) is not feasible since the validity of the CIT does not hold in all settings, or two paths cannot always be linked by a path of negligible contribution. Such issues are also pointed out in [START_REF] Asheim | A remedy for the failure of the numerical steepest descent method on a class of oscillatory integrals[END_REF]. In the latter work, the author shows that a polar change of variables can circumvent these issues in certain cases. Let us finally remark that recent contributions [START_REF] Wu | An efficient method for computing highly oscillatory physical optics integral[END_REF][START_REF] Wu | The modern high frequency methods for solving electromagnetic scattering problems[END_REF][START_REF] Wu | Computing highly oscillatory physical optics integral on the polygonal domain by an efficient numerical steepest descent path method[END_REF][START_REF] Wu | The numerical steepest descent path method for calculating physical optics integrals on smooth conducting quadratic surfaces[END_REF][START_REF] Wu | Frequency-independent approach to calculate physical optics radiations with the quadratic concave phase variations[END_REF][START_REF] Zhang | Efficient evaluation of the physical-optics integrals for conducting surfaces using the uniform stationary phase method[END_REF] developed with success the application of the numerical steepest descent method to the Physical Optics in high-frequency electromagnetic scattering, for both linear and quadratic phases.

In this work, we first consider the univariate version of the integral in [START_REF] Antoine | An introduction to operator preconditioning for the fast iterative integral equation solution of time-harmonic scattering problems[END_REF]. We identify two types of potential critical points on the integration interval that have to be circumvented by splitting the integration: firstly, stationary points of the phase function, and, secondly, points that give steepest descent path that run into one of the two possible singularities of the integrand in the complex plane. We give some explicit formulas for the determination of these splitting points and show that at most one of them can be present in a given bounded interval. We also provide explicit formulas for the steepest descent path. These paths turn out to differ depending on the location of their starting point with respect to a splitting point and depending on the incident direction of the wave with respect to the integration surface. In all settings, we can prove that paths that start on the same side relative to a splitting point allow the application of the steepest descent method, i.e. they can be connected by a path of negligible contribution to the integral and it is ensured that the Cauchy integral theorem applies on the resulting contour. In the investigation of the steepest descent paths in the context of the univariate case, we also explore the method to determine them by rewriting (3) as an ordinary differential equation. This method is further applicable when no explicit formula for (3) is available and therefore can be considered as an alternative to the Taylor's series expansion approach given in [START_REF] Asheim | Asymptotic analysis of numerical steepest descent with path approximations[END_REF]. For the integration over a two-dimensional surface we follow [START_REF] Huybrechs | The construction of cubature rules for multivariate highly oscillatory integrals[END_REF] and write (1) as nested univariate, oscillatory integrals. But in contrast to [START_REF] Huybrechs | The construction of cubature rules for multivariate highly oscillatory integrals[END_REF], where the general case is treated, we only consider a fixed type of integral and it turns out that in our case, the oscillatory components in the two perpendicular directions of integration can be decomposed as a product of two oscillatory integrals with a remaining coupled component that only varies slowly in both variables. The latter can hence can be treated by usual numerical methods. Numerical experiments at the end of the paper demonstrate the efficiency and accuracy of the method for univariate and two-dimensional oscillating integrals given by [START_REF] Antoine | An introduction to operator preconditioning for the fast iterative integral equation solution of time-harmonic scattering problems[END_REF].

The paper is organized as follows. In Section 2, we provide the setting and introduce some notations. Section 3 is devoted to the explicit calculation of the splitting points. The construction and properties of admissible steepest descent paths is analyzed in Section 4. Based on these results, we develop the integration technique in Section 5 for the univariate case first, and then for triangular surfaces. Numerical illustrations are provided in Section 6. Section 7 concludes the paper.

Settings and notations

Let us consider (1), for a triangle ∆ in R 3 and a given unitary vector θ in R 3 . Let us introduce a linear parametrization to map the unit triangle ∆ 1 defined by the corners (0, 0), (1, 0), (0, 1) to ∆ given by the affine mapping:

x → Ax + b, with x = (x, y) T , A := (A ) 1≤ ≤2 ∈ R 3×2 , A ∈ R 3 and b ∈ R 3 .
In this case, (1) writes as

I ∆ (k) = |A 1 × A 2 | 1 0 1-y 0 e ik( Ax+b-r +θ•(Ax+b)) Ax + b -r dxdy =: |A 1 × A 2 | I 1 ∆ 1 (k), (5)
where the index r in I r,∆ (k) is avoided to simplify the notations. Let us define P (x, y) := ||Ax + b -r|| 2 as the second-order polynomial with respect to the two variables x and y, which is also denoted by P (x) when y is temporarily frozen. We assume that the observation point r does not lie in ∆ which implies that P only takes positive real values and hence writes as

P (x) = c 0 (x -c)(x -c) = c 0 x 2 + c 0 |c| 2 -2c 0 xRe(c), (6)
with c 0 ∈ [0, ∞[ and c ∈ C \ R, where c and c are the roots of P . In this section, we only consider the one-dimensional inner integral of I 1 ∆ 1 (k) for some arbitrary but fixed y ∈ [0, 1], the functions and the coefficients being then y-dependent. The linear part in g is also written as a linear function x → θ • (Ax + b) =: qx + s, with real-valued q and s that vary as a function of y. Hence, g is given by g(x) = P (x) + qx + s [START_REF] Bruno | An o (1) integration scheme for three-dimensional surface scattering problems[END_REF] and the associated integral is of the following form

I(k, y, 0, 1 -y) = 1-y 0 e ikg(x) P (x) dx.
To keep a certain generality, we rather study integrals over real intervals [a, b] defined by

I(k, y, a, b) := b a e ikg(x) P (x) dx. ( 8
)
According to Definition 1, we directly see that a steepest descent path h x :

I = [0, T [→ C, with T ∈]0, ∞], starting at x ∈ [a, b] is defined by the equation g(h x (t)) = P (h x (t)) + qh x (t) + s = P (x) + qx + s + it, (9)
with h x (0) = x. In this early state of the work, √ . should denote an analytic continuation of a fixed definition of the square-root at P (x) along the path P (h x (t)) starting at P (x) = P (h x (0)). More precisely, let us define

Π := {z ∈ C | Re(z) = Re(c) and Im(z) ∈]Im(c), Im(c)[} and Π c := Re(c) + i (R \ Π).
One can directly check that the inverse image of ] -∞, 0] under P is Π c . Hence, if we consider paths that do not intersect Π c , one can use the main branch of the square-root in [START_REF] Chandler-Wilde | Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering[END_REF]. We will show how to ensure that the steepest descent paths given in the sequel satisfy the aforementioned condition. As a consequence, the definition of the complex square-root can be fixed to the main branch.

Explicit calculation of the splitting points

To compute (8) by the steepest descent method, we have to ensure that two steepest descent paths can be connected by a path of negligible contribution with regard to the integral (see Section 1). To this end, we need to determine the stationary point of g which is a first splitting point. In addition, we guarantee that the steepest descent paths are well-defined and that the closed curve emerging when connecting two such elementary paths does not enclose and/or intersect a root of P defined by [START_REF] Asheim | Asymptotic analysis of numerical steepest descent with path approximations[END_REF], which is a singularity of 1/ √ P and hence has to be avoided in [START_REF] Bruno | Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case[END_REF]. Hence, this results in two different classes of splitting points where the integration in (8) must be split. Here, explicit formulas for these points are derived according to the parameters q and c 0 in (6)- [START_REF] Bruno | An o (1) integration scheme for three-dimensional surface scattering problems[END_REF]. In Section 4, we also give explicit formulas for the steepest descent path for each case introduced in Section 2. This provides complete means for the computation of (8) thanks to the steepest descent method.

Let us start with the notion of splitting points.

Definition 2.

In the context of Section 2, there exist two splitting points 1. the stationary point c s of g defined by: g (c s ) = 0, 2. and the starting point c r of a path h cr passing through a root of P at some point t cr . Since the real part of g along a steepest descent path is constant, c r is given by

g(h cr (0)) = Re (g(h cr (t cr ))) = Re (g(c)) .
These splitting points can be explicitly computed depending on q (see Eq. ( 7)). A stationary point c s of the function g is solution to g (c s ) = 0. Considering the explicit definition (7) of g, this is equivalent to solving

c 0 (c s -Re(c)) + q P (c s ) = 0. ( 10
)
If q = 0, then we obtain: c s = Re(c). For q = ± √ c 0 , we see that the modulus of c s is infinite.

For infinite values of q, we have c s = Re(c)±Im(c), and then c s is bi-valued as solution to [START_REF] Chew | Fast and Efficient Algorithms in Computational Electromagnetics[END_REF], taking its values in {c, c}. The equation ( 10) is equivalent to: c s = Re(c) ± K cs , setting

K cs := Re(c) 2 + c 0 Re(c) 2 -q 2 |c| 2 c 0 -q 2 .
Injecting the above expression in (10) yields: c 0 K cs = ±q P (c s ). If |q| < √ c 0 , then the left hand side is positive. Thus, depending on the sign of q, we conclude for the exact formula of the stationary point. The case |q| > √ c 0 leads to complex values of c s and can be treated in a similar way.

The starting point c r of a path passing through a root of P is defined as solution to

g(h cr (0)) = Re (g(h cr (t cr ))) ,
which is equivalent to solving the equation

P (c r ) + qc r = Re( P (c) + qc).
Similarly to c s , we derive the exact formulas for c r , depending on q. Let us remark that, if q = 0, the term inside the square-root simplifies to -Im(c) 2 . Therefore, in this case, c r is also bi-valued and takes its values in {c, c}. Finally, the results are summarized in Table 1, setting

K cr := Re(c) 2 + q 2 Re(c) 2 -c 0 |c| 2 c 0 -q 2 . Case c s c r |q| = ∞ {c, c} Re(c) |q| = 0 Re(c) {c, c} |q| = √ c 0 |c s | = ∞ |c r | = ∞ |q| < √ c 0 q < 0 Re(c) + K cs Re(c) + K cr q > 0 Re(c) -K cs Re(c) -K cr |q| > √ c 0 q < 0 Re(c) + K cs Re(c) + K cr q > 0 Re(c) -K cs Re(c) -K cr
Table 1: Explicit formulas of the splitting points according to q.

The previous explicit formulas allow to investigate in details the complementary behaviors of c s and c r . All the possible cases are described in Figure 1. We can conclude that at most one splitting point exists in a given integration interval [a; b].

+i∞ + Re(c) -i∞ + Re(c) -∞ +∞ c s c r c c C +i∞ + Re(c) -i∞ + Re(c) -∞ +∞ c r c s c c C +i∞ + Re(c) -i∞ + Re(c) -∞ +∞ c r c s c c C +i∞ + Re(c) -i∞ + Re(c) -∞ +∞ c s c r c c C Figure 1:
Splitting points of c r and c s according to q. Upper left: q starts at -√ c 0 and tends to -∞. Upper right: q starts at -√ c 0 and tends to 0. Lower left: q starts at √ c 0 and tends to 0. Lower right: q starts at √ c 0 and tends to ∞.

Construction and properties of admissible steepest descent paths 4.1 Construction of admissible steepest descent paths

We now introduce the following definition.

Definition 3. A steepest descent path is called admissible steepest descent path if its starting point belongs to

R \ {c r , c s }.
Let us use (3) to derive the explicit expressions of the admissible steepest descent paths and its derivatives, which is needed in (4). For a path h x starting at x, we set:

K x (t) = P (x) + qx + it,
with t ≥ 0. By the definition [START_REF] Asheim | Asymptotic analysis of numerical steepest descent with path approximations[END_REF] giving P , we obtain

c 0 h x (t) 2 + c 0 |c| 2 -2c 0 Re(c)h x (t) = (K x (t) -qh x (t)) 2 . ( 11
)
Let us analyze the different situations according to q. For q = √ c 0 , (11) leads to

h x (t) = K x (t) 2 -c 0 |c| 2 2 K x (t) √ c 0 -c 0 Re(c) .
Similarly, the case q = -√ c 0 implies that

h x (t) = K x (t) 2 -c 0 |c| 2 -2(K x (t) √ c 0 + c 0 Re(c)) .
Now, let us consider the complementary case |q| = √ c 0 . Eq. ( 11) leads to

h x (t) = K 1 hx ± K 2 hx , (12)
with t ≥ 0, setting

K 1 hx := c 0 Re(c) -qK x (t) c 0 -q 2 , K 2 hx := K x (t) 2 -c 0 |c| 2 c 0 -q 2 - qK x (t) -c 0 Re(c) c 0 -q 2 2 .
Let us assume, without loss of generality, that Re(c) = 0. Indeed, if not, it only consists in a shift of the system. Then, formula [START_REF] Folland | Introduction to Partial Differential Equations[END_REF] becomes

h x (t) = -qK x (t) c 0 -q 2 ± √ c 0 |c 0 -q 2 | Φ(t), (13)
with Φ(t) = µ -t 2 + 2itν, ν = P (x) + xq and µ = ν 2 -|c| 2 (c 0 -q 2 ). Using these notations, it appears that Re(Φ(t)) = µ -t 2 is decreasing and that Im(Φ(t)) = 2tν is increasing with respect to t. We observe that the real part of the first term of ( 13) is constant. Hence, the variations of Re(h x (t)) are given by those of Re( √ Φ). A direct analysis of √ Φ shows that

• |Im( √ Φ)| is increasing, • if |q| < √ c 0 , then Re( √ Φ) is increasing, • if |q| > √ c 0 , Re( √ Φ) decays.
The term Φ(t) writes

Φ(t) = ± 1 2 (|Φ| + Re(Φ)) ± i 1 2 (|Φ| -Re(Φ)).
Amongst the four possible combinations of signs, the main branch of the complex square-root only provides those with a positive real part, i.e. (+, +) and (+, -). In the case |q| > √ c 0 , we know that the contribution of √ Φ to the real part of the path is still decreasing and positive. Let us recall that, for q > 0, all the paths have a positive imaginary part and an increasing real part. To respect these conditions, we must subtract the decaying real part, and then choose the operator -. Similarly, for q < 0, we choose the operator +. The case |q| < √ c 0 can be investigated in a similar way. The situation q = 0 follows exactly the same procedure, depending on the position of the starting point with respect to c s . Finally, the formulas for the derivatives of the paths are directly obtained by differentiating with respect to t. We summarize the results in Table 2 for a path h x starting at x ∈ R, according to q, where

ρ ± := 1 c 0 -q 2   -qK x (t) ± 2K x (t)K x (t) -2 (qK x (t) -c 0 Re(c)) qK x (t) c 0 -q 2 2K 2 hx    and τ ± := K x (t)K x (t) 2 √ c 0 ∓ 2K x (t)K x (t)c 0 Re(c) + K x (t)c 3/2 0 |c| 2 ±2 K x (t) √ c 0 ∓ c 0 Re(c) 2 . Case h x (t) h x (t) q = ± √ c 0 K x (t) 2 -c 0 |c| 2 ±2(K x (t) √ c 0 ∓ c 0 Re(c)) τ ± |q| < √ c 0 x < c s K 1 hx -K 2 hx ρ - x > c s K 1 hx + K 2 hx ρ + |q| > √ c 0 q < 0 K 1 hx + K 2 hx ρ + q > 0 K 1 hx -K 2 hx ρ -
Table 2: Explicit formulas for h x and h x according to q and x.

ODE reformulation and properties of the admissible steepest descent paths

Let us start with the following proposition.

Proposition 1. Let χ := {x ∈ C | g (x) = 0}.
h x (t)g (h x (t)) = i, t > 0, h x (0) = x. ( 14
)
ii) If g is specified as in [START_REF] Bruno | An o (1) integration scheme for three-dimensional surface scattering problems[END_REF] and x ∈ χ, then g is locally not invertible at x and (3) admits some non-unique solutions originating from x.

Proof. Let us start with i). For fixed x ∈ [a, b] and t ∈ [0, ∞[, the initial value problem is immediately derived from (3). Since g is supposed to be analytic on some domain X ⊂ C containing [a, b] and non-constant since we assume that there exists x ∈ [a, b]\χ, χ is a discrete set and F(t, y) → ig (y) -1 is smooth in the open set X \ χ. In particular, its derivative with respect to h x is continuous on [a, b]\χ. Hence, according to the Cauchy-Kovalevskaya theorem [START_REF] Folland | Introduction to Partial Differential Equations[END_REF] there exists some open interval of [0, ∞[ on which the problem admits a unique analytic solution, hence of class C 1 , and continuous at 0. Let us now focus on ii). This statement is given in [START_REF] Huybrechs | Highly oscillatory quadrature[END_REF]. However, it can be directly checked in our particular framework with [START_REF] Gao | A generalization of Filon-Clenshaw-Curtis quadrature for highly oscillatory integrals[END_REF]. If h x (t) / ∈ χ for all t ≥ 0 and using ( 6), ( 7) gives

     h x (t) = i g (h x (t)) = i P (h x (t)) c 0 (h x (t) -Re(c)) + q P (h x (t)) , h x (0) = x. (15)
Here, at some x ∈ χ there are exactly two solutions: one path in the upper complex half plane and another one in the lower half plane.

Solving the nonlinear system (3) amounts to the computation of the inverse function g -1 which might be a difficult task in general (and impossible at x ∈ χ). However, in many cases only a rough approximation of the exact path is needed which might be accomplished by solving [START_REF] Gao | A generalization of Filon-Clenshaw-Curtis quadrature for highly oscillatory integrals[END_REF].

Before formulating some important properties regarding the trace of steepest descent path, let us recall the following notations: Π := {z ∈ Γ c | Im(z) ∈ ]Im(c), Im(c)[ } and Π c = Γ c \ Π, with Γ c := Re(c) + iR. The value of q and the position of the starting point on the real axis are the two fundamental criteria that determine the compatibility between paths to apply in the steepest descent method. This fact is illustrated by Figure 2.

Theorem 1 below provides a rigorous justification of the behavior of the steepest descent paths which is observed in Figure 2. The proof needs the following proposition.

Proposition 2. Two admissible steepest descent paths can never intersect.

Proof. The result is a consequence of the fact that g must have a constant real part when applied along paths of the steepest descent.

We can now prove the following Theorem.

Theorem 1. In the setting introduced in Section 2, we consider an admissible steepest decent path h x and the parameter q. We denote by c sr a real number which equals the (real) splitting point, when it exists, and we define e := Re(h x (0)) -Re(c sr ). Then, the following results hold i) The absolute value of Im(h x (t)) tends to infinity as t → ∞. ii) The trace h * x does not intersect Π c .

(a) q < - √ c 0 (b) q = - √ c 0 (c) - √ c 0 < q < 0 (d) 0 < q < √ c 0 (e) q = √ c 0 (f) √ c 0 < q
iii) Re(h x ) is a monotonous function of t: it is increasing when q ≥ √ c 0 or for x > c s .

Otherwise, it is a decaying function.

iv) If | q |≥ √ c 0 , then sign(Im(h x (t))) = sign(q). Otherwise, we have:

sign(Im(h x (t))) =
sign(e).

v) If | q |= √ c 0 , then there exists T 0 ≥ 0 such that for all t ≥ T 0 , we have sign(q)Re(h x (t)) > sign(q)Re(c).

Otherwise, for all t ≥ T 0 , for some T 0 ≥ 0, we have: sign(e)Re(h x (t)) > sign(e)Re(c).

Proof. Let us prove the different points.

i) The imaginary part of the phase function g applied to h is bounded in domains with bounded imaginary part. Hence, since h is a steepest descent path, (3) holds and by contraposition Im(h(t)) tends to infinity as t → ∞.

ii) Let us show the statement for the case | q |> √ c 0 (see Figure 2 to sketch the idea), the other situations can be treated in a similar way. We consider a path h cr which satisfies Definition 1 but starting at c r . Such a path reaches the root of P of positive imaginary value at time t cr , and then continuously and monotonously describes Π c + := {z ∈ Π c | Im(z) > 0} (this is deduced by the explicit formulations which also show that for all t > t cr , the derivative of such a path is purely imaginary, with positive absolute value). In addition, h cr is not an admissible steepest descent path, but it however satisfies properties used in the proof of Proposition 2. Hence, it splits the upper complex half plane in two open and disjoint domains. This is the reason why the trace h x of some arbitrary steepest descent path starting at x ∈ R \ {c s , c r } never intersects Π c + . iii) The property is a consequence of the explicit formulation of the paths. iv) From ( 15), it turns out that the derivative of a path at its starting point is purely imaginary, and its sign depends on the position of the starting point and the value of q, as stated in iii). Furthermore, Proposition 2 ensures that no steepest descent path crosses the real axis. Hence, the sign of the imaginary value of the steepest descent paths remains constant.

v) It is a direct consequence of points i)-iv), together with Proposition 2.

Integration

Integration in one-dimension

We show in Theorem 2 below that two admissible steepest descent paths starting from points located on the same side of c r and c s can be connected in such a way that the steepest descent method is applicable. Otherwise, this is usually not true. Thus, the application of the steepest descent method on an interval The previous theorem ensures that if the integration in ( 8) is split at c s and c r as depicted in Figure 2, then the resulting closed curve, composed by real interval, admissible steepest descent path and linear connection between the two (for given finite argument t) fulfill the assumption of the CIT. Together with the following result, this proves the applicability of the steepest descent method to (8) provided that c s and c r are considered accordingly. As M can be arbitrary large, this implies that the integral over ϕ M tends to zero as t 1 , t 2 → ∞. Hence, we proved [START_REF] Huybrechs | On the evaluation of highly oscillatory integrals by analytic continuation[END_REF].

Theorem 2. Let a < b be on the same side of a splitting point, and h

Finally, we obtain that I(k, y, a, b) writes

I(k, y, a, b) = b a f (x)e ik(g(x)) dx = ha e ik √ P (z)+zq+s P (z) dz - h b e ik √ P (z)+zq+s P (z) dz = ∞ 0 e ik √ P (ha(t))+ha(t)q+s P (h a (t)) h a (t)dt - ∞ 0 e ik √ P (h b (t))+h b (t)q+s P (h b (t)) h b (t)dt

Integration over triangular surfaces

In this section, we move from the univariate integral in (8) to our initial objective, namely the numerical computation of the two-dimensional surface integral in [START_REF] Antoine | An introduction to operator preconditioning for the fast iterative integral equation solution of time-harmonic scattering problems[END_REF]. In [START_REF] Huybrechs | The construction of cubature rules for multivariate highly oscillatory integrals[END_REF], the authors provide theoretical results that deduce the application of the steepest descent method for multivariate oscillatory integrals for general functions that are analytic in some open complex neighbourhood of the domain of integration under some additional conditions [START_REF] Huybrechs | The construction of cubature rules for multivariate highly oscillatory integrals[END_REF]Theorem 4.5]. Based on a Filon type approach, the authors also derive a general cubature rule for those cases and quantify the corresponding accuracy. This method requires the computation of nested steepest descent path integrals with integrands that consist of a product of a multivariate monomial times a decaying term of the shape e -k p j , where the p j mean the non-negative integration variables, times the derivatives of the path functions. In this way, the non-oscillatory factor in the original integral vanishes in these nested path integrals and is instead taken into account in form of a truncated Taylor series [START_REF] Huybrechs | The construction of cubature rules for multivariate highly oscillatory integrals[END_REF]Section 5].

Here, we propose to handle the two dimensional integral in a different way. The crucial observation is that in our case, the oscillatory term in the outer integration variable of the nested two-dimensional integral can be factored out from the inner integral. As a consequence, we can avoid nested path integrals and the integration of the non-oscillatory part can be treated by means of usual numerical integration methods.

Let us decompose ∆ 1 in n y layers parallel to the x-axis and compute the value of the integral on each layer in the formula below. We assume that none of the layers in ∆ 1 contain a splitting point. Otherwise, the formula has to be split according to well-chosen sub-areas. The dependency of each parameter on the variables x and y is now specified, and we set s := s x + s y . Then, we have

I 1 ∆ 1 (k) =
P (h 1-y (t), y) h 1-y (t)dt   dy. ( 19
)
To each value y in [0, 1] corresponds a unique layer of ∆ 1 . Then, for y ∈ [0, 1] and a ∈ [0, 1-y], we define the path h a,y as the steepest descent path starting at (a, y) (i.e. the starting point a on the layer y). The dependency of a with respect to y is implicit. We remark now that the function

Λ a : y → ∞ 0 e -kt P (h a(y) (t), y) h a(y) (t)dt ( 20
)
is smooth and does not vary significantly. Hence it might be approximated by constant functions on elementary intervals [y j-1 , y j ] composing [0, 1]. This behaviour is illustrated for the configuration

A =    0 0 1 0 0 1    b =    0 0 0    r =    0.3 0.5 0    θ =    1 2 0    (21)
in Figure 3. This is the crucial point in the numerical calculation of I ∆ 1 (k). We now formalize this observation, showing that, under reasonable assumptions, the deviation of Λ in y is controlled by a Lipschitz constant which is independent of any oscillating term in the variable y. For all x, y in C, we define

µ : (x, y) → 1 P (x, y) , which is C ∞ on D µ := {(x, y) ∈ C 2 /P (x, y) / ∈] -∞, 0]}.
It is clear that on every compact set K ⊂ D µ , the function µ is Lipschitz continuous with some fixed constant L := L K > 0. Let us remind that we have some layers free of any splitting point, implying that we have paths along which the integrand is well-defined.

Then, we have the following theorem which confirms that the computational cost of the method only depends on the discretization but not on the wave number k. .

Proof. We consider y varying on the fixed interval [y j-1 , y j ]. Thus, for all t ∈ R and a ∈ [0, 1 -y], the Lipschitz condition above implies

|µ(h a,y (t), y) -µ(h a,ym j (t), y m j )| ≤ L h a,y (t) y - h a,ym j (t) y m j . ( 22
)
The exponential decay related to e -kt , for t ≥ 0, guarantees that there exists a threshold parameter

T 0 > 0 such that Λ a ≈ T 0 0 e -kt P (h a,y (t), y) h a,y (t)dt,
where ≈ formally means that the equality holds up to an a priori numerically negligible term. The continuity of the paths and their derivatives with respect to y can be deduced from the explicit formulations in Table 2. Hence, by compactness arguments, there exist strictly positive values M h and M h such that one obtains max

y∈[y j-1 ,y j ],t∈[0,T 0 ] h a,y (t) -h a,ym j (t) ≤ M h |y j -y j-1 | (23)
and max

y∈[y j-1 ,y j ],t∈[0,T 0 ] h a,y (t) -h a,ym j (t) ≤ M h |y j -y j-1 | . ( 24
)
Let y ∈ [y j 1 , y j ] and t ∈ [0, T 0 ]. Then, ( 22) and ( 23) lead to the sequence of inequalities

|µ(h a,y (t), y) -µ(h a,ym j (t), y m j )| ≤ L h a,y (t) -h a,ym j (t) 2 + y -y m j 2 1 2 ≤ L M 2 h (y j -y j-1 ) 2 + y j -y j-1 2 2 1 2 = L |y j -y j-1 | M 2 h + 1 4 1 2 . ( 25
)
Moreover, for t ∈ [0, ∞[, we have

|µ(h a,y (t), y)h a,y (t) -µ(h a,ym j (t), y m j )h a,ym j (t)| = µ(h a,y (t), y) h a,y (t) -h a,ym j (t) + h a,ym j (t) µ(h a,y (t), y) -µ(h a,ym j (t), y m j ) ≤ |µ(h a,y (t), y)| • h a,y (t) -h a,ym j (t) +|h a,ym j (t)| • µ(h a,y (t), y) -µ(h a,ym j (t), y m j ) . ( 26
)
Combining the above equation with ( 24) and ( 25) implies

|µ(h a,y (t), y)h a,y (t) -µ(h a,ym j (t), y m j )h a,ym j (t)| ≤ max y∈[y j-1 ,y j ],t∈[0,T 0 ] |µ(h a,y (t), y)| • M h |y j -y j-1 | + max y∈[y j-1 ,y j ],t∈[0,T 0 ] |h a,ym j (t)| • L |y j -y j-1 | M 2 h + 1 4 1 2 =: C y j ,y j-1 . |y j -y j-1 | , (27)
with C y j ,y j-1 a positive constant determined by the size of the layer. Finally, for all y ∈ [y j-1 , y j ], we deduce that

Λ a (y) -Λ a (y m j ) ≤ T 0 0 e -kt µ(h a,y (t), y)h a,y (t) -µ(h a,ym j (t), y m j )h a,ym j (t) dt ≤ C y j ,y j-1 |y j -y j-1 | e -kT 0 -1 =: C Λ y j ,y j-1 |y j -y j-1 | .
Assuming that Λ a is constant over each layer leads to an error parameter which is linearly varying with the layers size and is independent of k. We fix y at the middle value y m j , compute Λ a (y m j ) and factor this term out of the integral over [y j-1 , y j ]. For n y layers, we define Î∆ 1 ,ny (k) as the resulting approximation of I 1 ∆ 1 (k). Finally, one gets

Î∆ 1 ,ny (k) = ny j=1 Λ 0 (y m j ) y j y j-1 e ik √ P (0,y)+qyy+s dy -Λ 1-ym j (y m j ) y j y j-1 e ik √ P (1-y,y)+qx(1-y)+qyy+s dy . ( 28
)
The univariate phase functions in the remaining integrals only depend on y and have a similar form as those studied in the previous sections. Therefore, the same procedure applies and then it remains to apply the one-dimensional steepest descent method to each term. Finally, we end up with the evaluation of six steepest descent paths for each y-layer of the unit triangle.

Numerical results

In this section, we present some numerical experiments for the univariate case and triangular surfaces that show the accuracy and efficiency of the approach and formulas derived in the preceding sections.

The univariate case

For the univariate case, we numerically compute the integral in [START_REF] Bruno | Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case[END_REF] for different spatial configurations and wave numbers. We choose the standard Matlab integral function with allowed absolute error of 10 -16 as ground truth which has been deduced from a previous comparison with a brute force method based on the trapezoidal rule applied with an interval partitioning of resolution 10 -6 . In the steepest descent method, the fast decaying integrals over complex paths are computed by means of the Gauss-Laguerre rule with 160 nodes. The configurations are chosen as follows

A =    0 0 2 0 0 2    , b =    0 -1 0    , r = d    cos(α) sin(α) 0    +    0 sin(α) 0    , θ =    1 0 0    , ( 29
)
where α corresponds to 100 uniformly distributed angles from -π to π and d are 9 uniformly distributed distances from 0.2 to 3.8. The purpose of the second term in the determination of r is to bend the circles around [-1, 1] to avoid an intersection with this interval. These configurations lead to surface integrals as in [START_REF] Asheim | A Gaussian quadrature rule for oscillatory integrals on a bounded interval[END_REF], from where we obtain the univariate integrals as in [START_REF] Bruno | Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case[END_REF] and where we set y = 0, a = 0 and b = 1. Let us remind that in accordance with [START_REF] Bruno | Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case[END_REF], the univariate integrals are denoted by I(k, 0, 0, 1) in the cases treated here. Note that the bounds a = 0 and b = 1 in this notation are just a parametrisation and we actually integrate over the interval [-1, 1] in the physical setting. The configuration together with the relative error between the Matlab integral (I M ) and steepest descent integral (I SD ) is displayed in Figure 4 (in log 10 -scale). We took different wave numbers k ∈ {100, 500, 3000, 5000}, and for each we computed the integral for all 100•9 = 900 scenarios described above and measured the CPU time. The experiment covers the three types of possible splitting point configurations, i.e. no splitting point or c r or c s being present on the integration interval, respectively. For wavenumbers k ∈ {100, 500, 3000}, we achieve an average relative error (our method compared to Matlab) of less than 4.6 × 10 -7 with a standard deviation of 1.4 × 10 -11 . The maximum and minimum relative errors are 1.12 × 10 -4 and 3.26 × 10 -16 , respectively. The average absolute error is 2.2 × 10 -8 with a standard deviation of 2.3 × 10 -12 . The maximum and minimum absolute errors are 4.9 × 10 -6 and 7.2 × 10 -18 , respectively. For wave number k = 5000 we obtain the equal performance, except for 8 out of the 900 configurations, where the relative error is about 0.02 and the absolute error is about 3 × 10 -4 . The accuracy for those cases can however be improved to less than 10 -4 relative error by increasing the number of Gauss-Laguerre nodes to 1000. This shows that the method also works in this regime, but also underlines that some configurations require a careful adjustment of parameters. The average CPU time is reported in Figure 5. We observe that the computational cost grows linearly with respect to k for the Matlab integral function while the steepest descent method leads to a k-independent cost. For most configurations, the number of Gauss-Laguerre points could have been significantly reduced from 160 to about 40 with nearly no loss of accuracy. This provides one option to further accelerate the method. To give some explicit values, we fixed α = 0 and d = 0.6, which corresponds to a configuration with a splitting point, more precisely a stationary point in the middle of the interval. The numerical results for the corresponding univariate integrals are given in Table 3 3: of I(k, 0, 0, 1) with settings described above and different wave numbers k.

The case of triangular surfaces

For the surface integral in (5), we consider the same spatial configuration given in [START_REF] Wu | The numerical steepest descent path method for calculating physical optics integrals on smooth conducting quadratic surfaces[END_REF], where now α corresponds to 20 uniformly distributed angles from -π to π while d is given by 9 uniformly sampled distances between 0.2 and 3.8. Due to computational feasibility, we reduced the number of angles compared to the 1-D case, and choose wavenumbers k ≤ 20 to compare with Matlab integral2. The resolution in the y-direction in the implementation of ( 19) is set to 10 -2 and the decaying integrals over the complex paths are evaluated by the Gauss-Laguerre quadrature rule, taking 1000 nodes. The very large number of Gauss-Laguerre nodes is needed to show the performance in all configurations in the low wavenumber regime. At this point, we would like to mention that the computational costs of the overall method could be further reduced by an adaptive choice of parameters. The relative error in log 10 scale for all configurations is reported in Figure 6 and the absolute error is reported in Figure 7. We achieve an average relative error (our method compared to Matlab) of less than 6.4 × 10 -4 with a standard deviation of 4.5 × 10 -11 . The maximum and minimum relative errors are 4.3 × 10 -3 and 5.6 × 10 -6 , respectively. The average absolute error is 1.4 × 10 -7 with a standard deviation of 3.5 × 10 -11 . The maximum and minimum absolute errors are 1.2 × 10 -4 and 1.2 × 10 -8 , respectively.

To prove the computational efficiency, we run the method for all 9 • 20 = 180 scenarios as described above for each wave number k ∈ {100, 500, 3000, 5000}, and measure the CPU time. For this evaluation, the number of Gauss-Laguerre nodes has been reduced to 160, which is justified by the fact that the results deviate by less than 9.2 × 10 -6 and by 1.2 × 10 -7 on average from those computed with 1000 nodes. The CPU time is reported for the twodimensional steepest descent method only since the Matlab integral2 function could not give reliable results for the wave numbers considered here, and a brute force approach took several minutes to compute reliable results in this setting. The corresponding results are given in Figure 8. We again observe the efficiency of the method with respect to k. From an inspection of the Matlab profiler outcome, it turns out that the increasing CPU time for large wavenumbers in Figure 8 is solely due to increasing computational costs of Matlab integration functions needed to handle integration about splitting points. In agreement with the theory, the computational cost of the pure steepest descent method is independent of the wavenumber.

Finally, to give some concrete results, let us fix α = 0 and d = 0.6 and consider the results for k ∈ {5, 10, 20}. Again, the number of Gauss-Laguerre nodes was set to 1000 due to the low wavenumber regime. As comparison, we show integral values in Table [START_REF] Asheim | A remedy for the failure of the numerical steepest descent method on a class of oscillatory integrals[END_REF] computed by the integral2 function of Matlab with the paramater AbsTol being set to 10 -16 . For the computations, we approximate the unit triangle by 10000 rectangular domains, each of which has a height equal to 10 -4 and varying length equal to 1, 1 -10 -4 , 1 -2 × 10 -4 , ..., 10 -4 , respectively, and compute the integral for each of these domains separately with integral2. We see a good agreement between the computed values.

Conclusion

In this paper, we developed the application of the steepest descent method to the evaluation of the three-dimensional high-frequency acoustic single-layer operator over a surface Γ discretized by triangles. We derive formulas for both the splitting points and the admissible steepest descent paths, as well as some of their properties. We theoretically proof that the steepest descent method is applicable in every configuration as long the splitting points are treated accordingly. The integration scheme is next developed and numerical examples are reported to illustrate the efficiency and accuracy of the methodology, most particularly with respect to large wave numbers k. The approach should now be extended to other kinds of highly oscillatory integrals, and most particularly to the various standard surface potentials arising in three-dimensional acoustic and electromagnetic scattering.

Figure 2 :

 2 Figure 2: Steepest descent paths for several starting points in [-2.5, 2.5] for q = -1.5, -1, -0.7, 0.7, 1, 1.5 in (a), (b), (c), (d), (e), (f ), respectively, and throughout c 0 = 1, c = 0.5 + i, c = 0.5 -i. The red dot designates c r , the green one is c s and the black ones refer to c, c.

  [a, b] that contains c s or c r requires a decomposition of the integral in three terms I(k, y, a, b) = I(k, y, a, a 1 ) + I(k, y, a 1 , b 1 ) + I(k, y, b 1 , b), (16) with [a 1 , b 1 ] a small interval that contains either c s or c r . The integral I(k, y, a 1 , b 1 ) over the segment [a 1 , b 1 ] can be computed by standard quadrature rules.

  e kt P (h b (t)) h b (t)dt.

  y,y)+qy(1-y)+qxy+sP (h 1-y (t), y) h 1-y (t)dt P (h 0 (t), y) h 0 (t)dt -e ik √ P (1-y,y)+qy(1-y)+qxy+s ∞ 0 e -kt

Figure 3 :Theorem 3 .

 33 Figure 3: Variation of |Λ a | for the configuration given by (21) with a resolution of 10 -2 in the y-direction in[START_REF] Iserles | On quadrature methods for highly oscillatory integrals and their implementation[END_REF].

Figure 4 :Figure 5 :

 45 Figure 4: Relative error |I M -I SD |/|I M | (in log 10 -scale) for the univariate integral I(1000, 0, 0, 1) comparing our method (denoted by I SD ) with Matlab integral (given by I M ). The colored circle-like contours correspond to the different observation points, as described above, with the corresponding log 10 -relative error values as z-component. The black line segment describes the integration interval [-1, 1].

Figure 6 :

 6 Figure 6: Relative error |I M -I SD |/|I M | (in log 10 -scale) for the surface integral I ∆ (20) comparing our method (denoted by I SD ) with Matlab integral2 (given by I M ). The colored circle-like contours correspond to the different observation points, as described above, with the corresponding log 10 -relative error values as z-component. The black line segment describes the orthogonal projection of the triangle surface onto the xy-plane.

Figure 7 :

 7 Figure 7: Error |I M -I SD | for the surface integral I ∆ (20) comparing our method with Matlab integral2. The colored circle-like contours correspond to the different observation points, as described above, with the corresponding absolute error values as z-component. The black line segment describes the orthogonal projection of the triangle surface onto the xy-plane.

Figure 8 :

 8 Figure 8: Average CPU time of steepest descent method for the surface integral I ∆ (k) for the 180 configurations as described beford and different wavenumbers k = 100, 500, 1000, 3000, 5000. The standard deviation is indicated by error bars

  Proof. From Theorem 1 i), one gets that both Im(h a (t)) and Im(h b (t)) either tend to +∞ or -∞ as t → ∞. Thus, for an arbitrary large M > 0, we can find t 1 , t 2 > t 0 such that Im(h a (t 1 )) = Im(h b (t 2 )) and |Im(h a (t 1 ))| > M . Let ϕ M denotes the path that describes the linear segment connecting h a (t 1 ) and h b (t 2 ) and let φ M be the closed path along [a, b], h b , ϕ M and h a . Theorem 1 i), entails that neither h a nor h b intersect Γ c . This condition and the assumption of Theorem 2 ensure that e ikg(z) √ By construction, ϕ M is a line segment with constant imaginary part that is entirely located in either {z : Re(z) < Re(c)} or {z : Re(z) > Re(c)}. Therefore, one easily verifies that the imaginary part of the phase function g takes its maximum value over ϕ M at h a (t 1 ) or h b (t 2 ).

				P (z)	is holomorphic in an open neighborhood of the
	interior of φ M and hence by the CIT one gets
			φ M	e ikg(z) P (z)	dz = 0.	(18)
					a and h b be two
	admissible steepest descent paths, such that for some t 0 ≥ 0
	min(Re(h a (t)), Re(h b (t))) > Re(c) or max(Re(h a (t)), Re(h b (t))) < Re(c)
	when t > t 0 . Then, the following equality holds
	b a	e ikg(x) P (x)	dx =

∞ 0 e ikg(ha(t)) P (h a (t)) h a (t) dt -∞ 0 e ikg(h b (t)) P (h b (t)) h b (t) dt.

(17)

  .

	k	Steepest Descent	Matlab integral function	Rel. Error
	100	-0.0779921 -i0.1343569	-0.0779921 -i0.1343569 1.44 × 10 -15
	500	0.0481076 -i0.529484	0.0481076 -i0.0529484	1.15 × 10 -14
	1000	-0.0383622 -i0.0337473	-0.0383623 -i0.0337423	9.71 × 10 -5
	3000	-0.0237740 -i0.0180741	-0.0237740 -i0.0180741	5.11 × 10 -7
	5000 -0.01930278 -i0.01243698 -0.0193027 -i0.0124369	1.31 × -8

Table 4 :

 4 Evaluation of I ∆ (k) with the settings described above.

	1
	0.8
	0.6
	0.4
	0.2
	0	100	500	1000	1500	2000	2500	3000	3500	4000	4500	5000