Dr. Juan D. González
Acoustic Systems Division – Head –
Acoustic Propagation Department.
Argentinian Navy Research Office
UNIDEF (National Council of Scientific and Technical Research/ Ministry of Defense)
Laprida 555
Vicente Lopez, (1638)
Buenos Aires. Argentina

Acoording to a paper related to the Kirchoff approximation in acoustic scattering, the scattering from a convex object, can be assessed by computing a surface integral over a mesh of \textit{curved} triangles. In that work, integrals are performed by means of the Gauss-Krongrold Rule, in spite of its inefficiency for high oscilatory behavior integrands, we have obtained good results in comparisions with previously published works that belong to the area of computational acoustic.

To accelerate the numerical integration, high oscillatory quadrature rules can be used. This has lead us to study the steepest descent method, in particular, your published work [2].

Our main objective is to compute a double integral which has the following form

$$\int_{0}^{1}\int_{0}^{1-x} \tilde{P}(x,y)e^{i\omega Q(x,y)} dy dx,$$

where $P$,$Q$ are bivariate polynomials of degree 2. By doing a linear change of variables it is possible to re-write the previous integral as $$\int\int_\limits{\triangle_{ABC}} P(x,y)e^{i\omega(\beta_1x^2+\beta_2y^2)} dy dx,$$
where $\triangle_{ABC}$ is a triangle with vertices $A,B,C$ and both integrals differ by a multiplicative constant.

For the moment, suppose $\beta_1=\beta_2=1$.

Without loss of generality we can asumme that the integral will be

$$\int_a^b\int_{\ell_1(x)}^{\ell_2(x)} P(x,y)e^{i\omega(x^2+y^2)} dy dx,$$

$$\int_a^be^{i\omega x^2} \int_{\ell_1(x)}^{\ell_2(x)} P(x,y)e^{i\omega y^2} dy dx$$

The inner integral can be computed following Chapter 5 of reference [1], that leads to

$$\int_a^b e^{i\omega x^2} \left\{e^{i\omega \ell_1^2(x)}F(\ell_1(x),x,\omega)-e^{i\omega \ell_2^2(x)}F(\ell_2(x),x,\omega)+ 1_{\{\ell_1(x)\ell_2(x)<0\}}F_0(x,\omega)\right\} dx,$$
where

Here, we split the integral in three terms, namely

$$I_0 = \int_a^b e^{i\omega x^2} F_0(x,\omega) dx$$

$$I_1 = \int_a^b e^{i\omega (x^2 +\ell_1^2(x))} F(\ell_1(x),x,\omega) dx$$

and $I_2$ analogous to $I_1$.

Writing $(x^2 +\ell_1^2(x)) = \alpha(x-\xi)^2 + C$, we can rewrite $I_1$ as

$$ I_1 = e^{i\omega C} \int_{a^*}^{b^*} e^{i\tilde{\omega}x^2} F(\ell_1(x+\xi),x+\xi,\omega) dx$$

then, it would be possible to evaluate $I_1$ and $I_0$ through the techniques described in [1] (again).

We have implemented this method, however, a number of problems appeared, namely:

* $F(\ell_1(x),x,\omega)$ can be discontinuous if there is a change in the sign of $\ell_1(x)$ on the interval $(a,b)$. This could be fixed, by defining $x^*$ such that $\ell_1({x^*})=0$ and express $I_1$ as
$$I_1 = \int_{a}^{x^*} \dots dx+ \int_{x^*}^{b}\dots dx $$
But in this case, it is not clear for us how to properly define $F(\ell_1(x^*),x^*,\omega)$ because $\ell_1(x^*)=0$.

* We do not know what happen if $\ell_1(a)=\ell_2(a)$, the nature of the application of our interst is such that, both functions start in the same point.

* For the univariate case, we have observed that performance of the method get worse when the lower integral limit is next to the critical point. That is, for example, cases like $\int_{0.01}^b f(x)e^{ix^2}dx$ has poor convergence properties, we have to use 200 Gauss-Legendre quadrature rule points, in orther to get a good approximation.

When these cases do not happen, \footnote{for example \ell_1(x) and \ell_2(x) do not change the sign, and integration limits are far away from the critical point.} the method works properly, moreover, results are impressive, one can achieve relative error of order $10^{-14}$ with only 5 function evaluations. Nevertheless, in the application, we have to take into account the integral for a summatory of thousands of triangles, so that the problematic cases are very likely to happen in our real applicacion.