ХИМИЧЕСКАЯ ФИЗИКА, 2004, том 23, № 2, с. 13-19
ЭЛЕМЕНТАРНЫЕ ^^^^^^^^
ФИЗИКО-ХИМИЧЕСКИЕ ПРОЦЕССЫ
УДК 539.1
NUMERICAL EVALUATION OF CUSPOID AND BESSOID OSCILLATING INTEGRALS FOR APPLICATIONS
IN CHEMICAL PHYSICS
© 2004 r. J. N. L. Connor*, C. A. Hobbs**
* Department of Chemistry, University of Manchester, Manchester M13 9PL, UK **Department of Mathematical Sciences, School of Technology, Oxford Brookes University, Wheatley Campus, Oxford OX33 1HX, UK Received 16.11.2002
Oscillating integrals often arise in the theoretical description of phenomena in chemical physics, in particular in atomic and molecular collisions, and in spectroscopy. A computer code for the numerical evaluation of the oscillatory cuspoid canonical integrals and their first-order partial derivatives is described. The code uses a novel adaptive contour algorithm, which chooses a contour in the complex plane that avoids the violent oscillatory and exponential natures of the integrand and modifies its choice as necessary. Applications are made to the swallowtail canonical integral and to a bessoid integral.
1. INTRODUCTION
Oscillating integrals with coalescing saddle points often describe the scattering of atoms and molecules under short wavelength (or high frequency) conditions. They also arise in spectroscopic problems. Typical areas of application include the following:
• Rotational rainbows in atom-molecule and molecule-surface scattering.
• Analysis of vibrational transitions in molecular collisions.
• Theory of chemical reactions.
• Electron detachment in the scattering of negative ions.
• Charge transfer collisions.
• Penning ionisation.
• Pressure broadening of spectral lines and the appearance of satellite bands in spectra.
• Evaluation of Franck-Condon factors and the description of radiative transitions.
• Analysis of quasi-molecular orbital X-ray spectra.
• Atoms in magnetic fields and Stark spectroscopy.
• Autler-Townes effect for pulsed lasers.
More generally, oscillating integrals frequently arise in the theory of water, geophysical, electromagnetic and acoustic waves, as well as in the scattering of heavy nuclear ions. Many references to relevant research can be found in references [1-3].The evaluation of these os-
cillating integrals is often done using uniform asymptotic techniques [1, 2, 4]. Then an important problem is the numerical computation of certain canonical integrals and their first order partial derivatives. The simplest examples involve the cuspoid canonical integrals, which arise in the uniform asymptotic theory of oscillating integrals with 2, 3, 4, ..., n coalescing saddle points [1-4].
For two coalescing saddles, the canonical integral is the regular Airy function (or fold canonical integral):
Ai( x ) = ( 2n)J
exp
-.( 1 3 i I 3u + xu
du.
(1)
For three coalescing saddles, the canonical integral is the Pearcey integral (or cusp canonical integral) [5, 6]:
P(x, y) = J exp[i(u4 + xu2 + yu)]du.
(2)
In the case of four coalescing saddles, the canonical integral is
f 5 3 2
S(x, y, z) = I exp[i(u + xu + yu + zu)]du, (3)
which is called the swallowtail canonical integral [1-4].
Im u
Ci
R
Rn
Re u
Fig. 1. Contours used for the numerical evaluation of I+ (a). For simplicity, the superscript "+" has been omitted from C+ , C+ , C'+ , C+ , R+, R+ and M+.
The purpose of this paper is to describe a computer code for the numerical evaluation of integrals of the type (i)-(4), and to present some illustrative results. The computer code is called CUSPINT and it uses an algorithm in which the integration path along the real axis is replaced by a more convenient contour in the complex u plane, rendering the oscillatory integral more amenable to numerical quadrature.
Only an outline of the algorithm is presented since full details can be found in reference [3]. Numerical results are presented for \S(x, y, z)\. In addition, some results are shown for the bessoid integral [8]
r 4 2
J(x, y) = IJ0(yu)uexp[i(u + xu )]du (5)
where x and y are real and J0(...) is the Bessel function of order zero. The bessoid integral (5) characterizes cuspoid focusing when axial symmetry is present [8].
Section 2 outlines three general methods for the numerical evaluation of cuspoid and bessoid integrals. The adaptive contour algorithm used in CUSPINT is described in section 3. Results for \S(x, y, z)\ and J(x, y) are presented in section 4.
0
In the general case of n _ 1 coalescing saddles, it is necessary to compute numerically the cuspoid canonical integral
C
,(a) = | exp[ifn(a; u)]du,
(4)
a = (ai, a2, ..., an-2),
and its n - 2 first order partial derivatives
3Cn(a) 3Cn(a) dCn(a)
da1 ' da2 ' dan
where
fn( a; u) = u + ^
aku
k = i
with u and the ak real and n an integer greater than 2. The name of Cn(a) arises because fn(a; u) is a miniversal unfolding of a cuspoid singularity (also called an An _ i singularity) [1-4]. The importance of the integrals (i)-(4) can be seen from the fact that the forthcoming Digital Library of Mathematical Functions will contain a chapter entitled "Integrals with Coalescing Saddles" [7].
2. METHODS FOR THE NUMERICAL EVALUATION OF OSCILLATING INTEGRALS
There are three general methods available for the numerical evaluation of infinite oscillating integrals of the type (i)_(5).
a) Maclaurin series. Maclaurin series expansions of (i)_(5) converge for all values of x, y, ... However, they are cumbersome to use when \x\, \y\, ... are large because of cancellation and slow convergence.
b) Differential equations. This method derives a set of differential equations to which Cn(a) or J(x, y) is a solution and then solves the equations numerically. There are several advantages associated with this method: the derivatives are obtained automatically and it is an efficient way of generating grids of values for use in plotting. Disadvantages include: the derivation of the differential equations and their initial conditions is non-trivial; the method is difficult to implement on a computer for the case of general n; for certain values of x, y, . the independent solutions of the differential equations are exponentially increasing, thereby limiting the accuracy to which the integrals can be calculated.
For P(x, y) it is only necessary, in practice, to solve numerically one differential equation [9, 10] and this method has been used by Kaminski and Paris [11] to
2
2
n
30 20
Fig. 2. Grey shaded perspective plot of |S(x, y, z)| for x = 4.0.
30 20
Fig. 3. Grey shaded perspective plot of |S(x, y, z)| for x = 0.0.
study the zeroes of P(x, y) over a wide range of values of x and y (see also [12]). But for Cn(a) with n > 4, the disadvantages of the differential equation method become serious.
c) Contour Integral method. Since the integrand of Cn(a) is infinitely oscillating along the real axis, a direct numerical evaluation is not possible. However, by deforming the contour of integration into the complex u
30 20
Fig. 4. Grey shaded perspective plot of |S(x, y, z)| for x = -6.0.
Modulus 1.5-,
1.0
0.5 -
-10
0
x
-5
0
У
10
10 10
Fig. 5. Perspective and contour plots of |/(x, y)|. The broken curves are the branches of the caustic.
0
5
y
Fig. 6. Contour plot of argJ(x, y)/deg. The contours are -180(30)180. The thick full curves mark the phase discontinuities where argJ(x, y)/deg jumps in value from -180 to +180. The broken curves are the branches of the caustic.
plane we can make the integrand more amenable to a numerical quadrature. This method has the advantage that it is efficient, gives high accuracy results, is relatively easy to implement on a computer and can be generalized to other types of oscillating integrals such as J(x, y).
3. ADAPTIVE CONTOUR METHOD
The first step is to write the general cuspoid integral (4) in the form
Cn(a) = I+ (a) + I-(a), n = 3, 4, 5,.
where
b,( a ) = J exp [ ifn( a; ±u )] du,
a = (ai, a2, ..., an-2)•
(6)
We illustrate the adaptive contour method for In (a), as the procedure for I- (a) is similar [3].
Next we observe that the integrand of equation (6) is infinitely oscillating along the real axis making a direct numerical evaluation impossible. However, we can use the ray from 0 to ^ exp(in/2n) as a new contour of
integration. This does not change the value of I+ (a), as follows from an application of Cauchy's Theorem and Jordan's Lemma, using the fact that the integrand is entire. The new contour has the advantage that the integrand eventually becomes exponentially small, like exp(_tn) with t real, which suggests that a numerical approximation to I+ (a) should be possible.
However, there is still a problem [13]: for certain values of the coefficients, ak, the integrand can possess violent oscillations along the new contour, before it becomes exponentially small. This is a serious difficulty, which can prevent the accurate numerical
evaluation of I+ (a) along the direct ray from 0 to ^ exp(in/2n) [13].
We solve the difficulties discussed above by a compromise in the choice of contours [3]. Figure i shows
the three contours C1+ , C2+ , and C3+ that we employ.
I
0
The contour C+ proceeds from the origin along the real axis to a breakpoint, R+. The second contour C+ is (usually) a straight line, which joins the point R+ to the point R+exp(in/2n). The third contour C+ lies along the original direct ray and joins R+ exp(in/2n) to M+exp(in/2n) with M+ > R+. For suitable choices
of R+ and M+, the infinite integral I+ (a) can evidently be accurately evaluated, provided we can numerically compute the three finite integrals along C+, C+, and C+.
The quadratures along C+, C+, and C+ are performed in CUSPINT using specialist quadrature routines, especially suited to oscillating non-singular integrands. In particular, two versions of the code have been written: the first version uses the subroutine D01AKF present in the NAG Program Library [14], while the second version uses the subroutine DQAG and dependencies from the QUA
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.