Abstract:
Propagation of sound in an ocean is accurately predicted at all ranges from a source by a computational method that uses a newly created Decomposition Method to make the predictions at long range from the source. The predictions are done in a limited region of the ocean without the need to predict propagation of sound outside the region. The vertical boundaries of the region are parallel to a given direction of propagation of sound from the source. These boundaries are transparent, and sound penetrates them without reflection. The Decomposition Method is based on a special splitting of the sound field into primary and secondary components. The primary components are related to the sound field that would exist in the region if the vertical boundaries of the region were acoustically hard. They are found by integrating a system of coupled parabolic equations that are derived from the coupled mode equations for the sound field in the region. The primary components contain spurious reverberations that originate at the vertical boundaries of the region. The secondary components identify and cancel the spurious reverberations in the primary components. The secondary components are found by integrating standard parabolic equations subject to transparent boundary conditions at the vertical boundaries of the region. These boundary conditions also include source terms that are obtained from the primary components. When the primary and secondary components are summed, the spurious reverberations are cancelled in the sum. The actual sound field in the region is found from this sum. A modification of the Decomposition Method incorporates an extrapolation procedure to ensure accuracy of the predictions at long range from the source.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
   This application claims the benefit of PPA Ser. No. 60/478,573, filed 2003 Jun. 12 by the present inventor. 
   FEDERALLY SPONSORED RESEARCH 
   Not applicable 
   SEQUENCE LISTING OR PROGRAM 
   Not applicable 

   BACKGROUND OF THE INVENTION 
   1. Field of the Invention 
   The invention pertains to a method for predicting propagation of sound in an ocean. 
   2. Description of the Related Art 
   For the purpose of predicting propagation of sound in an ocean, the ocean can be modelled as a three-dimensional acoustic waveguide that is spatially irregular. A partial differential wave equation governs propagation of sound in the waveguide. A complex-valued distribution of pressure represents the sound field in this equation. However, the sound field also can be represented in terms of the local modes of the waveguide. This representation is useful because it reduces the dimensionality of the model. The equations that govern this representation are obtained by taking moments of the wave equation with respect to the local modes. These moments produce a system of coupled mode equations in two horizontal dimensions. The solution of these coupled mode equations can be used to reconstruct the sound field at any depth in the waveguide from the local mode representation. 
   In  FIG. 1  there is shown a circular region  10  in the horizontal plane with a source of sound at the center  12 . An article by A. T. Abawi, W. A. Kuperman and M. D. Collins (“The coupled mode parabolic equation”, Journal of the Acoustical Society of America, Vol. 102, No. 1, pp. 233-238, July 1997) discusses a numerical method for solving coupled mode equations to predict propagation of sound over region  10 . This method is an elaboration of an earlier numerical method that was developed by M. D. Collins ignoring coupling among the modes. Propagation of sound is predicted over all of region  10 . The preferred direction of propagation is the radially outward direction from the center  12 . 
   An article by M. D. Collins (“The adiabatic mode parabolic equation”, Journal of the Acoustical Society of America, Vol. 94, No. 4, pp. 2269-2278, October 1993) points out that it would be very efficient to make predictions just over an angular sector  14  of small angular width. The article notes that transparent boundary conditions would have to be imposed along the radial edges of sector  14  for this purpose, but it does not indicate how to formulate these conditions or how to adopt them in a solution method. The article by A. T. Abawi, W. A. Kuperman and M. D. Collins does not address these issues at all. 
   SUMMARY OF THE INVENTION 
   The invention permits propagation of sound in an ocean to be accurately predicted at all ranges from a source. At long ranges this prediction is done by a newly created Decomposition Method. The predictions are done for fully three-dimensional ocean environments. 
   The invention permits the predictions to be done in a limited region of the ocean without the need to predict propagation of sound outside the region. The vertical boundaries of the region are parallel to a given direction of propagation of sound from the source. These boundaries are transparent, and sound penetrates them without reflection. 
   The Decomposition Method is based on a special splitting of the sound field into primary and secondary components. The primary components are related to the sound field that would exist in the region if the vertical boundaries of the region were acoustically hard. They are found by integrating a system of coupled parabolic equations that are derived from the coupled mode equations for the sound field in the region. The primary components contain spurious reverberations that originate at the vertical boundaries of the region. 
   The secondary components identify and cancel the spurious reverberations in the primary components. The secondary components are found by integrating standard parabolic equations subject to transparent boundary conditions at the vertical boundaries of the region. These boundary conditions also include source terms that are obtained from the primary components. When the primary and secondary components are summed, the spurious reverberations are cancelled in the sun. The actual sound field in the region is found from this sum. 
   A modification of the Decomposition Method incorporates an extrapolation procedure to ensure that the predictions are accurate at long range from the source. 
   It is an object of the invention to provide an accurate method for predicting propagation of sound in an ocean. 
   Another object of the invention is to provide an accurate method for predicting propagation of sound from a source in a limited region of an ocean without the need to predict propagation of sound outside the region. 
   Yet another object of the invention is to provide an accurate method for predicting propagation of sound from a source in a limited region of an ocean without the need to predict propagation of sound outside the region when the length of the region is much greater than the width of the region. 
   Other objects, advantages and new features of the invention will be apparent from the following detailed description of the invention when considered in conjuction with the accompanying drawings. 

   
     DRAWINGS 
       FIG. 1  shows a circular region where prior art methods predict propagation of sound. 
       FIG. 2  shows a section of a waveguide that models the ocean environment. 
       FIG. 3  shows initial and boundary conditions in the horizontal plane of range and cross-range. 
       FIG. 4  shows a flowchart for the preferred embodiment of the method of the invention. 
   

   DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT 
   1. Coupled mode equations in slab geometry and in spherical geometry. 
   The propagation of sound in an irregular three-dimensional waveguide can be described mathematically by a partial differential equation subject to certain conditions at the boundaries of the waveguide. (By “irregular” we mean that the medium in the waveguide is spatially inhomogeneous. We also assume that this medium can be modelled as a lossless fluid.) Let the variables x, y and z denote range, cross-range and depth in the waveguide, respectively. Let the variables ω, c and ρ denote source frequency, local sound speed and local fluid density in the waveguide, respectively. We let the variable p denote the complex acoustic pressure. This quantity is defined by the property that the physical acoustic pressure is the real part of the product p exp(−iωt), where t denotes time. The partial differential equation in three dimensions that governs p is 
                     ∂               ∂   x       ⁢     (       1   ρ     ⁢       ∂   p       ∂   x         )       +         ∂               ∂   y       ⁢     (       1   ρ     ⁢       ∂   p       ∂   y         )       +         ∂               ∂   z       ⁢     (       1   ρ     ⁢       ∂   p       ∂   z         )       +       1   ρ     ⁢       ω   2       c   2       ⁢   p       =   0.           (   1.1   )             
 
   Referring now to  FIG. 2 , we assume that the pressure vanishes at a reflective top surface  20  of the waveguide, where z=0, and at a reflective bottom surface  22  of the waveguide, where z=−L:
 
p=0 at z=0 and at z=−L.  (1.2) 
 
We also assume that there is a single interface  24  between the surface and the bottom where the local sound speed and the local fluid density can be discontinuous. The depth of this interface can vary with range and with cross-range. We describe the interface by the equation
 
 S ( x,y,z )= H ( x,y )+ z= 0,  (1.3) 
 
where 0&lt;H(x, y)&lt;L. The acoustic pressure must be continuous across this interface. The component of acoustic fluid velocity that is normal to the interface must also be continuous across the interface. We can express these continuity conditions as follows. 
       p   ⁢           ⁢   and   ⁢           ⁢     1   ρ     ⁢     (           ∂   H       ∂   x       ⁢       ∂   p       ∂   x         +         ∂   H       ∂   y       ⁢       ∂   p       ∂   y         +       ∂   p       ∂   z         )         
 are continuous across z=−H(x, y).  (1.4)
 
   It is possible to convert the acoustic equation (1.1) and the conditions (1.2) and (1.4) on it to an equivalent infinite system of coupled partial differential equations in two dimensions. This system is based on the local modes of the waveguide, which are characterized by the eigenfunctions ƒ n (x, y, z) and by the squared wavenumbers ξ n   2 (X, y) that satisfy the differential equation 
               ∂               ∂   z       ⁢     (       1   ρ     ⁢       ∂     f   n         ∂   z         )       +       1   ρ     ⁢     (         ω   2       c   2       -     ξ   n   2       )     ⁢     f   n         =   0       
 if −L&lt;z&lt;0 and z≠−H(x,y),  (1.5)
 
subject to the conditions that
 
     ƒ n =0 at z=−L and at z=0,  (1.6)          f   n     ⁢           ⁢   and   ⁢           ⁢     1   ρ     ⁢       ∂     f   n         ∂   z           
 are continuous across z=−H(x, y).  ( 1.7 ) 
   At each position (x,y) these equations describe a nonsingular selfadjoint characteristic value problem. The squared wavenumbers ξ n   2 (x, y) are real and countably infinite in number. They can be indexed to decrease monotonically to −∞ as the index n increases through the positive integers. The local modes can be computed with numerical methods that are well-known in the art. 
   At each position (x, y), σƒ n (x, y 0)/σz=0 if and only if ƒ n (x, y, z) vanishes identically for −L≦z≦0. Therefore, we can assume that σƒ n (x, z, y, 0)/σz≠0 and that ƒ n  (x, y ,z) does not vanish identically. Furthermore, we can assume that ƒ n (x, y, z) is real-valued. Since the density ρ is positive, the eigenfunctions are orthogonal in the following sense: 
                 ∫     -   L     0     ⁢       1   ρ     ⁢     f   m     ⁢     f   n     ⁢     ⅆ   z         =     {                   d   n     ⁡     (     x   ,   y     )       &gt;     0   ⁢           ⁢   if   ⁢           ⁢   m       =   n     ,                 0   ⁢           ⁢   if   ⁢           ⁢   m     ≠     n   .                       (   1.8   )             
 
   We also assume that the eigenfunctions are normalized so that the signs of the partial derivatives σƒ n (x, y, 0)/θz and the values of the integrals d n (x, y) are independent of the position (x, y). 
   The functions p, σp/σx and σp/σy have the following eigenfunction expansions: 
               p   =         ∑     n   =   1     ∞     ⁢       a   n     ⁢       f   n     /     d   n       ⁢           ⁢   where   ⁢           ⁢     a   n         =       ∫     -   L     0     ⁢       1   ρ     ⁢     pf   n     ⁢     ⅆ   z             ,           (   1.9   )                     ∂   p       ∂   x       =         ∑     n   =   1     ∞     ⁢       b   n     ⁢       f   n     /     d   n       ⁢           ⁢   where   ⁢           ⁢     b   n         =       ∫     -   L     0     ⁢       1   ρ     ⁢       ∂   p       ∂   x       ⁢     f   n     ⁢     ⅆ   z             ,           (   1.10   )                   ∂   p       ∂   y       =         ∑     n   =   1     ∞     ⁢       c   n     ⁢       f   n     /     d   n       ⁢           ⁢   where   ⁢           ⁢     c   n         =       ∫     -   L     0     ⁢       1   ρ     ⁢       ∂   p       ∂   y       ⁢     f   n     ⁢       ⅆ   z     .                   (   1.11   )             
 
The coefficients a n , b n  and c n  in these expansions are functions of the position (x,y). They are governed by the following system of coupled partial differential equations. 
                     ∂     b   n         ∂   x       +       ∂     c   n         ∂   y       +       ξ   n   2     ⁢     a   n         =       ∑     m   =   1     ∞     ⁢       (         b   m     ⁢     B   mn       +       c   m     ⁢     C   mn         )     ⁢     1     d   m             ,           (   1.12   )                   b   n     =         ∂     a   n         ∂   x       +       ∑     m   =   1     ∞     ⁢       a   m     ⁢     B   nm     ⁢     1     d   m               ,           (   1.13   )                 c   n     =         ∂     a   n         ∂   y       +       ∑     m   =   1     ∞     ⁢       a   m     ⁢     C   nm     ⁢       1     d   m       .                   (   1.14   )             
 
These equations constitute the fundamental first-order formulation of coupled mode theory for an irregular three-dimensional waveguide.
 
   The coupling coefficients B mn  and C mn  can be computed as follows. Let 
               g   n     =       1   ρ     ⁢       ∂     f   n         ∂   z                 (   1.15   )             
 
and define the following matrix elements: 
                 D   mn     =       ∫     -   L     0     ⁢       [           ∂               ∂   x       ⁢     (       1   ρ     ⁢       ω   2       c   2         )     ⁢     f   m     ⁢     f   n       +         ∂   ρ       ∂   x       ⁢     g   m     ⁢     g   n         ]     ⁢     ⅆ   z           ,           (   1.16   )                   E   mn     =       ∫     -   L     0     ⁢         ∂               ∂   x       ⁢     (     1   ρ     )     ⁢     f   m     ⁢     f   n     ⁢     ⅆ   z           ,           (   1.17   )                   F   mn     =       ∫     -   L     0     ⁢       [           ∂               ∂   y       ⁢     (       1   ρ     ⁢       ω   2       c   2         )     ⁢     f   m     ⁢     f   n       +         ∂   ρ       ∂   y       ⁢     g   m     ⁢     g   n         ]     ⁢     ⅆ   z           ,           (   1.18   )                   G   mn     =       ∫     -   L     0     ⁢         ∂               ∂   y       ⁢     (     1   ρ     )     ⁢     f   m     ⁢     f   n     ⁢     ⅆ   z           ,           (   1.19   )                   J   mn     =           [       1   ρ     ⁢       ω   2       c   2         ]         -   H     -   0         -   H     +   0       ⁢       f   m     ⁡     (     -   H     )       ⁢       f   n     ⁡     (     -   H     )         +         [   ρ   ]         -   H     -   0         -   H     +   0       ⁢       g   m     ⁡     (     -   H     )       ⁢       g   n     ⁡     (     -   H     )             ,           (   1.20   )                 K   mn     =         [     1   ρ     ]         -   H     -   0         -   H     +   0       ⁢       f   m     ⁡     (     -   H     )       ⁢         f   n     ⁡     (     -   H     )       .               (   1.21   )             
 
If m≠n then 
                 B   mn     =     -         (       ξ   m   2     -     ξ   n   2       )       -   1       ⁡     [       D   mn     -       E   mn     ⁢     ξ   n   2       +         ∂   H       ∂   x       ⁢     (       J   mn     -       K   mn     ⁢     ξ   n   2         )         ]           ,           (   1.22   )                   C   mn     =     -         (       ξ   m   2     -     ξ   n   2       )       -   1       ⁡     [       F   mn     -       G   mn     ⁢     ξ   n   2       +         ∂   H       ∂   y       ⁢     (       J   mn     -       K   mn     ⁢     ξ   n   2         )         ]           ,           (   1.23   )             
 
and 
                 B   nn     =       -     1   2       ⁢     (       E   nn     +         ∂   H       ∂   x       ⁢     K   nn         )         ,           (   1.24   )                 C   nn     =       -     1   2       ⁢       (       G   nn     +         ∂   H       ∂   y       ⁢     K   nn         )     .               (   1.25   )             
 
   In practice, it is necessary to work with a finite number of coupled mode equations. For this reason, we truncate the coupled mode equations (1.12-1.14) as follows. Fixing a positive integer M, we assume that the expansion coefficients a n , b n  and c n  vanish if n&gt;M. Thus, if n=1, . . . , M then 
                     ∂     b   n         ∂   x       +       ∂     c   n         ∂   y       +       ξ   n   2     ⁢     a   n         =       ∑     m   =   1     M     ⁢       (         b   m     ⁢     B   mn       +       c   m     ⁢     C   mn         )     ⁢     1     d   m             ,           (   1.26   )                   b   n     =         ∂     a   n         ∂   x       +       ∑     m   =   1     M     ⁢       a   m     ⁢     B     n   ⁢           ⁢   m       ⁢     1     d   m               ,           (   1.27   )                 c   n     =         ∂     a   n         ∂   y       +       ∑     m   =   1     M     ⁢       a   m     ⁢     C     n   ⁢           ⁢   m       ⁢       1     d   m       .                   (   1.28   )             
 
   We can reduce this system of equations to an equivalent system of second-order partial differential equations for just the coefficients a n . If n=1, . . . , M then 
                       ∂   2     ⁢     a   n         ∂     x   2         +         ∂   2     ⁢     a   n         ∂     y   2         +       ξ   n   2     ⁢     a   n         =       ∑     m   =   1     M     ⁢       (         α     n   ⁢           ⁢   m       ⁢     a   m       +       β     n   ⁢           ⁢   m       ⁢       ∂     a   m         ∂   x         +       γ     n   ⁢           ⁢   m       ⁢       ∂     a   m         ∂   y           )     ⁢     1     d   m             ,           (   1.29   )             
 
where 
                 a     n   ⁢           ⁢   m       =         ∑     l   =   1     M     ⁢       (         B     l   ⁢           ⁢   m       ⁢     B     l   ⁢           ⁢   n         +       C     l   ⁢           ⁢   m       ⁢     C     l   ⁢           ⁢   n           )     ⁢     1     d   l           -       ∂     B     n   ⁢           ⁢   m           ∂   x       -       ∂     C     n   ⁢           ⁢   m           ∂   y           ,           (   1.30   )             
 
   The waveguide considered in the preceding discussion is a slab between two parallel planes. Now let us consider a waveguide that is a shell between two concentric spheres. Let the radius of the outer sphere be R 0 , and let the radius of the inner sphere be R 0 −L. We pick a spherical coordinate system for this problem anticipating that our analysis of sound propagating in the shell will be restricted to a neighborhood of a great circle E on the outer sphere. We identify E as the equator of this sphere. Let M denote another great circle on the outer sphere that intersects E at right angles and near the source of sound in the waveguide. We identify M as the prime meridian of the outer sphere. 
   Let P denote a point in the shell. We let τ be the radial distance of P from the common center of the spheres. We let ψ be the angle of latitude of P relative to the equator E. We let φ be the angle of longitude of P relative to the prime meridian M. For definiteness let φ=0 where E intersects M near the source. 
   It is convenient to introduce a new system of coordinates (x, y, z).
 
x=R 0 φ, y=R 0 ln(sec ψ+tan ψ), z=r−R 0 .  (1.32) 
 
Note that y≈R 0 ψ if |ψ|&lt;&lt;1. In the new coordinate system the partial differential equation that governs the acoustic pressure p in the waveguide is 
                     ∂               ∂   x       ⁢     (       1   ρ     ⁢       ∂   p       ∂   x         )       +         ∂               ∂   y       ⁢     (       1   ρ     ⁢       ∂   p       ∂   y         )       +       cos   2     ⁢   ψ   ⁢       ∂               ∂   z       ⁢     (       1   ρ     ⁢       r   2       R   0   2       ⁢       ∂   p       ∂   x         )       +       cos   2     ⁢   ψ   ⁢     1   ρ     ⁢       ω   2       c   2       ⁢       r   2       R   0   2       ⁢   p       =   0.           (   1.33   )             
 
   We assume that the pressure vanishes on the outer sphere, where z=0, and on the inner sphere, where z=−L:
 
p=0 at z=0 and at z=−L.  (1.34) 
 
We also assume that there is a single interface between these spheres where the local sound speed and the local fluid density in the shell can be discontinuous. The radial coordinate of this interface can vary with latitude and with longitude. We describe the interface in the new coordinate system by the equation
 
 S ( x,y,z )= H ( x,y )+ z= 0,  (1.35) 
 
where 0&lt;H(x, y)&lt;L. We can express the continuity conditions at the interface as follows. 
               p   ⁢           ⁢   and   ⁢           ⁢     1   ρ     ⁢     (           ∂   H       ∂   x       ⁢       ∂   p       ∂   x         +         ∂   H       ∂   y       ⁢       ∂   p       ∂   y         +       cos   2     ⁢   ψ   ⁢       r   2       R   0   2       ⁢       ∂   p       ∂   z           )       ⁢     
     ⁢       are   ⁢           ⁢   continuous   ⁢           ⁢   across   ⁢           ⁢   z     =     -       H   ⁡     (     x   ,   y     )       .                 (   1.36   )             
 
   The local modes of the spherical shell waveguide are characterized by the eigenfunctions ƒ n (x, y, z) and by the squared wavenumbers ξ n   2 (x, y) that satisfy the differential equation 
                     ∂     ∂   z       ⁢     (       1   ρ     ⁢       r   2       R   0   2       ⁢       ∂     f   n         ∂   z         )       +       1   ρ     ⁢     (           ω   2       c   2       ⁢       r   2       R   0   2         -     ξ   n   2       )     ⁢     f   n         =   0     ⁢     
     ⁢         if   ⁢           -   L     &lt;   z   &lt;     0   ⁢           ⁢   and   ⁢           ⁢   z     ≠     -     H   ⁡     (     x   ,   y     )           ,             (   1.37   )             
 
subject to the conditions that
 
ƒ n =0 at z=−L and at z=0,  (1.38) 
 
                 f   n     ⁢           ⁢   and   ⁢           ⁢     1   ρ     ⁢       ∂     f   n         ∂   z       ⁢           ⁢   are   ⁢           ⁢   continuous   ⁢           ⁢   across   ⁢           ⁢   z     =     -       H   ⁡     (     x   ,   y     )       .               (   1.39   )             
 
At each position (x, y) these equations describe a nonsingular selfadjoint characteristic value problem. The squared wavenumbers ξ n   2 (x, y) are real and countably infinite in number. We can assume that they are indexed to decrease monotonically to −∞ as the index n increases through the positive integers. We can also assume that σƒ n (x, y, 0)/σz≠0 and that ƒ n (x, y, z) does not vanish identically. Furthermore, we can assume that ƒ n (x, y, z) is real-valued. Since the density ρ is positive, the eigenfunctions are orthogonal in the following sense: 
                 ∫     -   L     0     ⁢       1   ρ     ⁢     f   m     ⁢     f   n     ⁢     ⅆ   z         =     {               d   n     ⁡     (     x   ,   y     )       &gt;   0               if   ⁢           ⁢   m     =   n     ,             0           if   ⁢           ⁢   m     ≠     n   .                       (   1.40   )             
 
We assume that the eigenfunctions are normalized so that the signs of the partial derivatives σƒ n (x, y, 0)/σz and the values of the integrals d n (x, y) are independent of the position (x, y).
 
   The functions p, σp/σx and σp/σy have the following eigenfunction expansions: 
               p   =         ∑     n   =   1     ∞     ⁢       a   n     ⁢       f   n     /     d   n       ⁢           ⁢   where   ⁢           ⁢     a   n         =       ∫     -   L     0     ⁢       1   ρ     ⁢   p   ⁢           ⁢     f   n     ⁢     ⅆ   z             ,           (   1.41   )                     ∂   p       ∂   x       =         ∑     n   =   1     ∞     ⁢       b   n     ⁢     f   n     ⁢     d   n     ⁢           ⁢   where   ⁢           ⁢     b   n         =       ∫     -   L     0     ⁢       1   ρ     ⁢       ∂   p       ∂   x       ⁢     f   n     ⁢     ⅆ   z             ,           (   1.42   )                   ∂   p       ∂   y       =         ∑     n   =   1     ∞     ⁢       c   n     ⁢       f   n     /     d   n       ⁢           ⁢   where   ⁢           ⁢     c   n         =       ∫     -   L     0     ⁢       1   ρ     ⁢       ∂   p       ∂   y       ⁢     f   n     ⁢       ⅆ   z     .                   (   1.43   )             
 
The coefficients a n , b n  and c n  in these expansions are functions of the position (x,y). They are governed by the following system of partial differential equations. 
                     ∂     b   n         ∂   x       +       ∂     c   n         ∂   y       +       ξ   n   2     ⁢     cos   2     ⁢   ψ   ⁢           ⁢     a   n         =       ∑     m   =   1     ∞     ⁢       (         b   m     ⁢     B     m   ⁢           ⁢   n         +       c   m     ⁢     C     m   ⁢           ⁢   n           )     ⁢     1     d   m             ,           (   1.44   )                   b   n     =         ∂     a   n         ∂   x       +       ∑     m   =   1     ∞     ⁢       a   m     ⁢     B     m   ⁢           ⁢   n       ⁢     1     d   m               ,           (   1.45   )                 c   n     =         ∂     a   n         ∂   y       +       ∑     m   =   1     ∞     ⁢       a   m     ⁢     C     n   ⁢           ⁢   m       ⁢       1     d   m       .                   (   1.46   )             
 
These equations constitute the fundamental order formulation of coupled mode theory for an irregular three-dimensional waveguide in spherical geometry.
 
   The coupling coefficients B mn  and C mn  can be computed as follows. Let 
               g   n     =       1   ρ     ⁢       ∂     f   n         ∂   z                 (   1.47   )             
 
and define the following matrix elements: 
                     ⁢         D     m   ⁢           ⁢   n       =       ∫     -   L     0     ⁢       [         ∂     ∂   x       ⁢     (       1   ρ     ⁢       w   2       c   2         )     ⁢     f   m     ⁢     f   n       +         ∂   ρ       ∂   x       ⁢     g   m     ⁢     g   n         ]     ⁢       r   2       R   0   2       ⁢     ⅆ   z           ,             (   1.48   )                       ⁢         E     m   ⁢           ⁢   n       =       ∫     -   L     0     ⁢       ∂     ∂   x       ⁢     (     1   ρ     )     ⁢     f   m     ⁢     f   n     ⁢     ⅆ   z           ,             (   1.49   )                       ⁢         F     m   ⁢           ⁢   n       =       ∫     -   L     0     ⁢       [         ∂     ∂   x       ⁢     (       1   ρ     ⁢       ω   2       c   2         )     ⁢     f   m     ⁢     f   n       +         ∂   ρ       ∂   y       ⁢     g   m     ⁢     g   n         ]     ⁢       r   2       R   0   2       ⁢     ⅆ   z           ,             (   1.50   )                       ⁢         G     m   ⁢           ⁢   n       =       ∫     -   L     0     ⁢       ∂     ∂   y       ⁢     (     1   ρ     )     ⁢     f   m     ⁢     f   n     ⁢     ⅆ   z           ,             (   1.51   )                   J     m   ⁢           ⁢   n       =       {           [       1   ρ     ⁢       ω   2       c   2         ]         -   H     -   0         -   H     +   0       ⁢       f   m     ⁡     (     -   H     )       ⁢       f   n     ⁡     (     -   H     )         +         [   ρ   ]         -   H     -   0         -   H     +   0       ⁢       g   m     ⁡     (     -   H     )       ⁢       g   m     ⁡     (     -   H     )           }     ⁢           ⁢       (     1   -     H   /     R   0         )     2         ,           (   1.52   )                       ⁢       K     m   ⁢           ⁢   n       =         [     1   ρ     ]         -   H     -   0         -   H     +   0       ⁢       f   m     ⁡     (     -   H     )       ⁢         f   n     ⁡     (     -   H     )       .                 (   1.53   )             
 
These definitions differ in form from the corresponding definitions for the slab waveguide only by the radial factors in the formulas for the matrix elements D mn , F mn  and J mn .
 
   If m≠n then 
                 B     m   ⁢           ⁢   n       =     -         (       ξ   m   2     -     ξ   n   2       )       -   1       ⁡     [       D     m   ⁢           ⁢   n       -       E     m   ⁢           ⁢   n       ⁢     ξ   n   2       +         ∂   H       ∂   x       ⁢     (       J     m   ⁢           ⁢   n       -       K     m   ⁢           ⁢   n       ⁢     ξ   n   2         )         ]           ,           (   1.54   )                   C     m   ⁢           ⁢   n       =     -         (       ξ   m   2     -     ξ   n   2       )       -   1       ⁡     [       F     m   ⁢           ⁢   n       -       G     m   ⁢           ⁢   n       ⁢     ξ   n   2       +         ∂   H       ∂   y       ⁢     (       J     m   ⁢           ⁢   n       -       K     m   ⁢           ⁢   n       ⁢     ξ   n   2         )         ]           ,           (   1.55   )             
 
and 
                 B     n   ⁢           ⁢   n       =       -     1   2       ⁢     (       E     n   ⁢           ⁢   n       +         ∂   H       ∂   x       ⁢     K     n   ⁢           ⁢   n           )         ,           (   1.56   )                 C     n   ⁢           ⁢   n       =       -     1   2       ⁢       (       G     n   ⁢           ⁢   n       +         ∂   H       ∂   y       ⁢     K     n   ⁢           ⁢   n           )     .               (   1.57   )             
 
These equations are formally identical to the corresponding equations for the slab waveguide.
 
   Let us truncate eqs.(1.44-1.46) by fixing a positive integer M and assuming that the expansion coefficients a n , b n  and c n  vanish if n&gt;M. Thus, if n=1, . . . , M then 
                     ∂     b   n         ∂   x       +       ∂     c   n         ∂   y       +       ξ   n   2     ⁢     cos   2     ⁢   ψ   ⁢           ⁢     a   n         =       ∑     m   =   1     M     ⁢       (         b   m     ⁢     B     m   ⁢           ⁢   n         +       c   m     ⁢     C     m   ⁢           ⁢   n           )     ⁢     1     d   m             ,           (   1.58   )                   b   n     =         ∂     a   n         ∂   x       +       ∑     m   =   1     M     ⁢       a   m     ⁢     B     n   ⁢           ⁢   m       ⁢     1     d   m               ,           (   1.59   )                 c   n     =         ∂     a   n         ∂   y       +       ∑     m   =   1     M     ⁢       a   m     ⁢     C     n   ⁢           ⁢   m       ⁢       1     d   m       .                   (   1.60   )             
 
We can reduce this system of equations to an equivalent system of second-order partial differential equations for just the coefficients a n . If n=1, . . . , M then 
                       ∂   2     ⁢     a   n         ∂     x   2         +         ∂   2     ⁢     a   n         ∂     y   2         +       ξ   n   2     ⁢     cos   2     ⁢   ψ   ⁢           ⁢     a   n         =       ∑     m   =   1     M     ⁢       (         α     n   ⁢           ⁢   m       ⁢     a   m       +       β     n   ⁢           ⁢   m       ⁢       ∂     a   m         ∂   x         +       γ     n   ⁢           ⁢   m       ⁢       ∂     a   m         ∂   y           )     ⁢     1     d   m             ,           (   1.61   )             
 
where 
                 α     n   ⁢           ⁢   m       =         ∑     l   =   1     M     ⁢       (         B     l   ⁢           ⁢   m       ⁢     B     l   ⁢           ⁢   n         +       C     l   ⁢           ⁢   m       ⁢     C     l   ⁢           ⁢   n           )     ⁢     1     d   l           -       ∂     B     n   ⁢           ⁢   m           ∂   x       -       ∂     C     n   ⁢           ⁢   m           ∂   y           ,           (   1.62   )             
 
2. Paraxial approximations for coupled mode equations.
 
   Since the coupled mode equations (1.29,1.61) are elliptic partial differential equations, they can model propagation in all directions at the same time. In many problems of practical interest, however, it is sufficient to consider a single direction of propagation at a time. In slab geometry we can take this direction to be the direction of increasing range, and in spherical geometry we can take it to be the direction of increasing longitude along the equator. In each case, we can formulate parabolic equations from the coupled mode equations to model propagation of sound in the given direction. We shall describe how to formulate these parabolic equations from the coupled mode equations in slab geometry. The approach is identical for the coupled mode equations in spherical geometry. 
   We let the normalization integral d m =1 for every index m. Hence, the coupled mode equations take the form 
                       ∂   2     ⁢     a   m         ∂     x   2         +         ∂   2     ⁢     a   m         ∂     y   2         +       ξ   m   2     ⁢     a   m         =       ∑     n   =   1     M     ⁢     (         α     n   ⁢           ⁢   m       ⁢     a   n       +       β     m   ⁢           ⁢   n       ⁢       ∂     a   n         ∂   x         +       γ     m   ⁢           ⁢   n       ⁢       ∂     a   n         ∂   y           )         ⁢     
     ⁢       (       m   =   1     ,   …   ⁢           ,   M     )     .             (   2.1   )             
 
We introduce a positive reference wavenumber k m (x) for each index m at every range x. We define the corresponding reference phase σ m (x) by 
                 ϑ   m     ⁡     (   x   )       =       ∫   0   x     ⁢         k   m     ⁡     (   u   )       ⁢       ⅆ   u     .                 (   2.2   )             
 
We define the complex amplitude φ m (x, y) of the expansion coefficient a m (x, y) as follows.
 
φ m ( x,y )=a m ( x,y ) exp{− iσ   m (x)}.  (2.3) 
 
Noting that dσ m (x)/dx=k m (x) and letting k m (x)=dk m (x)/dx, we see that 
                   ∂     a   m         ∂   x       =       (         ∂     φ   m         ∂   x       +     ⅈ   ⁢           ⁢     k   m     ⁢     φ   m         )     ⁢     ⅇ     ⅈ   ⁢           ⁢   ϑ   ⁢           ⁢   m           ,           (   2.4   )                     ∂   2     ⁢     a   m         ∂     x   2         =       (           ∂   2     ⁢     φ   m         ∂     x   2         +     2   ⁢           ⁢   ⅈ   ⁢           ⁢     k   m     ⁢       ∂     φ   m         ∂   x         +     ⅈ   ⁢                 ⁢   k     .     m     ⁢     φ   m       -       k   m   2     ⁢     φ   m         )     ⁢       ⅇ     ⅈ   ⁢           ⁢   ϑ   ⁢           ⁢   m       .               (   2.5   )             
 
Now we assume that the following paraxial approximations are valid: 
                        ∂   2     ⁢     φ   m         ∂     x   2              ⪡     2   ⁢           ⁢     k   m     ⁢            ∂     φ   m         ∂   x            ⁢           ⁢       (       m   =   1     ,   …   ⁢           ,   M     )     .               (   2.6   )             
 
These assumptions allow us to approximate eq.(2.5) by deleting the first term on its right-hand side (RHS). 
                     ∂   2     ⁢     a   m         ∂     x   2         =       (       2   ⁢   ⅈ   ⁢           ⁢     k   m     ⁢       ∂     φ   m         ∂   x         +     ⅈ   ⁢           ⁢       k   .     m     ⁢     φ   m       -       k   m   2     ⁢     φ   m         )     ⁢     ⅇ     ⅈ   ⁢           ⁢   ϑ   ⁢           ⁢   m           ,           (   2.7   )             
 
Combining eqs.(2.3,2.4,2.7) with eq.(2.1), we get the following coupled parabolic equations. 
                 2   ⁢           ⁢   ⅈ   ⁢           ⁢     k   m     ⁢       ∂     φ   m         ∂     x   2           +     ⅈ   ⁢           ⁢       k   .     m     ⁢     φ   m       +         ∂   2     ⁢     φ   m         ∂     y   2           =         (       k   m   2     -     ξ   m   2       )     ⁢     φ   m       +       ∑     n   =   1     M     ⁢       [         (       α     m   ⁢           ⁢   n       +       ⅈβ     m   ⁢           ⁢   n       ⁢     k   n         )     ⁢     φ   n       +       β     m   ⁢           ⁢   n       ⁢       ∂     φ   n         ∂   x         +       γ     m   ⁢           ⁢   n       ⁢       ∂     φ   n         ∂   y           ]     ⁢     ⅇ     ⅈ   ⁡     (       ϑ   m     -     ϑ   n       )                       (   2.8   )             
 
for m=1, . . . , M.
 
   Note the appearance of the partial derivatives σφ n /σx on the RHS of eq.(2.8). Let us rewrite the coupled parabolic equations as follows. 
               2   ⁢           ⁢     ⅈ   ⁡     [         k   m     ⁢       ∂     φ   m         ∂   x         +       i   2     ⁢       ∑     n   =   1     M     ⁢       β     m   ⁢           ⁢   n       ⁢     ⅇ       -   ⅈ     ⁢           ⁢     (       ϑ   m     -     ϑ   n       )         ⁢       ∂     φ   n         ∂   x               ]         =         -   ⅈ     ⁢           ⁢     k   m     ⁢     φ   m       -         ∂   2     ⁢     φ   m         ∂     y   2         +       (       k   m   2     -     ξ   m   2       )     ⁢     φ   m       +       ∑     n   =   1     M     ⁢       [         (       α     m   ⁢           ⁢   n       +     ⅈ   ⁢           ⁢     β     m   ⁢           ⁢   n       ⁢     k   n         )     ⁢     φ   n       +       γ     m   ⁢           ⁢   n       ⁢       ∂     φ   n         ∂   y           ]     ⁢     ⅇ     -     ⅈ   ⁡     (       ϑ   m     -     ϑ   n       )                         (   2.9   )             
 
for m=1, . . . , M. This system of equations can be solved uniquely for the partial derivatives σφ n /σx if and only if the matrix 
       M   =           ⁢       (   2.10   )     ⁢     
     [           ⁢           k   1             (     i   /   2     )     ⁢     β   12     ⁢     e     -     i   ⁡     (       ϑ   1     -     ϑ   2       )                 …           (     i   /   2     )     ⁢     β     1   ⁢   M       ⁢     e     -     i   ⁡     (       ϑ   1     -     ϑ   M       )                         (     i   /   2     )     ⁢     β   21     ⁢     e     -     i   ⁡     (       ϑ   2     -     ϑ   1       )                   k   2         …           (     i   /   2     )     ⁢     β     2   ⁢           ⁢   M       ⁢     e     i   ⁡     (       ϑ   2     -     ϑ   M       )                   ⋮       ⋮       ⋰       ⋮               (     i   /   2     )     ⁢     β   M1     ⁢     e     -     i   ⁡     (       ϑ   M     -     ϑ   1       )                     (     i   /   2     )     ⁢     β     M   ⁢           ⁢   2       ⁢     e     -     i   ⁡     (       ϑ   M     -     ϑ   2       )                 …         k   M           ]         
 
is nonsingular. (Note that B mm =0 for every m.)
 
   We call M the reference matrix. Since the coupling coefficients B mn  are real-valued and skew-symmetric (that is, B nm =−B mn  for every m and n), this matrix is Hermitian. Note that M is a function of position (x, y). We shall assume that M is positive definite at every position (x, y). A sufficient condition for M to be positive definite is that this matrix be strictly diagonally dominant. 
                   k   n     &gt;       (     1   /   2     )     ⁢       ∑     m   ≠   n       ⁢            β     m   ⁢           ⁢   n            ⁢           ⁢   for   ⁢           ⁢   n           =   1     ,   …   ⁢           ,     M   .             (   2.11   )             
 
   There is a useful alternative formulation of the coupled parabolic equations. First, define 
                     ⁢         ψ   m     =         k   m     ⁢     φ   m       +       (     i   /   2     )     ⁢       ∑     n   =   1     M     ⁢       β     m   ⁢           ⁢   n       ⁢     φ   n     ⁢     e     -     i   ⁡     (       ϑ   m     -     ϑ   n       )                     ,             (   2.12   )                 q   m     =         (       k   m   2     -     ξ   m   2       )     ⁢     φ   m       +             (   2.13   )                       ⁢       ∑     n   =   1     M     ⁢     [         (       α     m   ⁢           ⁢   n       +         k   m       2   ⁢           ⁢     k   m         ⁢     β     m   ⁢           ⁢   n         +     i   ⁢           ⁢     k   m     ⁢     β     m   ⁢           ⁢   n         -       ∂     β     m   ⁢           ⁢   n           ∂   x         )     ⁢     φ   n       +                                           ⁢         γ     m   ⁢           ⁢   n       ⁢       ∂     φ   n         ∂   y         +           ⁢       i     2   ⁢     k   m         ⁢       ∂   2       ∂     y   2         ⁢     (       β     m   ⁢           ⁢   n       ⁢     φ   n       )         ]     ⁢       e       -   i     ⁢           ⁢     (       ϑ   m     -     ϑ   n       )         .                           
 
Now the following equation is equivalent to eq.(2.8). 
                 2   ⁢   i   ⁢           ⁢     k   m     ⁢       ∂     ψ   m         ∂   x         -     i   ⁢           ⁢     k   m     ⁢     ψ   m       +         ∂   2     ⁢     ψ   m         ∂     y   2           =       k   m     ⁢       q   m     .               (   2.14   )             
 
Although eqs.(2.8,2.14) are similar, they differ in one important way: no partial derivative with respect to x of any dependent variable appears on the RHS of eq.(2.14).
 
   Note that we can solve the system of equations (2.12)(m=1, . . . , M) for the complex amplitudes φ m  because the reference matrix M is nonsingular. Define the complex column vectors
 
φ=[φ 1 , φ 2 , . . . , φ M ] T ,  (2.15) 
 
ψ=[ψ 1 , ψ 2 , . . . , ψ M ] T ,  (2.16) 
 
where the superscript T indicates transpose. Then, ψ=Mφ, so φ=M −1 ψ.
 
   Solutions of the coupled parabolic equations (2.8) obey a conservation law involving two quantities that are closely related to power flow in the waveguide. These quantities are 
                 J   P     =         ∑     m   =   1     M     ⁢       k   m     ⁢            φ   m          2         +       (     i   /   2     )     ⁢       ∑     m   ,     n   =   1       M     ⁢       β     m   ⁢           ⁢   n       ⁢       φ   _     m     ⁢     φ   n     ⁢     e     -     i   ⁡     (       ϑ   m     -     ϑ   n       )                     ,           (   2.17   )                   K   P     =       i   2     ⁡     [         ∑     m   =   1     M     ⁢     (         φ   m     ⁢       ∂     φ   m         ∂   y         -         ∂     φ   m         ∂   y       ⁢     φ   m         )       +       ∑     m   ,     n   =   1       M     ⁢       γ     m   ⁢           ⁢   n       ⁢     φ   m     ⁢     φ   n     ⁢     e     -     i   ⁡     (       ϑ   m     -     ϑ   n       )                 ]         ,           (   2.18   )             
 
where {overscore (φ)} is the complex conjugate of φ. The conservation law states that 
                   ∂     J   P         ∂   x       +       ∂     K   P         ∂   y         =   0.           (   2.19   )             
 
   Note that J p =φ H Mφ, where φ H   32  |{overscore (φ)} 1 , {overscore (φ)} 2 , . . . , {overscore (φ)} M | is the Hermitian transpose of φ. If φ≠0, then J p &gt;0 because M is positive definite. 
   Now let us discuss power flow in the waveguide and show how the quantities J p  and K p  are related to it. First, consider a vertical strip of infinitesimal width Δy that extends from the bottom of the waveguide at z=−L to the surface of the waveguide at z=0 and that is perpendicular to the x direction. The time-averaged power through the strip is (2ω) −1 JΔy, where 
             J   =       1     2   ⁢   i       ⁢       ∫     -   L     0     ⁢       1   ρ     ⁢     (         p   _     ⁢       ∂   p       ∂   x         -     p   ⁢       ∂     p   _         ∂   x           )     ⁢       ⅆ   z     .                   (   2.20   )             
 
   Second, consider a vertical strip of infinitesimal width Δx that extends from the bottom of the waveguide to the surface of the waveguide and that is perpendicular to the y direction. The time-averaged power through this strip is (2w) −1 KΔx, where 
             K   =       1     2   ⁢   i       ⁢       ∫     -   L     0     ⁢       1   ρ     ⁢     (         p   _     ⁢       ∂   p       ∂   y         -     p   ⁢       ∂     p   _         ∂   y           )     ⁢       ⅆ   z     .                   (   2.21   )             
 
   Finally, consider a vertical cylinder that extends from the bottom of the waveguide to the surface of the waveguide and that has a simply-connected cross-section A. Using Green&#39;s theorem, we find that the time-averaged power out of this cylinder is 
                   1     2   ⁢           ⁢   ω       ⁢       ∮     ∂   A       ⁢     (       J   ⁢           ⁢     ⅆ   y       -     K   ⁢           ⁢     ⅆ   x         )         =       1     2   ⁢           ⁢   ω       ⁢     ∫       ∫   A     ⁢       (         ∂   J       ∂   x       +       ∂   K       ∂   y         )     ⁢     ⅆ   x     ⁢     ⅆ   y               ,           (   2.22   )             
 
where the boundary curve σA has the usual counter-clockwise orientation. Note that the partial derivatives σJ/σx and σK/σy have the physical dimensions of power flux, that is, power/area×time. For this reason, we call J and K power flux integrals.
 
   Suppose that we evaluate these power flux integrals for a solution of the acoustic equation (1.1) that satisfies eq.(1.2) at the boundaries and eq.(1.4) at the interface. We can use these equations and the divergence theorem to show that the integral on the left-hand side (LHS) of eq.(2.22) vanishes. Clearly, the integral on the RHS of this equation vanishes as well. Since the cross-section A is arbitrary, it follows that the integrand of this integral must vanish. Thus, 
                   ∂   J       ∂   x       +       ∂   K       ∂   y         =   0           (   2.23   )             
 
for exact solutions of eqs.(1.1,1.2,1.4).
 
   If, however, we evaluate the power flux integrals J and K for an approximate solution of these equations that is constructed from complex amplitude solutions φ m  of the coupled parabolic equations (2.8), then we find that 
               J   =         1     2   ⁢   i       ⁢       ∑     m   =   1     M     ⁢     (           φ   _     m     ⁢       ∂     φ   m         ∂   x         -       φ   m     ⁢       ∂       φ   _     m         ∂   x           )         +     J   P         ,           (   2.24   )               K   =       K   P     .             (   2.25   )             
 
In this case, eq.(2.19) implies that eq.(2.23) is not valid. Instead, we find that 
                   ∂   J       ∂   x       +       ∂   K       ∂   y         =       1     2   ⁢   i       ⁢       ∑     m   =   1     M     ⁢       (           φ   _     m     ⁢         ∂   2     ⁢     φ   m         ∂     x   2           -       φ   m     ⁢         ∂   2     ⁢       φ   _     m         ∂     x   2             )     .                 (   2.26   )             
 
This equation agrees with eq.(2.23) only to the extent that the paraxial approximations (2.6) are satisfied.
 
3. Adaptive selection of reference wavenumbers to conserve power.
 
   There may be a cross-range W such that we can neglect the acoustic pressure in sections of the waveguide where |y|&gt;W. For example, this can happen if the source produces a bean. In this case, we can assume that a pair of hard vertical walls in the planes through y=±W transversely bound the waveguide. This assumption makes the waveguide a duct, and it imposes the following Neumann boundary conditions on the complex amplitudes φ m . 
                   ∂     ∂   y       ⁢       φ   m     ⁡     (     x   ±   W     )         =       0   ⁢           ⁢   for   ⁢           ⁢   m     =   1       ,   …   ⁢           ,     M   .             (   3.1   )             
 
In this case we also can assume that the coupling coefficients γ mn  vanish at the walls.
 
γ mn (x,±W)=0 for m,n=1, . . . , M.  (3.2) 
 
Under these assumptions, eqs.(2.18,2.25) imply that the power flux integrals K=K p  vanish at the walls.
 
K(x,±W)=K p (x,±W)=0.  (3.3) 
 
Thus, if we integrate the conservation law (2.19) from y=−W to y=W, then we find that 
                 ∂     ∂   x       ⁢       ∫     -   W     W     ⁢       J   P     ⁢     ⅆ   y           =   0.           (   3.4   )             
 
This equation implies that the integral ∫ −W   W J p dy is independent of the range.
 
   If eq.(2.23) were valid, then the longitudinal power (2w) −1 ∫ −W   W J dy through the duct would be independent of the range as well. Although eq.(2.23) is not valid for the coupled parabolic equations, we can adaptively select the reference wavenumbers to make the longitudinal power through the duct be independent of the range. If we integrate eq. (2.24) from y=−W to y=W, then we get 
                 ∫     -   W     W     ⁢     J   ⁢     ⅆ   y         =         ∑     m   =   1     M     ⁢       1     2   ⁢   i       ⁢       ∫     -   W     W     ⁢       (           φ   _     m     ⁢       ∂     φ   m         ∂   x         -       φ   m     ⁢       ∂       φ   _     m         ∂   x           )     ⁢     ⅆ   y             +       ∫     -   W     W     ⁢       J   P     ⁢       ⅆ   y     .                   (   3.5   )             
 
Since the integral ∫ −W   W J p dy is independent of the range, this equation implies that the longitudinal power through the duct also will be independent of the range if the sum on the RHS is constant. The best way to accomplish this is to require that each summand vanish, which yields the integral constraints 
                   ∫     -   W     W     ⁢       (           φ   _     m     ⁢       ∂     φ   m         ∂   x         -       φ   m     ⁢       ∂       φ   _     m         ∂   x           )     ⁢     ⅆ   y         =       0   ⁢           ⁢   for   ⁢           ⁢   m     =   1       ,   …   ⁢           ,     M   .             (   3.6   )             
 
   These constraints are equivalent to a system of coupled quadratic equations in the reference wavenumbers k m . Let k=[k 1 , k 2 , . . . , k M ] T . For m=1, . . . , M define the following functions g m (k). 
                       g   m     ⁡     (   k   )       =       ⁢         k   m   2     ⁢       ∫     -   W     W     ⁢              φ   m          2     ⁢     ⅆ   y           +       ∑     n   =   1     M     ⁢       k   n     ⁢           ⁢   Re   ⁢       ∫     -   W     W     ⁢       φ   m     ⁢   i   ⁢           ⁢     β     m   ⁢           ⁢   n       ⁢     φ   n     ⁢     ⅇ     -     ⅈ   ⁡     (       ϑ   m     -     ϑ   n       )           ⁢     ⅆ   y             +                     ⁢         ∫     -   W     W     ⁢       (                ∂     φ   m         ∂   y            2     -       ξ   m   2     ⁢            φ   m          2         )     ⁢     ⅆ   y         +                     ⁢       ∑     n   =   1     M     ⁢     Re   ⁢           ⁢       ∫     -   W     W     ⁢         φ   m     ⁡     (         α     m   ⁢           ⁢   n       ⁢     φ   n       +       β     m   ⁢           ⁢   n       ⁢       ∂     φ   n         ∂   x         +       γ     m   ⁢           ⁢   n       ⁢       ∂     φ   n         ∂   y           )       ⁢     ⅇ     -     ⅈ   ⁡     (       ϑ   m     -     ϑ   n       )           ⁢       ⅆ   y     .                           (   3.7   )             
 
Then the constraints (3.6) are satisfied if and only if
 
g m (k)=0 for m=1, . . . , M.  (3.8) 
 
   It is convenient to write this equation in vector form. Let g(k)=[g 1 (k), . . . , g M (k)] T . Then eq.(3.8) is equivalent to the equation g(k)=0. In practice we can use Newton&#39;s method to solve this equation iteratively. The Jacobian matrix σg/σk that is associated with g(k) has the elements 
                 ∂     g   m         ∂     k   n         =     {           2   ⁢     k   m     ⁢       ∫     -   W     W     ⁢              φ   m          2     ⁢     ⅆ   y                     if   ⁢           ⁢   m     =   n     ,     (       since   ⁢           ⁢     β     m   ⁢           ⁢   m         =   0     )                 Re   ⁢       ∫     -   W     W     ⁢         φ   _     m     ⁢   i   ⁢           ⁢     β   mn     ⁢     φ   n     ⁢     ⅇ     -     ⅈ   ⁡     (       θ   m     -     θ   n       )           ⁢     ⅆ   y                   if   ⁢           ⁢   m     ≠     n   .                       (   3.9   )             
 
Given an estimate k 0  of the solution vector k, we can compute a correction Δk to it as the solution of the linear equation 
                   ∂   g       ∂   k       ⁢     |     k   0       ⁢       ·   Δ     ⁢           ⁢   k       =     -       g   ⁡     (     k   0     )       .               (   3.10   )             
 
Then we replace k 0  with k 0 +Δk. This completes one iteration of Newton&#39;s method. The Jacobian matrix σg/σk is real symmetric. If no amplitude φ m  vanishes identically and if the reference matrix M is positive definite at every position, then the Jacobian matrix σg/σk also is positive definite.
 
   The solutions of eq.(3.8) are optimal, but they are expensive to compute. There are suboptimal ways to estimate these reference wavenumbers cheaply that can be useful. For example, if we neglect the coupling coefficients in eq.(3.7), then we get the following approximate solution of eq.(3.8). 
                 k   m   2     =       ∫     -   W     W     ⁢       (         ξ   m   2     ⁢            φ   m          2       -              ∂     φ   m         ∂   y            2       )     ⁢       ⅆ   y     /       ∫     -   W     W     ⁢              φ   m          2     ⁢     ⅆ   y                 ⁢           ⁢     
     ⁢         for   ⁢           ⁢   m     =   1     ,   …   ⁢           ,     M   .               (   3.11   )             
 
4. Exact integrals of the coupled parabolic equations.
 
   We can use Fourier cosine expansions to integrate the alternative coupled parabolic equations (2.14) subject to Neumann boundary conditions at the boundaries y=±W, 
                   ∂     ∂   y       ⁢       φ   m     ⁡     (     x   ,     ±   W       )         =       0   ⁢           ⁢   for   ⁢           ⁢   m     =   1       ,   …   ⁢           ,     M   .             (   4.1   )             
 
It is convenient to introduce the following operator notation. Let g(y) be a continuously differentiable function over the interval −W≦y≦W such that dg(±W)/dy=0. The Fourier cosine coefficients ĝ of g are defined by the integrals 
                   g   ^     n     =         ε   ^     n     ⁢       ∫   0     2   ⁢   W       ⁢       g   ⁡     (     s   -   W     )       ⁢   cos   ⁢     {       (     n   ⁢           ⁢     π   /   2     ⁢   W     )     ⁢   s     }     ⁢       ⅆ   s     /   W             ⁢           ⁢     
     ⁢         for   ⁢           ⁢   n     =   0     ,   1   ,   2   ,   …   ⁢           ,             (   4.2   )             
 
where {circumflex over (ε)} 0 =½ and {circumflex over (ε)} n =1 for n=1, 2, . . . . Then the Fourier cosine expansion 
               g   ⁡     (     s   -   W     )       =       ∑     n   =   0     ∞     ⁢         g   ^     n     ⁢   cos   ⁢     {       (     n   ⁢           ⁢     π   /   2     ⁢   W     )     ⁢   s     }     ⁢           ⁢     (     0   ≤   s   ≤     2   ⁢   W       )                 (   4.3   )             
 
is uniformly convergent. Let ĝ denote the sequence {ĝ 0 , ĝ 1 , . . . , ĝ n  . . . }. The function g and the sequence ĝ uniquely determine each other. We define the operators ℑ and ℑ −1  by the relations
 
ℑg=ĝ, ℑ −1 ĝ=g.  (4.4) 
 
   Next, we introduce the spectral parameter η. This parameter takes on the discrete values
 
η=nπ/2W for n=0,1,2, . . .  (4.5) 
 
Given a function h(η) and a sequence ℑg, we define the new sequence h(η)ℑg as follows.
 
h(η)ℑg={h(0)ĝ 0 , h(π/2W)ĝ 1 , . . . , h(ηπ/2W)ĝ n , . . . }.  (4.6) 
 
   The integral form of eq.(2.14) is 
                 ψ   m     ⁡     (     x     j   +   1       )       =           (         k   m     ⁡     (     x     j   +   1       )           k   m     ⁡     (     x   j     )         )       1   2       ⁢       -   1       ⁢     (     exp   ⁢     {       -     (     i   ⁢           ⁢       η   2     /   2       )       ⁢       ∫     x   j       x     j   +   1         ⁢       ⅆ   u         k   m     ⁡     (   u   )             }     ⁢   ⁢           ⁢       ψ   m     ⁡     (     x   j     )         )       -       (     i   /   2     )     ⁢       ∫     x   j       x     j   +   1         ⁢         (         k   m     ⁡     (     x     j   +   1       )           k   m     ⁡     (   x   )         )       1   2       ⁢       -   1       ⁢             (     exp   ⁢     {       -     (     i   ⁢           ⁢       η   2     /   2       )       ⁢       ∫     x   j       x     j   +   1         ⁢       ⅆ   u         k   m     ⁡     (   u   )             }     ⁢   ⁢           ⁢       q   m     ⁡     (   x   )         )     ⁢       ⅆ   x     .                         (   4.7   )             
 
This equation, which holds for m=1, . . . , M, is the basis of our method for integrating the coupled parabolic equations numerically.
 
5. Discrete approximation of the exact integrals.
 
   We shall introduce three numerical approximations in eq.(4.7) that lead to a system of linear equations for the values of the complex amplitudes (m at discrete mesh points in the x, y plane. First, we use the trapezoidal rule to approximate the integral from x=x j  to x=x j+1  on the RHS of this equation. Letting Δx=x j+1 −x j  and neglecting terms of ο((Δx) 3 ), we get 
                   ψ   m     ⁡     (     x     j   +   1       )       +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢       q   m     ⁡     (     x     j   +   1       )           =         (         k   m     ⁡     (     x     j   +   1       )           k   m     ⁡     (     x   j     )         )       1   2       ×               -   1       ⁢       (     exp   ⁢     {       -     (     i   ⁢           ⁢       η   2     /   2       )       ⁢       ∫     x   j       x     j   +   1         ⁢       ⅆ   u         k   m     ⁡     (   u   )             }     ⁢     ⁡     [         ψ   m     ⁡     (     x   j     )       +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢       q   m     ⁡     (     x   j     )           ]         )     .                   (   5.1   )             
 
   We also can use the trapezoidal rule to approximate the argument of the exponential function in eq.(5.1). Neglecting terms of ο((Δx) 3 ), we get 
                 ∫     x   j       x     j   +   1         ⁢       ⅆ   u         k   m     ⁡     (   u   )           =       (         k   m     -   1       ⁡     (     x   j     )       +       k   m     -   1       ⁡     (     x     j   +   1       )         )     ⁢   Δ   ⁢           ⁢     x   /   2.               (   5.2   )             
 
We can approximate the phases σ m (x j+1 ) the same way. Neglecting terms of ο((Δx) 3 ), we get 
                       ϑ   m     ⁡     (     x     j   +   1       )       =       ⁢         ϑ   m     ⁡     (     x   j     )       +       ∫     x   j       x     j   +   1         ⁢         k   m     ⁡     (   u   )       ⁢     ⅆ   u                       =       ⁢         ϑ   m     ⁡     (     x   j     )       +       (         k   m     ⁡     (     x   j     )       +       k   m     ⁡     (     x     j   +   1       )         )     ⁢   Δ   ⁢           ⁢     x   /   2.                       (   5.3   )             
 
   The second numerical approximation that we introduce is to discretize the operators ℑ and ℑ −1  using discrete cosine transforms. Fix a positive integer N, and let g={g 0 , g 1 , . . . , g N ) be a sequence of N+1 complex numbers. Define the discrete Fourier cosine coefficients {overscore (g)} n  as follows. 
                   g   _     n     =             ε   _     n     ⁡     (     2   /   N     )       ⁢       ∑     k   =   0     N     ⁢         ε   _     k     ⁢     g   k     ⁢     cos   ⁡     (     π   ⁢           ⁢     kn   /   N       )       ⁢           ⁢   for   ⁢           ⁢   n         =   0       ,   1   ,   …   ⁢           ,   N   ,     
     ⁢       where   ⁢           ⁢       ε   _     0       =         ε   _     N     =         1   2     ⁢           ⁢   and   ⁢           ⁢       ε   _     n       =       1   ⁢           ⁢   for   ⁢           ⁢   n     =   1           ,   …   ⁢           ,     N   -     1.   ⁢           ⁢   Now               (   5.4   )             
 for n=0, 1, . . . , N,  (5.4)
 
where {overscore (ε)} 0 ={overscore (ε)} N =½ and {overscore (ε)} n =1 for n=1, . . . , N−1. Now 
                 g   k     =         ∑     n   =   0     N     ⁢         g   _     n     ⁢     cos   ⁡     (     π   ⁢           ⁢     kn   /   N       )       ⁢           ⁢   for   ⁢           ⁢   k       =   0       ,   1   ,   …   ⁢           ,   N           (   5.5   )             
 for k=0, 1, . . . , N  (5.5)
 
if and only if the coefficients {overscore (g)} n  are given by eq.(5.4). Let {overscore (g)}={{overscore (g)} 0 , {overscore (g)} 1  , . . . , {overscore (g)} N } denote the sequence of these coefficients. We define the operators ℑ N  and ℑ N   −1  by the relations
 
ℑ N g={overscore (g)}, ℑ N   −1 {overscore (g)}=g  (5.6) 
 
   We replace the operators ℑ and ℑ −1  in eq.(5.1) by the operators ℑ N  and ℑ n   −1 . Since these operators transform sequences to sequences, however, we must also replace the functions ψ m ±(iΔx/4)q m  in this equation with the sequences of their values at the points y k , which we write as {ψ m (y k )±(iΔx/4)q m (y k )} k=0   N . Under these approximations, eq.(5.1) becomes 
                 {         ψ   m     ⁡     (       x     j   +   1       ,     y   k       )       +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢       q   m     ⁡     (       x     j   +   1       ,     y   k       )           }       k   =   0     N     =         (         k   m     ⁡     (     x     j   +   1       )           k   m     ⁡     (     x   j     )         )       1   2       ×             N     -   1       ⁢       (     exp   ⁢     {       -     (     i   ⁢           ⁢       η   2     /   2       )       ⁢       ∫     x   j       x     j   +   1         ⁢       ⅆ   u         k   m     ⁡     (   u   )             }     ⁢     N     ⁢       {         ψ   m     ⁡     (       x   j     ,     y   l       )       -       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢       q   m     ⁡     (       x   j     ,     y   l       )           }       l   =   0     N       )     .                   (   5.7   )             
 
Note that the spectral parameter η in this equation is bounded because η=ηπ/2W for n=0, 1, . . . , N. Hence, 0≦η≦π/Δy.
 
   We must approximate the partial derivatives σφ n /σy and σ 2 (β mn φ n )/σy 2  that appear in eq.(2.13) for q m  before we can evaluate the sequences in eq.(5.7). Since the variable y in this equation is restricted to the uniformly spaced points y k , our third approximation is to replace the partial derivatives at these points with central differences. At the lowest order of approximation, we use the three-point formulas 
                     ∂               ∂   y       ⁢       φ   n     ⁡     (     y   k     )         =         1     2   ⁢   Δ   ⁢           ⁢   y       ⁡     [         φ   n     ⁡     (     y     k   +   1       )       -       φ   n     ⁡     (     y     k   -   1       )         ]       +     O   ⁡     (       (     Δ   ⁢           ⁢   y     )     2     )           ,           (   5.8   )                       ∂   2     ⁢               ∂     y   2         ⁢     (       β   mn     ⁢     φ   n       )     ⁢     (     y   k     )       =         1       (     Δ   ⁢           ⁢   y     )     2       ⁡     [         (       β   mn     ⁢     φ   n       )     ⁢     (     y     k   +   1       )       -     2   ⁢     (       β   mn     ⁢     φ   n       )     ⁢     (     y   k     )       +       (       β   mn     ⁢     φ   n       )     ⁢     (     y     k   -   1       )         ]       +       O   ⁡     (       (     Δ   ⁢           ⁢   y     )     2     )       .               (   5.9   )             
 
Now we combine eqs.(5.8,5.9) with eq.(2.13). Rearranging terms in the result, we find that 
                 q   m     ⁡     (     y   k     )       =         (       k   m   2     -       ξ   m   2     ⁡     (     y   k     )         )     ⁢       φ   m     ⁡     (     y   k     )         +       ∑     n   =   1     M     ⁢       [         α   mn     ⁡     (     y   k     )       +       (           k   .     m       2   ⁢     k   m         +     ik   m     -     i         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         )     ⁢       β   mn     ⁡     (     y   k     )         -         ∂               ∂   x       ⁢       β   mn     ⁡     (     y   k     )           ]     ⁢     e     -     i   ⁡     (       θ   m     -     θ   n       )           ⁢       φ   n     ⁡     (     y   k     )           +       ∑     n   =   1     M     ⁢       (         1     2   ⁢   Δ   ⁢           ⁢   y       ⁢       γ   mn     ⁡     (     y   k     )         +       i     2   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (     y     k   +   1       )           )     ⁢     e     -     i   ⁡     (       θ   m     -     θ   n       )           ⁢       φ   n     ⁡     (     y     k   +   1       )           +       ∑     n   =   1     M     ⁢       (         -     1     2   ⁢   Δ   ⁢           ⁢   y         ⁢       γ   mn     ⁡     (     y   k     )         +       i     2   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (     y     k   -   1       )           )     ⁢     e     -     i   ⁡     (       θ   m     -     θ   n       )           ⁢       φ   n     ⁡     (     y     k   -   1       )           +       O   ⁡     (       (     Δ   ⁢           ⁢   y     )     2     )       .               (   5.10   )             
 
   Therefore, neglecting terms of ο(Δx(Δy) 2 ), we find that 
                   ψ   m     ⁡     (       x     j   +   1       ,     y   k       )       +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢       q   m     ⁡     (       x     j   +   1       ,     y   k       )           =                 [       k   m     +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢     (       k   m   2     -       ξ   m   2     ⁡     (       x     j   +   1       ,     y   k       )         )         ]     ⁢       φ   m     ⁡     (       x     j   +   1       ,     y   k       )         +       (       i   2     +       i   ⁢       k   .     m     ⁢   Δ   ⁢           ⁢   x       8   ⁢     k   m         -         k   m     ⁢   Δ   ⁢           ⁢   x     4     +       Δ   ⁢           ⁢   x       4   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2           )     ⁢       ∑     n   =   1     M     ⁢         β   mn     ⁡     (       x     j   +   1       ,     y   k       )       ⁢     e     -     i   ⁡     (       θ   m     -     θ   n       )           ⁢       φ   n     ⁡     (       x     j   +   1       ,     y   k       )             +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢       ∑     n   =   1     M     ⁢       (         α   mn     ⁡     (       x     j   +   1       ,     y   k       )       -         ∂               ∂   x       ⁢       β   mn     ⁡     (       x     j   +   1       ,     y   k       )           )     ⁢     e     -     i   ⁡     (       θ   m     -     θ   n       )           ⁢       φ   n     ⁡     (       x     j   +   1       ,     y   k       )             +       ∑     n   =   1     M     ⁢       (           i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢   Δ   ⁢           ⁢   y       ⁢       γ   mn     ⁡     (       x     j   +   1       ,     y   k       )         -         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (       x     j   +   1       ,     y     k   +   1         )           )     ⁢     e     -     i   ⁡     (       θ   m     -     θ   n       )           ⁢       φ   n     ⁡     (       x     j   +   1       ,     y     k   +   1         )           +       ∑     n   =   1     M     ⁢       (         -       i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢   Δ   ⁢           ⁢   y         ⁢       γ   mn     ⁡     (       x     j   +   1       ,     y   k       )         -         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (       x     j   +   1       ,     y     k   -   1         )           )     ⁢     e     -     i   ⁡     (       θ   m     -     θ   n       )           ⁢       φ   n     ⁡     (       x     j   +   1       ,     y     k   -   1         )             ,               (   5.11   )             
 
where k n =k n (x j+1 ) and σ n =σ n (x j+1 ) for n=1, . . . , M.
 
   Similarly, neglecting terms of ο(Δx(Δy) 2 ), we find that 
                     ψ   m     ⁡     (       x   j     ,     y   l       )       -       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢       q   m     ⁡     (       x   j     ,     y   l       )           =         [           ⁢       k   m     -       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢     (       k   m   2     -       ξ   m   2     ⁡     (       x   j     ,     y   l       )         )         ]     ⁢           ⁢       φ   m     ⁡     (       x   j     ,     y   l       )         +       (       i   2     -         ik   m     ⁢   Δ   ⁢           ⁢   x       8   ⁢     k   m         +         k   m     ⁢   Δ   ⁢           ⁢   x     4     -       Δ   ⁢           ⁢   x       4   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2           )     ⁢           ⁢       ∑     n   =   1     M     ⁢         β   mn     ⁡     (       x   j     ,     y   l       )       ⁢           ⁢     e     -     i   ⁡     (       θ   m     -     θ   n       )           ⁢           ⁢       φ   n     ⁡     (       x   j     ,     y   l       )             -       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢       ∑     n   =   1     M     ⁢       (         α   mn     ⁢     (       x   j     ,     y   l       )       -         ∂               ∂   x       ⁢       β   mn     ⁡     (       x   j     ,     y   l       )           )     ⁢     e     -     i   ⁡     (       θ   m     -     θ   n       )           ⁢       φ   n     ⁡     (       x   j     ,     y   l       )             +       ∑     n   =   1     M     ⁢       (         -       i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢   Δ   ⁢           ⁢   y         ⁢     γ   mn     ⁢     (       x   j     ,     y   l       )       +         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (       x   j     ,     y     l   +   1         )           )     ⁢           ⁢     e     -     i   ⁡     (       θ   m     -     θ   n       )           ⁢           ⁢       φ   n     (           ⁢       x   j     ,     y     l   +   1         )         +       ∑     n   =   1     M     ⁢       (         -       i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢   Δ   ⁢           ⁢   y         ⁢     γ   mn     ⁢     (       x   j     ,     y   l       )       +         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢           ⁢       β   mn     (           ⁢       x   j     ,           ⁢     y     l   -   1         )         )     ⁢           ⁢     e     -     i   ⁡     (       θ   m     -     θ   n       )           ⁢           ⁢       φ   n     (           ⁢       x   j     ,     y     l   -   1         )             ,           (   5.12   )             
 
where k n =k n (x j ) and σ n =σ n (x j ) for n=1, . . . , M.
 
   We need additional information to interpret these equations at the boundaries y=±W. Let us require that eq.(4.1) be consistent with eq (3.1). It follows that for k=0, N and for all x j ,
 
φ n (x j ,y k−1 )+φ n (x j ,y k+1 ) for n=1, . . . ,M,  (5.13) 
 
β mn (x j ,y k−1 )=β mn (x j ,y k+1 ) for m,n=1, . . . ,M.  (5.14) 
 
   Equations (5.7,5.11), subject to eqs (5.13,5.14), determine a system of linear equations for the values φ m (x j+1 , y k ) at a range step x j+1 , which in turn depend linearly on the values φ m (x j , y 1 ) at the previous range step x j  through eq.(5.12). In order to solve these equations numerically, it is helpful to put them in matrix form. First, we write eq.(5.7) as follows. 
                   {       Ψ   m     ⁡     (         x     j   +   1       -   0     ,     y   k       )       }       k   =   0     N     =         (         k   m     ⁡     (     x     j   +   1       )           k   m     ⁡     (     x   j     )         )       1   2       ×     N     -   1       ⁢           ⁢     (     exp   ⁢     {       -     (     i   ⁢           ⁢       η   2     /   2       )       ⁢       ∫     x   i       x     i   +   1         ⁢       ⅆ   u         k   m     ⁡     (   u   )             }     ⁢     N     ⁢       {       Ψ   m     ⁡     (         x   j     +   0     ,     y   l       )       }       l   =   0     N       )         ,           (   5.15   )             
 
where
 
Ψ m (x j+1 −0,y k )=ψ m (x j+1 ,y k )+(iΔx/4)q m (x j+1 ,y k ),  (5.16) 
 
Ψ m (x j +0,y 1 )=ψ m (x j ,y 1 )−(iΔx/4)q m (x j ,y 1 ).  (5.17) 
 
Second, for k=0,1, . . . , N we define the column vectors
 
φ k =[φ 1 (x j+1 ,y k ), . . . , φ M (x j+1 ,y k )] T ,  (5.18) 
 
Ψ k   − =[Ψ 1 (x j+1 −ο,y k ), . . . , Ψ M (x j+1 −ο, y k )] T .  (5.19) 
 
Third, for k=0,1, . . . , N we define the M×M matrices A kk , A kk+1  and A kk−1  as follows. If k n =k n (x j+1 ) and σ n =σ n (x j+1 ) for n=1, . . . , M, then 
                 (     A   kk     )     mn     =     {               k   n     +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢     (       k   n   2     -       ξ   n   2     ⁡     (       x     j   +   1       ,     y   k       )       +       α   nn     ⁡     (       x     j   +   1       ,     y   k       )         )         ,           (       if   ⁢           ⁢   m     =   n     )               [       (       i   2     +       i   ⁢       k   .     m     ⁢   Δ   ⁢           ⁢   x       8   ⁢     k   m         -         k   m     ⁢   Δ   ⁢           ⁢   x     4     +       Δ   ⁢           ⁢   x       4   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2           )     ⁢       β   mn     ⁡     (       x     j   +   1       ,     y   k       )                                     +     (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )       ⁢     (         α   mn     ⁡     (       x     j   +   1       ,     y   k       )       -         ∂               ∂   x       ⁢       β   mn     ⁡     (       x     j   +   1       ,     y   k       )           )     ⁢     e     -     i   ⁡     (       θ   m     -     θ   n       )             ,           (       if   ⁢           ⁢   m     ≠   n     )                     (   5.20   )                   (     A     kk   ±   1       )     mn     =       (         ±       i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢           ⁢   Δ   ⁢           ⁢   y         ⁢       γ   mn     ⁡     (       x     j   +   1       ,     y   k       )         -         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (       x     j   +   1       ,     y     k   ±   1         )           )     ⁢       e     -     i   ⁡     (       θ   m     -     θ   n       )           .               (   5.21   )             
 
Finally, we define the M×M matrices B 1  and B N−1  as follows. For k=1, N−1, 
                 (     B   k     )     mn     =       -       Δ   ⁢           ⁢   x       4   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2           ⁢       β   mn     ⁡     (       x     j   +   1       ,     y   k       )       ⁢       e     -     i   ⁡     (       θ   m     -     θ   n       )           .               (   5.22   )             
 
Note that B 1 =A ο−1 +A ο1  and B N−1 =A NN−1 +A NN+1 . Now eqs.(5.11,5.16), subject to eqs.(5.13,5.14), are equivalent to the following matrix equation. 
                 [           A   00           B   1                                                   A   10           A   11           A   12                                                 ⋰       ⋰       ⋰                                                 A     N   -     1   ⁢   N     -   2             A     N   -     1   ⁢   N     -   1             A     N   -     1   ⁢   N                                                       B     N   -   1             A   NN           ]     ·     [           φ   0               φ   1             ⋮             φ     N   -   1                 φ   N           ]       =       [           Ψ   0   -               Ψ   1   -             ⋮             Ψ     N   -   1     -               Ψ   N   -           ]     .             (   5.23   )             
 
   We recommend solving this equation numerically by Gaussian elimination with partial pivoting (GEPP) followed by one round of iterative refinement. The block tri-diagonal system matrix on the LHS of eq.(5.23) is also a band matrix of order M(N+1) that has 2M−1 lower diagonals and 2M−1 upper diagonals. Since M&lt;&lt;N in practice, it is most efficient to do the GEPP computations with subroutines that take advantage of band, storage. 
   Let us write eq.(5.23) as A·Φ=Ψ − , and let Φ 0  be the numerical solution of this equation computed in finite precision arithmetic by GEPP. Iterative refinement is a technique to improve the accuracy of this solution. One round of this technique has three steps. 
   (1) Compute the residual τ 0 =Ψ − −A·Φ 0 . 
   (2) Using the LU decomposition of A from the original solution of eq. (5.23) by GEPP, solve the equation A·δ 0 =τ 0  for the correction δ 0 . 
   (3) Refine the solution: Φ t =Φ 0 +δ 0 . 
   The computations in these steps are done at the same arithmetic precision as the original solution. 
   The system matrix on the LHS of eq.(5.23) is nonsingular if it is strictly diagonally dominant by columns. We see from eqs.(5.20,5.21) that this condition requires the following strict inequalities to be satisfied for the block indices k=2, . . . , N−2: 
                          k   n     +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢     (       k   n   2     -       ξ   n   2     ⁡     (       x     j   +   1       ,     y   k       )       +       a   nn     ⁡     (       x     j   +   1       ,     y   k       )         )              &gt;         ∑       m   =   1       (     m   ≠   n     )       M     ⁢              (       i   2     +         ik   m     ⁢   Δ   ⁢           ⁢   x       8   ⁢     k   m         -         k   m     ⁢   Δ   ⁢           ⁢   x     4     +       Δ   ⁢           ⁢   x       4   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2           )     ⁢       β     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )         +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢     (         a     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )       -         ∂               ∂   x       ⁢       β     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )           )                +       ∑     m   =   1     M     ⁢                i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢           ⁢   Δ   ⁢           ⁢   y       ⁢       γ     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y     k   -   1         )         -         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )                  +       ∑     m   =   1     M     ⁢                  i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢           ⁢   Δ   ⁢           ⁢   y       ⁢       γ     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y     k   +   1         )         -         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )                ⁢           ⁢   for   ⁢           ⁢   n           =   1     ,   …   ⁢           ,     M   .             (   5.24   )             
 
These inequalities are modified slightly for the remaining indices k=0, 1, N−1, N. This condition is more restrictive than the similar condition (2.11) for the reference matrix M to be positive definite at the points y k . If this condition is satisfied, then Gaussian elimination requires no pivoting and is stable.
 
   Note the factors k m Δx, Δx/Δy and Δx/k m (Δy) 2  that appear in this equation. We recommend choosing the increments Δx and Δy so that k m Δx≦ο(1) and k m Δy≧ο(1). These restrictions imply that the ratios Δx/Δy and Δr/k m (Δy) 2  are at most ο(1). These restrictions also limit the errors that we make by applying the trapezoidal rule to eq.(4.7) to obtain eq.(5.1). 
   Finally, we note that higher order central difference formulas can be used in place of eqs.(5.8,5.9). For example, at the next order of approximation, we would use five-point formulas. They lead to a matrix formulation like eq.(5.23) in which the system matrix is block penta-diagonal. 
   6. Discrete transparent boundary conditions. 
   The most simple standard parabolic equation is 
                   2   ⁢     ik   0     ⁢       ∂   φ       ∂   x         +         ∂   2     ⁢   φ       ∂     y   2           =   0     ,           (   6.1   )             
 
where the wavenumber k 0  is a positive constant. Let us assume
 
   (1) the solution φ(x, y) is defined in a region of the x, y plane that includes the upper right quadrant where x≧0 and y≧0, 
   (2) the solution φ(x, y) and its partial derivatives are bounded in this quadrant, 
   (3) φ(0, y)=0 for y≧0.
 
We can translate these assumptions into a transparent boundary condition that the solution φ(z, y) satisfies on the right-half x-axis (x&gt;0, y=0). This condition makes the half-axis transparent in the sense that the wavefield described by the solution φ(x, y) propagates across this half-axis into the upper right quadrant without reflection. 
               φ   ⁡     (     x   ,   0     )       =       -         e     i   ⁢           ⁢     π   /   4         ⁡     (     2   ⁢   π   ⁢           ⁢     k   0       )           -   1     /   2         ⁢       ∫   0   x     ⁢         ∂               ∂   y       ⁢     φ   ⁡     (     u   ,   0     )       ⁢       (     x   -   u     )         -   1     /   2       ⁢       ⅆ   u     .                   (   6.2   )             
 
If we change the sign of the partial derivative σ/σy in this equation, then we get a new condition that makes the right-half x-axis transparent in the opposite sense. Now the wavefield described by the solution φ(x, y) propagates across this half-axis into the lower right quadrant without reflection. 
               φ   ⁡     (     x   ,   0     )       =           e     i   ⁢           ⁢     π   /   4         ⁡     (     2   ⁢   π   ⁢           ⁢     k   0       )           -   1     /   2       ⁢       ∫   0   x     ⁢         ∂               ∂   y       ⁢     φ   ⁡     (     u   ,   0     )       ⁢       (     x   -   u     )         -   1     /   2       ⁢       ⅆ   u     .                   (   6.3   )             
 
Let us replace the reference wavenumber k 0  in eq.(6.1) with a function k(x)&gt;0, 
                 2   ⁢     ik   ⁡     (   x   )       ⁢       ∂   φ       ∂   x         +         ∂   2     ⁢   φ       ∂     y   2           =   0.           (   6.4   )             
 
The change of variable 
               τ   ⁡     (   x   )       =       ∫   0   x     ⁢       ⅆ   u       k   ⁡     (   u   )                   (   6.5   )             
 
transforms this equation into a standard parabolic equation, 
                 2   ⁢   i   ⁢       ∂   φ       ∂   τ         +         ∂   2     ⁢   φ       ∂     y   2           =   0.           (   6.6   )             
 
Let {x 0 ,x 1 , . . . ,xj, . . . } be an increasing sequence of discrete range steps, with x 0 =0. Define ψ j (y)=φ(x, y) and τ j =τ(x j ) for j=0,1, . . . . Using the trapezoidal rule to integrate eq.(6.6) from τ j  to τ j+1  and neglecting the ο((τ j+1 −τ j ) 3 ) remainder, we find that 
                 2   ⁢     i   ⁡     [         ψ     j   +   1       ⁡     (   y   )       -       ψ   j     ⁡     (   y   )         ]         +       1   2     ⁢       (       τ     j   +   1       -     τ   j       )     ⁡     [           ∂   2       ∂     y   2         ⁢       ψ   j     ⁡     (   y   )         +         ∂   2       ∂     y   2         ⁢       ψ     j   +   1       ⁡     (   y   )           ]           =   0.           (   6.7   )             
 
Similarly, 
                 τ     j   +   1       -     τ   j       =         ∫     x   j       x     j   +   1         ⁢       d   ⁢           ⁢   u       k   ⁡     (   u   )           =       1   2     ⁢         (       x     j   +   1       -     x   j       )     ⁡     [         k     -   1       ⁡     (     x   j     )       +       k     -   1       ⁡     (     x     j   +   1       )         ]       .                 (   6.8   )             
 
where we evaluated the integral by the trapezoidal rule and dropped the ο((x j+1 −x j ) 3 ) remainder. We get a semi-discrete parabolic equation from eqs.(6.7,6.8). 
                     ∂   2       ∂     y   2         ⁢       ψ   j     ⁡     (   y   )         +         ∂   2       ∂     y   2         ⁢       ψ     j   +   1       ⁡     (   y   )         +     i   ⁢       8       (       x     j   +   1       -     x   j       )     ⁡     [         k     -   1       ⁡     (     x   j     )       +       k     -   1       ⁡     (     x     j   +   1       )         ]         ⁡     [         ψ     j   +   1       ⁡     (   y   )       -       ψ   j     ⁡     (   y   )         ]           =   0.           (   6.9   )             
 
   Let us assume 
   (1) the functions ψ j (y) and their derivatives are uniformly bounded for y≧0, 
   (2) ψ 0 (y)=0 for y≧0.
 
We can translate these assumptions into discrete transparent boundary conditions that the functions ψ j+1 (y)(j=0,1, . . . ) satisfy at y=0. A newly created algorithm computes these conditions recursively. First, define 
                 p     j   +   1       =             e       -   i     ⁢           ⁢     π   /   4         ⁡     (     8       (       x     j   +   1       -     x   j       )     ⁡     [         k     -   1       ⁡     (     x   j     )       +       k     -   1       ⁡     (     x     j   +   1       )         ]         )         1   2       ⁢           ⁢   for   ⁢           ⁢   j     =   0       ,   1   ,   …           (   6.10   )             
 
Second, define 
                 g   k   j     =           p   k     -     p   j           p   k     +     p   j         ⁢           ⁢   for   ⁢           ⁢   j       ,     k   =   1     ,   2   ,   …           (   6.11   )             
 
Note that g k   j  is real and |g k   j  |&lt;1 for all j, k.
 
   If j=1, then we let 
               a   1   1     =         1   2     ⁢     p   1     -   2       ⁢     ∂     ∂   y       ⁢       ψ   1     ⁡     (   0   )       ⁢           ⁢   and   ⁢           ⁢     a   2   1       =       -     1   2       ⁢     p   1     -   2       ⁢     ∂     ∂   y       ⁢         ψ   1     ⁡     (   0   )       .                 (   6.12   )             
 
If j&gt;1, then we compute the coefficients a 1   j , . . . , a j+1   j  from the coefficients a 1   j−1 , . . . , a j   j−1  and from the derivatives σψ j−1 (ο)/σy and σψ j (ο)/σy in two steps.
 
   (1) Compute ā 1   j−1 , . . . , ā j   j−1  and {overscore (b)} 1   j−1 , . . . , {overscore (b)} j   j−1  from the following recurrence formulas.
 
ā 1   j−1 =a 1   j−1 ,  (6.13) 
 
ā k   j−1 =ā k   j−1 −g k−1   j ā k−1   j−1  for k=2, . . . ,j,  (6.14) 
 
{overscore (b)} 1   j−1 =g 1   j ā 1   j−1 ,  (6.15) 
 
{overscore (b)} k   j−1 =g k   j ā k   j−1 +ā k−1   j−1 −g k−1   j {overscore (b)} k−1   j−1  for k=2, . . . ,j.  (6.16) 
 
   (2). Compute the coefficients a 1   j , . . . , a j+1   j  as follows. 
                 a   1   j     =         -     1   2       ⁢     a   1     j   -   1         -       1   2     ⁢     g   1   j     ⁢       b   _     1     j   -   1             ,           (   6.17   )                   a   k   j     =           -     1   2       ⁢     a   k     j   -   1         -       1   2     ⁢     g   k   j     ⁢       b   _     k     j   -   1         -       1   2     ⁢       b   _       k   -   1       j   -   1       ⁢           ⁢   for   ⁢           ⁢   k       =   2       ,   …   ⁢           ,     j   -   1     ,           (   6.18   )             
 
   (Note: omit the preceding equation if j=2.) 
                 a   j   j     =         1   2     ⁢       p   j     -   2       ⁡     [           ∂               ∂   y       ⁢       ψ     j   -   1       ⁡     (   0   )         +         ∂               ∂   y       ⁢       ψ   j     ⁡     (   0   )           ]         -       1   2     ⁢     a   j     j   -   1         -       1   2     ⁢     a   j     j   -   1         -       1   2     ⁢       b   _       j   -   1       j   -   1             ,           (   6.19   )                 a     j   +   1     j     =         1   2     ⁢       p   j     -   2       ⁡     [           ∂               ∂   y       ⁢       ψ     j   -   1       ⁡     (   0   )         +         ∂               ∂   y       ⁢       ψ   j     ⁡     (   0   )           ]         -       1   2     ⁢         b   _     j     j   -   1       .                 (   6.20   )             
 
The coefficients ā j+1   j  are used to compute the discrete transparent boundary conditions as follows. 
                       ∂               ∂   y       ⁢       ψ   j     ⁡     (   0   )         +       p     j   +   1       ⁢       ψ   j     ⁡     (   0   )           =   0     ,           (   6.21   )                         ∂               ∂   y       ⁢       ψ     j   +   1       ⁡     (   0   )         +       p     j   +   1       ⁢       ψ     j   +   1       ⁡     (   0   )           =           -       ∂               ∂   y         ⁢       ψ   j     ⁡     (   0   )         -       p     j   +   1       ⁢       ψ   j     ⁡     (   0   )         +     2   ⁢     p     j   +   1     2     ⁢       a   _       j   +   1     j     ⁢           ⁢   for   ⁢           ⁢   j       =   1       ,   2   ,   …           (   6.22   )             
 
   This algorithm can be greatly simplified if the quantities p j  are independent of the index j. This happens when the range increments x j+1 −x j  and the reference wavenumber values k(x j ) are independent of the index j. Let us assume that x j =jΔx and k(x j )=k 0  for j=0,1, . . . , where the increment Δx and the wavenumber k 0  are fixed. Then it follows from eq. (6.10) that p j =p 0  for every j, where
 
 p   0   =e   ix/4 2( k   0   /Δx ) 1/2 .  (6.23) 
 
In addition, eq.(6.11) implies that g k   j =0 for all j,k. Thus, eqs.(6.13,6.14) become
 
ā k   j−1 =a k   j−1  for k=1, . . . ,j,  (6.24) 
 
eq.(6.15) becomes {overscore (b)} 1   j−1 =0 and eq.(6.16) becomes
 
{overscore (b)} k   j−1 =ā k−1   j−1  for k=2, . . . ,j.  (6.25) 
 
   We can use these equations to eliminate the coefficients ā k   j−1  and {overscore (b)} k   j−1  from eqs.(6.17-6.20) if we extend the definition of the coefficients a k   j  as follows. Let
 
a −1   j =0 for j=0,1, . . . ,  (6.26) 
 
a 0   j =0 for j=0,1 . . .  (6.27) 
 
Using these extensions, we get the following simplified forms of eqs.(6.17-6.20). For j&gt;1, 
                 a   k   j     =         -     1   2       ⁢     a   k     j   -   1         -       1   2     ⁢     a   k     j   -   1         -       1   2     ⁢     a     k   -   2       j   -   1             ,       for   ⁢           ⁢   k     =   1     ,   …   ⁢           ,     j   -   1     ,           (   6.28   )                   a   j   j     =         1   2     ⁢       p   0     -   2       ⁡     [         ∂     ∂   y       ⁢       ψ     j   -   1       ⁡     (   0   )         +       ∂     ∂   y       ⁢       ψ   j     ⁡     (   0   )           ]         -     a   j     j   -   1       -       1   2     ⁢     a     j   -   2       j   -   1             ,           (   6.29   )                 a     j   +   1     j     =         -     1   2       ⁢       p   0     -   2       ⁡     [         ∂     ∂   y       ⁢       ψ     j   -   1       ⁡     (   0   )         +       ∂     ∂   y       ⁢       ψ   j     ⁡     (   0   )           ]         -       1   2     ⁢       a     j   -   1       j   -   1       .                 (   6.30   )             
 
   We also define a 1   0 =0. This extension has the following useful consequences. 
   (1) We can let j=1 in eqs.(6.29,6.30), which become eq.(6.12) because σψ 0 (0)/σy=0. 
   (2) Because of eq.(6.24), we can get eq.(6.21) from eq.(6.22) by letting j=0. Therefore, we can write these discrete transparent boundary conditions as follows. 
                     ∂     ∂   y       ⁢       ψ     j   +   1       ⁡     (   0   )         +       p   0     ⁢       ψ     j   +   1       ⁡     (   0   )           =         -     ∂     ∂   y         ⁢       ψ   j     ⁡     (   0   )         -       p   0     ⁢       ψ   j     ⁡     (   0   )         +     2   ⁢     p   0   2     ⁢     a     j   +   1     j           ⁢     
     ⁢         for   ⁢           ⁢   j     =   0     ,   1   ,   …             (   6.31   )             
 
   (3) An alternative form of these conditions is obtained by iterating eq.(6.31) backwards. The initial conditions ψ 0 (0)=0 and σψ 0 (0)/σy=0 are used to stop the iteration. 
                     ∂     ∂   y       ⁢       ψ     j   +   1       ⁡     (   0   )         +       p   0     ⁢       ψ     j   +   1       ⁡     (   0   )           =     2   ⁢     p   0   2     ⁢       ∑     k   =   0     j     ⁢         (     -   1     )       j   -   k       ⁢     a     k   +   1     k             ⁢           ⁢     
     ⁢         for   ⁢           ⁢   j     =   0     ,   1   ,   …             (   6.32   )             
 
   When eqs. (6.286.30) are evaluated in finite-precision computer arithmetic, roundoff errors in the additions tend to accumulate. However, we can use compensated summation to identify these errors and to cancel them to working precision. Let us describe how this is done in a simple case. 
   Consider the following sequence of differential boundary values σψ j (0)/σy. 
                 ∂     ∂   y       ⁢       ψ   j     ⁡     (   0   )         =     {           2   ⁢     p   0   2                 if   ⁢           ⁢   j     =   1     ,             0           if   ⁢           ⁢   j     ≠   1.                     (   6.33   )             
 
Since a 1   0 =0, it follows from eq. (6.32) that ψ 1 (0)=−2p 0  and 
                   ψ     j   +   1       ⁡     (   0   )       =       2   ⁢     p   0     ⁢       ∑     k   =   1     j     ⁢         (     -   1     )       j   -   k       ⁢     a     k   +   1     k     ⁢           ⁢   for   ⁢           ⁢   j         =   1       ,   2   ,   …           (   6.34   )             
 
Thus, ψ j+1 (0)=−2p 0 σ j  where 
                 σ   j     =       -       ∑     k   =   1     j     ⁢         (     -   1     )       j   -   k       ⁢     a     k   +   1     k     ⁢           ⁢   for   ⁢           ⁢   j         =   1       ,   2   ,   …   ⁢           ,           (   6.36   )             
 
which can be calculated recursively as follows.
 
σ j =−a j+1   j −σ j−1  for j=2,3, . . .  (6.36) 
 
   We can compute the coefficients a k   j  for j=1, 2, 3 exactly from eqs.(6.28-6.30). Using the input values (6.33), we find that 
                             a   1   1     =   1               a   2   1     =     -   1       ,                             a   1   2     =     -     1   2                 a   2   2     =   2               a   3   2     =     -     3   2         ,                             a   1   3     =     1   4               a   2   3     =     -   1               a   3   3     =     7   4               a   4   3     =     -   1.                   }           (   6.37   )             
 
Next, using the input values (6.33), we write eqs.(6.28-6.30) for j≧4 as follows. 
                 a   1   j     =       -     1   2       ⁢     a   1     j   -   1           ,           (   6.38   )                   a   2   j     =       -     1   2       ⁢     a   2     j   -   1           ,           (   6.39   )                   a   k   j     =           -     1   2       ⁢     a   k     j   -   1         -       1   2     ⁢     a     k   -   2       j   -   1       ⁢           ⁢   for   ⁢           ⁢   k       =   3       ,   …   ⁢           ,     j   -   1     ,           (   6.40   )                   a   j   j     =       -     a   j     j   -   1         -       1   2     ⁢     a     j   -   2       j   -   1             ,           (   6.41   )                 a     j   +   1     j     =       -     1   2       ⁢       a     j   -   1       j   -   1       .               (   6.42   )             
 
We use the values of a 1   3 , a 2   3 , a 3   3 , a 4   3  from eq.(6.37) to start these equations when j=4. We use eq.(6.36) to compute σ j . We start this equation with σ 3 =½, which follows from eqs.(6.35, 6.37).
 
   Compensated summation is based on the following principle. Let û and {circumflex over (υ)} denote real floating-point numbers, and let fl(û ±{circumflex over (υ)}) denote the floating-point sum and difference of û and {circumflex over (υ)}. Define the following floating-point numbers:
 
ŝ=fl(û+{circumflex over (υ)}),  (6.43) 
 
û′=fl(ŝ−{circumflex over (υ)}),  (6.44) 
 
{circumflex over (υ)}″=fl(ŝ−{circumflex over (υ)}û′),  (6.45) 
 
{circumflex over (δ)}û′(û)=fl(û−û′),  (6.46) 
 
{circumflex over (δ)}″({circumflex over (υ)})=fl({circumflex over (υ)}−{circumflex over (υ)}″),  (6.47) 
 
Δ(û,{circumflex over (υ)})=fl({circumflex over (δ)}′(û)+{circumflex over (δ)}″({circumflex over (υ)})).  (6.48) 
 
The quantity Δ(û,{circumflex over (υ)}) is a very good estimate of the roundoff error(û+{circumflex over (υ)})−ŝ. Note that the arithmetic operations in eqs.(6.43-6.48) carry over directly to complex floating-point arithmetic because the operations on the real parts of the complex numbers and the operations on the imaginary parts of these numbers are identical but independent.
 
   First, let us apply compensated summation to eq.(6.40). We want to compute 
                 a   k   j     =         -     (         1   2     ⁢     a   k     j   -   1         +       1   2     ⁢     a     k   -   2       j   -   1           )       ⁢           ⁢   for   ⁢           ⁢   k     =   3       ,   …   ⁢           ,     j   -   1     ,           (   6.49   )             
 
but we actually compute 
                 a   k   j     =         -     fl   ⁡     (         1   2     ⁢     a   k     j   -   1         +       1   2     ⁢     a     k   -   2       j   -   1           )         ⁢           ⁢   for   ⁢           ⁢   k     =   3       ,   …   ⁢           ,     j   -   1.             (   6.50   )             
 
Let 
                 δ   k   j     =       (         1   2     ⁢       a   ^     k     j   -   1         +       1   2     ⁢       a   ^       k   -   2       j   -   1           )     -     fl   ⁡     (         1   2     ⁢       a   ^     k     j   -   1         +       1   2     ⁢       a   ^       k   -   2       j   -   1           )           ⁢           ⁢     
     ⁢         for   ⁢           ⁢   k     =   3     ,   …   ⁢           ,     j   -   1     ,             (   6.51   )             
 
and
 
ε k   j =a k   j −â k   j  for k=1, . . . ,j+1.  (6.52) 
 
Now, eqs.(6.49-6.52) imply that 
                     a   k   j     =       -     (         1   2     ⁢     a   k     j   -   1         +       1   2     ⁢     a     k   -   2       j   -   1           )       +     δ   k   j                   =       -     (         1   2     ⁢     a   k     j   -   1         -       1   2     ⁢     ε   k     j   -   1         +       1   2     ⁢     a     k   -   2       j   -   1         -       1   2     ⁢     ε     k   -   2       j   -   1           )       +     δ   k   j                   =       -     (         1   2     ⁢     a   k     j   -   1         +       1   2     ⁢     a     k   -   2       j   -   1           )       +       1   2     ⁢     ε   k     j   -   1         +       1   2     ⁢     ε     k   -   2       j   -   1         +     δ   k   j                     =         a   k   j     +       1   2     ⁢     ε   k     j   -   1         +       1   2     ⁢     ε     k   -   2       j   -   1         +       δ   k   j     ⁢           ⁢   for   ⁢           ⁢   k       =   3       ,   …   ⁢           ,     j   -   1.                   (   6.53   )             
 
Hence, 
                 ε   k   j     =         (         1   2     ⁢     ε   k     j   -   1         +       1   2     ⁢     ε     k   -   2       j   -   1         +     δ   k   j       )     ⁢           ⁢   for   ⁢           ⁢   k     =   3       ,   …   ⁢           ,     j   -   1.             (   6.54   )             
 
In practice we compute an approximation to ε k   j  that we call {circumflex over (ε)} k   j : 
                   ε   ^     k   j     =         -     fl   ⁡     (         1   2     ⁢       ε   ^     k     j   -   1         +       1   2     ⁢       ε   ^       k   -   2       j   -   1         +     Δ   k   j       )         ⁢           ⁢   for   ⁢           ⁢   k     =   3       ,   …   ⁢           ,     j   -   1     ,           (   6.55   )             
 
where {circumflex over (Δ)} k   j  is the estimate for δ k   j  that we get with 
         u   ^     =       1   2     ⁢       a   ^     k     j   -   1             
 
and 
         v   ^     =       1   2     ⁢     a     k   -   2       j   -   1             
 
in eq. (6.48). Note that the RHS of eq.(6.55) is ambiguous because the associative law does not hold for floating-point addition. However, we shall ignore the small influence that the order of these additions may have on the numerical value of the total sum.
 
   Next, we apply compensated summation to eq.(6.41). We want to compute 
                 a   j   j     =     -     (       a   j     j   -   1       +       1   2     ⁢     a     j   -   2       j   -   1           )         ,           (   6.56   )             
 
but we actually compute 
                 a   ^     j   j     =       -   f     ⁢           ⁢       l   ⁡     (         a   ^     j     j   -   1       +       1   2     ⁢       a   ^       j   -   2       j   -   1           )       .               (   6.57   )             
 
Let 
               δ   j   j     =       (         a   ^     j     j   -   1       +       1   2     ⁢       a   ^       j   -   2       j   -   1           )     -       fl   ⁡     (         a   ^     j     j   -   1       +       1   2     ⁢       a   ^       j   -   2       j   -   1           )       .               (   6.58   )             
 
Now, eqs.(6.52,6.56-6.58) imply that 
                       a   ^     j   j     =       (         a   ^     j     j   -   1       +       1   2     ⁢       a   ^       j   -   2       j   -   1           )     +     δ   j   j                   =       -     (       a   j     j   -   1       -     ε   j     j   -   1       +       1   2     ⁢     a     j   -   2       j   -   1         -       1   2     ⁢     ε     j   -   2       j   -   1           )       +     δ   j   j                   =       -     (       a   j     j   -   1       +       1   2     ⁢     a     j   -   2       j   -   1           )       +     ε   j     j   -   1       +       1   2     ⁢     ε     j   -   2       j   -   1         +     δ   j   j                   =       a   j   j     +     ε   j     j   -   1       +       1   2     ⁢     ε     j   -   2       j   -   1         +       δ   j   j     .                     (   6.59   )             
 
Hence, 
               ε   j   j     =     -       (       ε   j     j   -   1       +       1   2     ⁢     ε     j   -   2       j   -   1         +     δ   j   j       )     .               (   6.60   )             
 
In practice we compute an approximation to ε j   j  that we call {circumflex over (ε)} j   j : 
                   ε   ^     j   j     =     -     fl   ⁡     (         ε   ^     j     j   -   1       +       1   2     ⁢       ε   ^       j   -   2       j   -   1         +       Δ   ^     j   j       )           ,           (   6.61   )             
 
where {circumflex over (Δ)} j   j  is the estimate for δ j   j  that we get with û=â j   j−1  and 
         v   ^     =       1   2     ⁢     a     j   -   2       j   -   1             
 
in eq.(6.48). The remarks that follow eq.(6.55) apply to eq. (6.61) as well.
 
   The treatment of eqs.(6.38,6.39,6.42) is simple because they involve no summations. Since we want to compute 
           a   1   j     =       -     1   2       ⁢     a   1     j   -   1           ,       a   2   j     =       -     1   2       ⁢     a   2     j   -   1           ,       a     j   +   1     j     =       -     1   2       ⁢     a     j   -   1       j   -   1           ,       
 
the computed values of these coefficients satisfy the same equations: 
                   a   ^     1   j     =       -     1   2       ⁢       a   ^     1     j   -   1           ,         a   ^     2   j     =       -     1   2       ⁢       a   ^     2     j   -   1           ,         a   ^       j   +   1     j     =       -     1   2       ⁢         a   ^       j   -   1       j   -   1       .                 (   6.62   )             
 
Now the errors ε 1   j , ε 2   j , ε j+1   j  also satisfy these equations: 
                 ε   1   j     =       -     1   2       ⁢     ε   1     j   -   1           ,       ε   2   j     =       -     1   2       ⁢     ε   2     j   -   1           ,       ε     j   +   1     j     =       -     1   2       ⁢       ε     j   -   1       j   -   1       .                 (   6.63   )             
 
Hence, we use the same equations to compute the approximations {circumflex over (ε)} 1   j , {circumflex over (ε)} 2   j , {circumflex over (ε)} j+1   j  to these errors: 
                   ε   ^     1   j     =       -     1   2       ⁢       ε   ^     1     j   -   1           ,         ε   ^     2   j     =       -     1   2       ⁢       ε   ^     2     j   -   1           ,         ε   ^       j   +   1     j     =       -     1   2       ⁢         ε   ^       j   -   1       j   -   1       .                 (   6.64   )             
 
   The approximate corrections {circumflex over (ε)} k   3  vanish because the computed coefficients â k   3  are exact. We use the values {circumflex over (ε)} 1   3 =0, {circumflex over (ε)} 2   3 =0, {circumflex over (ε)} 3   3 =0, {circumflex over (ε)} 4   3 =0 to start eqs.(6.55,6.61,6.64) when j=4. Since {circumflex over (ε)} 1   3 =0 and {circumflex over (ε)} 2   3 =0, the first two parts of eq.(6.64) imply that {circumflex over (ε)} 1   j =0 and {circumflex over (ε)} 2   j =0 for all j≧4. Therefore, we need only the last part of eq.(6.64). 
   Finally, we consider {circumflex over (σ)} j  and develop an approximate correction to it. We want to compute
 
σ j =−(a j+1   j +σ j−1 ) for j≧4,  (6.65) 
 
with σ 3 =½. In practice, we compute
 
{circumflex over (σ)} j =−fl(â j+1   j +{circumflex over (σ)} j−1 ) for j≧4,  (6.66) 
 
with {circumflex over (σ)} 3 =½. Define
 
e j =σ j −{circumflex over (σ)} j  for j≧3.  (6.67) 
 
By the same approach that lead to eqs.(6.54,6.60) we find that
 
e j =−(ε j+1   j +e j−1 +δ j ) for j≧4,  (6.68) 
 
where
 
δ j =â j+1   j +{circumflex over (σ)} j−1 −fl(â j+1   j +{circumflex over (σ)} j−1 ).  (6.69) 
 
Although we could use eq.(6.48) to estimate δ j  we shall neglect it against ε j+1   j . Therefore, we compute an approximate correction ê j  to {circumflex over (σ)} j  by the following recurrence formula. 
 
{circumflex over (ε)} j =−fl({circumflex over ( 68  )} j+1   j +ê j−1 ) for j≧4,  (6.70) 
 
where ê 3 =0.
 
   Let us summarize the computational forms of eqs.(6.36,6.38-6.42) for j≧4 that use compensated summation to identify roundoff errors and to cancel them to working precision. 
                           a   ^     1   j     =       -     1   2       ⁢       a   ^     1     j   -   1           ,                     a   ^     2   j     =       -     1   2       ⁢       a   ^     1     j   -   1           ,                   for   ⁢           ⁢   k     =   3     ,   …   ⁢           ,       j   -   1     :                       a   ^     k   j     =         -     fl   ⁡     (         1   2     ⁢       a   ^     k     j   -   1         +       1   2     ⁢       a   ^       k   -   2       j   -   1           )         ⁢           ⁢   and   ⁢           ⁢       ε   ^     k     -   j         =     -     fl   ⁡     (         1   2     ⁢       ε   ^     k     j   -   1         +       1   2     ⁢       ε   ^       k   -   2       j   -   1         +       Δ   ^     k   j       )             ,                     a   ^     j   j     =         -     fl   ⁡     (         a   ^     j     j   -   1       +       1   2     ⁢       a   ^       j   -   2       j   -   1           )         ⁢           ⁢   and   ⁢           ⁢       ε   ^     j     -   j         =     -     fl   ⁡     (         ε   ^     j     j   -   1       +       1   2     ⁢       ε   ^       j   -   2       j   -   1         +       Δ   ^     j   j       )             ,                     a   ^       j   +   1     j     =         1   2     ⁢       a   ^       j   -   1       j   -   1       ⁢           ⁢   and   ⁢           ⁢       ε   ^       j   +   1       -   j         =       -     1   2       ⁢       ε   ^       j   -   1       j   -   1             ,                   σ   ^     j     =         -     fl   ⁡     (         a   ^       j   +   1     j     +       σ   ^       j   -   1         )         ⁢           ⁢   and   ⁢           ⁢       e   ^     j       =     -       fl   ⁡     (         ε   ^       j   +   1       -   j       +       e   ^       j   -   1         )       .                 }           (   6.71   )             
 
   We approximate σ j  by the sum {circumflex over (σ)} j +ê j . 
   The number of arithmetic operations in eq. (6.28) of the simple algorithm is proportional to j−1, which increases with range. We can limit the size of this operation count by truncating the algorithm. First, we choose an integer L≧3, which will be called the length of the truncated algorithm. If j≧L, then we discard all the coefficients a k   j−1  such that k≦j−L. We effectively equate these coefficients to 0, which lets us start the recursion in eq.(6.28) for the truncated algorithm at k=j−L+1. Since the coefficient a j−L+1   j  is not used in subsequent steps, however, there is no need to evaluate it. For this reason, if j≧L, then we can start the recursion in eq.(6.28) for the truncated algorithm at k=j−L+2. Now the operation count is proportional to L−2, which is fixed. 
   We shall state the truncated algorithm in terms of new coefficients a l   j−1  that are defined for j=1,2, . . . and l=1, . . . ,L as follows. 
   (1) Let
 
a l   0 =0 for l=1, . . . ,L.  (6.72) 
 
   (2) Let
 
a 0   j =0 for j=0,1, . . .  (6.73) 
 
   (3) For j=1,2, . . . , 
                 α   l   j     =         -     1   2       ⁢     α     l   +   1       j   -   1         -       1   2     ⁢     α     l   -   1       j   -   1             ,       for   ⁢           ⁢   l     =   1     ,   …   ⁢           ,     L   -   2     ,           (   6.74   )                   α     L   -   1     j     =         1   2     ⁢       p   0     -   2       ⁡     [         ∂     ∂   y       ⁢       ψ     j   -   1       ⁡     (   0   )         +       ∂     ∂   y       ⁢       ψ   j     ⁡     (   0   )           ]         -     α   L     j   -   1       -       1   2     ⁢     α     L   -   2       j   -   1             ,           (   6.75   )                 α   L   j     =         1   2     ⁢       p   0     -   2       ⁡     [         ∂     ∂   y       ⁢       ψ     j   -   1       ⁡     (   0   )         +       ∂     ∂   y       ⁢       ψ   j     ⁡     (   0   )           ]         -       1   2     ⁢       α     L   -   1       j   -   1       .                 (   6.76   )             
 
The transparent boundary conditions associated with this algorithm are 
                     ∂     ∂   y       ⁢       ψ     j   +   1       ⁡     (   0   )         +       p   0     ⁢       ψ     j   +   1       (   L   )       ⁡     (   0   )           =           -     ∂     ∂   y         ⁢       ψ   j     ⁡     (   0   )         -       p   0     ⁢       ψ   j     (   L   )       ⁡     (   0   )         +     2   ⁢     p   0   2     ⁢     α   L   j     ⁢           ⁢   for   ⁢           ⁢   j       =   0       ,   1   ,   …           (   6.77   )             
 
The superscript “L” indicates that the sequence of boundary values {ψ j   (L) (0)} j=0   ∞  is produced from the sequence of boundary derivatives {σψ j (0)/σy} j=0   ∞  by the truncated algorithm of length L. These boundary values are only approximate.
 
   An alternative form of these conditions is obtained by iterating eq.(6.77) backwards. The initial conditions ψ 0   (L) (0)=0 and σψ 0 (0)/σy=0 are used to stop the iteration. The result is like eq.(6.32). 
                     ∂     ∂   y       ⁢       ψ     j   +   1       ⁡     (   0   )         +       p   0     ⁢       ψ     j   +   1       (   L   )       ⁡     (   0   )           =       2   ⁢     p   0   2     ⁢       ∑     k   =   0     j     ⁢         (     -   1     )       j   -   k       ⁢     α   L   k     ⁢           ⁢   for   ⁢           ⁢   j         =   0       ,   1   ,   …           (   6.78   )             
 
   Suppose the sequence of differential boundary values σψ j (0)/σy in eq.(6.33) is used as input to the truncated algorithm. The errors in the approximate boundary values ψ j   (L) (0) relative to the exact boundary values ψ j (0) can be characterized as follows. Fix ξ&gt;0, and let {j 1 , j 2 , . . . j l , . . . } be a sequence of integers such that j l ˜ξl 2  as l→∞. Then the relative error in ψ jL   (L) (0) converges to a finite limit as L→∞ that is exponentially small if ξ is small. 
             0   &lt;       lim     L   →   ∞       ⁢     (     1   -         ψ   jL     (   L   )       ⁡     (   0   )       /       ψ   jL     ⁡     (   0   )           )       &lt;     2   ⁢       c       -   2     /   ɛ       .               (   6.79   )             
 
   Suppose the index j and the length L are given, and let ξ=jL −2 . A practical way to interpret eq.(6.79) is that 1−ψ j   (L) (0)/ψ j (0)≈2 exp(−2/ξ) if L &gt;&gt;1. 
   Let u denote the unit roundoff in a finite-precision floating-point computer arithmetic. Two real numbers a and b are represented by the same floating-point number if |1−a/b|&lt;u. These numbers are equivalent in this computer arithmetic. If L&gt;[ln(2/u)j/2] 1/2 , then 2 exp(−2/ξ)&lt;u. Hence, ψ j   (L) (0) and ψ j (0) are equivalent in the computer arithmetic if L&gt;&gt;1. 
   Now let J be an integer, and let L J  be the smallest integer greater than [ln(2/u)J/2] 1/2 .
 
L j =[|ln(2/u)J/2) 1/2 ]=ο(J 1/2 ).  (6.80) 
 
If L J &gt;&gt;1, then ψ j   (L     J     ) (0) and ψ j (0) are equivalent in the computer arithmetic for all j≦J. The number of arithmetic operations to compute each ψ j   (L     J     ) (0) is proportional to L J , but the number of arithmetic operations to compute each ψ j (0) is proportional to j. The ratio ρ(J) of the total operation count for computing the subsequence {ψ j   (L     J     ) (0)} j=1   J  to the total operation count for computing the subsequence {ψ j (0)} j=1   J  is approximately 
                 ρ   ⁡     (   J   )       ≈     J   ⁢           ⁢       L   J     /     1   2       ⁢     J   2         =       2   ⁢       L   J     /   J       =       O   ⁡     (     J       -   1     /   2       )       .               (   6.81   )             
 
   Therefore, if L J &gt;&gt;1, then the truncated algorithm of length L J  is more efficient for computing the first J boundary values to working precision than the complete algorithm. As an example of practical importance, consider IEEE double-precision floating-point arithmetic. Since u=2 −53  in this arithmetic,
 
L J ≈4.326J 1/2  for IEEE double-precision floating-point arithmetic.  (6.82) 
 
These criteria are valid for all sequences of differential boundary values σψ j (0)/σy. However, let us continue to use the sequence in eq.(6.33) to explain how to program the truncated algorithm with compensated summation.
 
   Assume that L&gt;4. We can calculate the coefficients a l   j  for j=1,2,3 exactly. 
                             α   l   1     =       0   ⁢           ⁢   for   ⁢           ⁢   l     ≤     L   -   2                 α     L   -   1     1     =   1               α   L   1     =     -   1       ,                             α   l   2     =       0   ⁢           ⁢   for   ⁢           ⁢   l     ≤     L   -   3                 α     L   -   2     2     =     -     1   2                 α     L   -   1     2     =   2               α   L   2     =     -     3   2         ,                             α   l   3     =       0   ⁢           ⁢   for   ⁢           ⁢   l     ≤     L   -   4                 α     L   -   3     3     =     1   4               α     L   -   2     3     =     -   1               α     L   -   1     3     =         7   4     ⁢           ⁢     α   L   3       =     -   1.                     }           (   6.83   )             
 
Next, we use eqs.(6.33,6.73) to write eqs.(6.74-6.76) for j≧4 as follows. 
                 α   1   j     =       -     1   2       ⁢     α   2     j   -   1           ,           (   6.84   )                   α   l   j     =           -     1   2       ⁢     α     l   +   1       j   -   1         -       1   2     ⁢     α     l   -   1       j   -   1       ⁢           ⁢   for   ⁢           ⁢   l       =   2       ,   …   ⁢           ,     L   -   2     ,           (   6.85   )                   α     L   -   1     j     =       -     α   L     j   -   1         -       1   2     ⁢     α     L   -   2       j   -   1             ,           (   6.86   )                 α   L   j     =       -     1   2       ⁢       α     L   -   1       j   -   1       .               (   6.87   )             
 
   Since a L   0 =0, it follows from eq.(6.78) that ψ 1   (L) (0)=−2p 0  and 
                   ψ     j   +   1       (   L   )       ⁡     (   0   )       =       2   ⁢           ⁢     p   0     ⁢       ∑     k   =   1     j     ⁢         (     -   1     )       j   -   k       ⁢     α   L   k     ⁢           ⁢   for   ⁢           ⁢   j         =   1       ,   2   ,     …   ⁢           .             (   6.88   )             
 
Thus, ψ j+1   (L) (0)=−2p 0 σ j   (L)  where 
                 σ   j     (   L   )       =       -       ∑     k   =   1     j     ⁢         (     -   1     )       j   -   k       ⁢     α   L   k     ⁢           ⁢   for   ⁢           ⁢   j         =   1       ,   2   ,   …   ⁢           ,           (   6.89   )             
 
which can be calculated recursively as follows.
 
σ j   (L) =a L   j −σ j−1   (L)  for j=2,3, . . .  (6.90) 
 
We use the values of a l   3  for l=1, . . . ,L from eq.(6.83) to start eqs.(6.84-6.87) when j=4. We compute σ j   (L)  from eq.(6.90). We start this equation with σ 3   (L) ={fraction ( 1 / 2 )}, which follows from eqs.(6.83,6.89).
 
   If eqs.(6.84-6.87,6.90) are programmed directly, then roundoff errors tend to accumulate. We can use compensated summation to counteract these roundoff errors in the same way that we use it to counteract roundoff errors for the complete algorithm. Let us denote the computed values of a k   j  and σ j   (L)  by â k   j  and {circumflex over (σ)} j   (L) , respectively. We list the corrected computational forms of eqs.(6.84-6.87,6.90) next. If j≧4, then 
                       α   ^     1   j     =         -     1   2       ⁢       α   ^     2     j   -   1       ⁢           ⁢   and   ⁢           ⁢       ε   ^     1   j       =       -     1   2       ⁢       ε   ^     2     j   -   1             ,                   for   ⁢           ⁢   l     =   2     ,   …   ⁢           ,     L   -     2   ⁢     :                           α   ^     l   j     =         -     fl   ⁡     (         1   2     ⁢       α   ^       l   +   1       j   -   1         +       1   2     ⁢       α   ^       l   -   1       j   -   1           )         ⁢           ⁢   and   ⁢           ⁢       ε   ^     l   j       =     -     fl   ⁡     (         1   2     ⁢       ε   ^       l   +   1       j   -   1         +       1   2     ⁢       ε   ^       l   -   1       j   -   1         +       Δ   ^     l   j       )             ,                     α   ^       L   -   1     j     =         -     fl   ⁡     (         α   ^     L     j   -   1       +       1   2     ⁢       α   ^       L   -   2       j   -   1           )         ⁢           ⁢   and   ⁢           ⁢       ε   ^       L   -   1     j       =     -     fl   ⁡     (         ε   ^     L     j   -   1       +       1   2     ⁢       ε   ^       L   -   2       j   -   1         +       Δ   ^       L   -   1     j       )             ,                     α   ^     L   j     =         -     1   2       ⁢       α   ^       L   -   1       j   -   1       ⁢           ⁢   and   ⁢           ⁢       ε   ^     L   j       =       -     1   2       ⁢       ε   ^       L   -   1       j   -   1             ,                   σ   ^     j     (   L   )       =         -     fl   ⁡     (         α   ^     L   j     +       σ   ^       j   -   1       (   L   )         )         ⁢           ⁢   and   ⁢           ⁢       ε   ^     j     (   L   )         =     -       fl   ⁡     (         ε   ^     L   j     +       ε   ^       j   -   1       (   L   )         )       .                 }     ⁢           ⁢     (   6.91   )         
 
Using eqs.(6.43-6.48),we compute {circumflex over (Δ)} l   j  for l=2, . . . ,L−2 as {circumflex over (Δ)}(û,{circumflex over (υ)}) with 
         u   ^     =       1   2     ⁢       α   ^       l   +   1       j   -   1             
 
and 
           υ   ^     =       1   2     ⁢       α   ^       l   -   1       j   -   1           ,       
 
and we compute {circumflex over (Δ)} L−1   j  as {circumflex over (Δ)}(û,{circumflex over (υ)}) with û=â L   j−1  and 
         υ   ^     =       1   2     ⁢         α   ^       L   -   2       j   -   1       .           
 
   These equations are developed like the ones that are collected as eq.(6.71). We start them when j=4 with the following initial values: 
                       α   ^     l   3     =       ⁢       0   ⁢           ⁢   for   ⁢           ⁢   l     ≤     L   -   4                       α   ^       L   -   3     3     =       ⁢     1   4                     α   ^       L   -   2     3     =       ⁢     -   1                     α   ^       L   -   1     3     =       ⁢     7   4                     α   ^     L   3     =       ⁢     -   1                       σ   ^     3     (   L   )       =       ⁢     1   2       ,                   ε   ^     l   3     =       ⁢       0   ⁢           ⁢   for   ⁢           ⁢   l     ≤   L                     e   ^     3     (   L   )       =       ⁢   0.                 (   6.92   )             
 
We approximate σ j   (L)  by the sum {circumflex over (σ)} j   (L) , which we denote by {circumflex over (σ)} j   (L).  
 
7. A Chebyshev tau method for semi-discrete parabolic equations.
 
   If x j =jΔx and k(x j )=k 0  for j=0,1, . . . , where the increment Δx and the wavenumber k 0  are fixed, then eq.(6.9) takes the following simple form. 
                   ψ     j   +   1       ⁡     (   y   )       -       ψ   j     ⁡     (   y   )         =           (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4     ⁢     k   0       )     ⁡     [           ∂   2       ∂     y   2         ⁢       ψ   j     ⁡     (   y   )         +         ∂   2       ∂     y   2         ⁢       ψ     j   +   1       ⁡     (   y   )           ]       ⁢           ⁢   for   ⁢           ⁢   j     =           0   ,   1   ,   …                 (   7.1   )             
 
We can use Chebyshev interpolation polynomials to integrate this equation numerically over an interval −W≦y≦W. Let T k (ξ) denote the k-th Chebyshev polynomial (of the first kind). For k=0,1, . . . these polynomials are defined by the identity
 
T k (ξ)=cos(kθ), where ξ=cos θ.  (7.2) 
 
Let N be a positive integer. A Chebyshev interpolation polynomial of degree N over the interval −W≦y≦W is a sum of the form 
               f   ⁡     (   y   )       =       ∑     k   =   0     N     ⁢           ⁢       f   k     ⁢         T   k     ⁡     (     y   /   W     )       .                 (   7.3   )             
 
   We shall approximate each function ψ j (y) by a polynomial ƒ j (y) of degree N, which we represent uniquely as a Chebyshev interpolation polynomial. 
                 f   j     ⁡     (   y   )       =           ∑     k   =   0     N     ⁢           ⁢       f     j   ,   k       ⁢       T   k     ⁡     (     y   /   W     )       ⁢           ⁢   for       -   W     ≤   y   ≤     W   .               (   7.4   )             
 
The first and second derivatives of ƒ j (y) with respect to y are polynomials of degrees N−1 and N−2. We also represent these polynomials uniquely as Chebyshev interpolation polynomials. 
                   ∂     ∂   y       ⁢       f   j     ⁡     (   y   )         =           1   W     ⁢       ∑     k   =   0       N   -   1       ⁢           ⁢       f     j   ,   k       (   1   )       ⁢       T   k     ⁡     (     y   /   W     )       ⁢           ⁢   for         -   W     ≤   y   ≤   W       ,           (   7.5   )             
                   ∂   2       ∂     y   2         ⁢       f   j     ⁡     (   y   )         =           1     W   2       ⁢       ∑     k   =   0       N   -   2       ⁢           ⁢       f     j   ,   k       (   2   )       ⁢       T   k     ⁡     (     y   /   W     )       ⁢           ⁢   for         -   W     ≤   y   ≤     W   .               (   7.6   )             
 
   Assume that the coefficients in these equations satisfy the following Chebyshev tau equations. 
                 f       j   +   1     ,   k       -     f     j   ,   k         =       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4     ⁢     k   0     ⁢     W   2       )     ⁢     (       f       j   +   1     ,   k       (   2   )       +     f     j   ,   k       (   2   )         )     ⁢           ⁢   for   ⁢     {             j   =   0     ,   1   ,     …   ⁢           ⁢   and                   k   =   0     ,   1   ,   …   ⁢           ,     N   -   2.                         (   7.7   )             
 
Under this assumption, the polynomials ƒ j (y) satisfy equations that are similar to eq.(7.1). For j=0,1, . . . and −W≦y≦W, 
                     f     j   +   1       ⁡     (   y   )       -       f   j     ⁡     (   y   )         =         (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4     ⁢     k   0       )     ⁡     [           ∂   2       ∂     y   2         ⁢       f     j   +   1       ⁡     (   y   )         +         ∂   2       ∂     y   2         ⁢       f   j     ⁡     (   y   )           ]       +       e       j   +   1     ,     N   -   1         ⁢       T     N   -   1       ⁡     (     y   /   W     )         +       e       j   +   1     ,   N       ⁢       T   N     ⁡     (     y   /   W     )             ,           (   7.8   )             
 
where
 
e j+1,N−1 =ƒ j+1,N−1  and e j+1,N =ƒ j+1,N −ƒ j,N .  (7.9) 
 
   We also define the following dimensionless and real-valued parameter.
 
q=4k 0 W 2 /Δz.  (7.10) 
 
Now the Chebyshev tau equations can be written in the following equivalent form. 
                 f       j   +   1     ,   k       (   2   )       +     iq   ⁢           ⁢     f       j   +   1     ,   k           =       g     j   ,   k       ⁢           ⁢   for   ⁢     {             j   =   0     ,   1   ,     …   ⁢           ⁢   and                   k   =   0     ,   1   ,   …   ⁢           ,     N   -   2     ,                       (   7.11   )             
 
where
 
g j,k =−ƒ j,k   (2) +iqƒ j,k .  (7.12) 
 
The quantities g j,k  can be computed recursively without the coefficients ƒ j+1,k   (2) . 
                     g       j   +   1     ,   k       =       ⁢       -     f       j   +   1     ,   k       (   2   )         +     iq   ⁢           ⁢     f       j   +   1     ,   k                       =       ⁢       -     (       f       j   +   1     ,   k       (   2   )       +     iq   ⁢           ⁢     f       j   +   1     ,   k           )       +     2   ⁢   iq   ⁢           ⁢     f       j   +   1     ,   k                       =       ⁢       -     g     j   ,   k         +     2   ⁢   iq   ⁢           ⁢       f       j   +   1     ,   k       .                       (   7.13   )             
 
   It is possible also to eliminate the coefficients ƒ j+1,k   (2)  from the Chebyshev tau equations. This leads to an alternative formulation of these equations in terms of the coefficients ƒ j+1,k  alone. The coefficients ƒ j+1,k  and ƒ j+1,k   (2)  in eqs.(7.4,7.6) satisfy eq.(7.7) if and only if. 
                   f       j   +   1     ,   N       +       1   4     ⁢     1     N   ⁡     (     N   -   1     )         ⁢   iq   ⁢           ⁢     f       j   +   1     ,     N   -   2             =     h     j   ,   N         ,           (   7.14   )                     f       j   +   1     ,     N   -   1         +       1   4     ⁢     1       (     N   -   1     )     ⁢     (     N   -   2     )         ⁢   iq   ⁢           ⁢     f       j   +   1     ,     N   -   3             =     h     j   ,     N   -   1           ,           (   7.15   )                       (     1   -       1   2     ⁢     1       (     N   -   1     )     ⁢     (     N   -   3     )         ⁢   iq       )     ⁢     f       j   +   1     ,     N   -   2           +       1   4     ⁢     1       (     N   -   2     )     ⁢     (     N   -   3     )         ⁢   iq   ⁢           ⁢     f       j   +   1     ,     N   -   4             =     h     j   ,     N   -   2           ,           (   7.16   )                       (     1   -       1   2     ⁢     1       (     N   -   2     )     ⁢     (     N   -   4     )         ⁢   iq       )     ⁢     f       j   +   1     ,     N   -   3           +       1   4     ⁢     1       (     N   -   3     )     ⁢     (     N   -   4     )         ⁢   iq   ⁢           ⁢     f       j   +   1     ,     N   -   5             =     h     j   ,     N   -   3           ,           (   7.17   )                       1   4     ⁢     1       (     k   +   1     )     ⁢   k       ⁢   iq   ⁢           ⁢     f       j   +   1     ,     k   +   2           +       (     1   -       1   2     ⁢     1       k   2     -   1       ⁢   iq       )     ⁢     f       j   +   1     ,   k         +       1   4     ⁢     1     k   ⁡     (     k   -   1     )         ⁢   iq   ⁢           ⁢     f       j   +   1     ,     k   -   2             =         h     j   ,   k       ⁢           ⁢   for   ⁢           ⁢   k     =     N   -   4         ,   …   ⁢           ,   3   ,           (   7.18   )             
 for k=N−4, . . . ,3,  (7.18) 
                     1   24     ⁢   iq   ⁢           ⁢     f       j   +   1     ,   4         +       (     1   -       1   6     ⁢   iq       )     ⁢     f       j   +   1     ,   2         +       1   4     ⁢   iq   ⁢           ⁢     f       j   +   1     ,   0           =     h     j   ,   2         ,           (   7.19   )             
 
for j=0,1, . . . , where 
                 h     j   ,   N       =       1   4     ⁢     1     N   ⁡     (     N   -   1     )         ⁢     g     j   ,     N   -   2             ,           (   7.20   )                 h     j   ,     N   -   1         =       1   4     ⁢     1       (     N   -   1     )     ⁢     (     N   -   2     )         ⁢     g     j   ,     N   -   3     ,                 (   7.21   )                 h     j   ,     N   -   2         =         -     1   2       ⁢     1       (     N   -   1     )     ⁢     (     N   -   3     )         ⁢     g     j   ,     N   -   2           +       1   4     ⁢     1       (     N   -   2     )     ⁢     (     N   -   3     )         ⁢     g     j   ,     N   -   4     ,                   (   7.22   )                 h     j   ,     N   -   3         =         -     1   2       ⁢     1       (     N   -   2     )     ⁢     (     N   -   4     )         ⁢     g     j   ,     N   -   3           +       1   4     ⁢     1       (     N   -   3     )     ⁢     (     N   -   4     )         ⁢     g     j   ,     N   -   5     ,                   (   7.23   )                   h     j   ,   k       =           1   4     ⁢     1       (     k   +   1     )     ⁢   k       ⁢     g     j   ,     k   +   2           -       1   2     ⁢     1       k   2     -   1       ⁢     g     j   ,   k         +       1   4     ⁢     1     k   ⁡     (     k   -   1     )         ⁢     g     j   ,     k   -   2         ⁢           ⁢   for   ⁢           ⁢   k       =     N   -   4         ,   …   ⁢           ,   3   ,           (   7.24   )                 h     j   ,   2       =         1   24     ⁢     g     j   ,   4         -       1   6     ⁢     g     j   ,   2         +       1   4     ⁢       g     j   ,   0       .                 (   7.25   )             
 
From this point on we assume that N is even. Note that
         (1) there are N−1 equations in eqs.(7.14-7.19); ,   (2) the even-indexed coefficients ƒ j+1,N , ƒ j+1,N−2 , . . . ƒ j+1 , ƒ j+1,0  are uncoupled from the odd-indexed coefficients ƒ j+1,N−1 , ƒ j+1,N−3 , . . . ,ƒ j+1,1  in eqs.(7.14-7.19);   (3) the coefficient f j+1,N  appears only in eq.(7.14), and the coefficient ƒ j+1,N−1  appears only in eq.(7.15).       

   We adopt conditions at the boundaries y=±W that combine the discrete transparent boundary conditions (6.32) with source terms −ƒ j+1   (S) (±W) that are independent of the solutions f k (y) for k≦j+1. Thus, 
                     f     j   +   1       ⁡     (     ±   W     )       ±       1     p   0       ⁢     ∂     ∂   y       ⁢       f     j   +   1       ⁡     (     ±   W     )           =           g   j     ⁡     (     ±   W     )       ⁢           ⁢   for   ⁢           ⁢   j     =   0       ,   1   ,   …   ⁢           ,           (   7.26   )             
 j=0,1, . . . ,  (7,26)
 
where 
                 g   j     ⁡     (     ±   W     )       =       -       f     j   +   1       (   S   )       ⁡     (     ±   W     )         +     2   ⁢     p   0     ⁢       ∑     k   =   0     j     ⁢         (     -   1     )       j   -   k       ⁢         a     k   +   1     k     ⁡     (     ±   W     )       .                     (   7.27   )             
 
The coefficients a k+1   k  (±W) in this equation are computed recursively. First, a 1   0 (±W)=0. Second, for j=1,2, . . . , 
                   a   k   j     ⁡     (     ±   W     )       =         -     1   2       ⁢       a   k     j   -   1       ⁡     (     ±   W     )         -       1   2     ⁢       a     k   -   2       j   -   1       ⁡     (     ±   W     )             ,     for   ⁢           ⁢         k   =   1     ,   …   ⁢           ,     j   -   1         ︸       Skip   ⁢           ⁢   this   ⁢           ⁢   when   ⁢           ⁢   j     =   1.           ,           (   7.28   )                     a   j   j     ⁡     (     ±   W     )       =         ±           ⁢     1   2       ⁢       p   0     -   2       ⁡     [         ∂     ∂   y       ⁢       f     j   -   1       ⁡     (     ±   W     )         +       ∂     ∂   y       ⁢       f   j     ⁡     (     ±   W     )           ]         -       a   j     j   -   1       ⁡     (     ±   W     )       -       1   2     ⁢       a     j   -   2       j   -   1       ⁡     (     ±   W     )             ,           (   7.29   )                     a     j   +   1     j     ⁡     (     ±   W     )       =         ∓     1   2       ⁢       p   0     -   2       ⁡     [         ∂     ∂   y       ⁢       f     j   -   1       ⁡     (     ±   W     )         +       ∂     ∂   y       ⁢       f   j     ⁡     (     ±   W     )           ]         -       1   2     ⁢       a     j   -   1       j   -   1       ⁡     (     ±   W     )             ,           (   7.30   )             
 
where
 
a −1   j (±W)=0 for j=0,1, . . . ,  (7.31) 
 
a 0   j (±W)=0 for j=0,1, . . .  (7.32) 
 
If ƒ 0 (y) is specified, then the N+1 equations in eqs.(7.14-7.19,7.26)) determine the polynomials ƒ j+1 (y) for j=0,1, . . . uniquely.
 
   It is easy to evaluate ƒ j+1 (±W) from eq.(7.4) because T k (±1)=(±1) k . 
                 f     j   +   1       ⁡     (     ±   W     )       =       ∑     k   =   0     N     ⁢         (     ±   1     )     k     ⁢       f       j   +   1     ,   k       .                 (   7.33   )             
 
It is more convenient to evaluate σƒ j+1 (±W)/σy by differentiating eq.(7.4) directly than by using eq.(7.5) because T k ′(±1)=(±1) k+1 k 2 . Thus, 
                   ±     ∂     ∂   y         ⁢     f     j   +   1         ±     (   W   )       =       1   W     ⁢       ∑     k   =   0     N     ⁢         (     ±   1     )     k     ⁢     k   2     ⁢       f       j   +   1     ,   k       .                   (   7.34   )             
 
Now
 
p 0 W=e −iπ/4 2W(k 0 Δx) 1/2 =e −iπ/4 q 1/2 .  (7.35) 
 
Therefore, it follows from eqs.(7.26,7.33-7.35) that 
                     ∑     k   =   0     N     ⁢     f       j   +   1     ,   k         +       ⅇ     ⅈ   ⁢           ⁢     π   /   4         ⁢     q       -   1     /   2       ⁢       ∑     k   =   0     N     ⁢       k   2     ⁢     f       j   +   1     ,   k               =       g   j     ⁡     (   W   )         ,           (   7.36   )                     ∑     k   =   0     N     ⁢         (     -   1     )     k     ⁢     f       j   +   1     ,   k           +       ⅇ     ⅈ   ⁢           ⁢     π   /   4         ⁢     q       -   1     /   2       ⁢       ∑     k   =   0     N     ⁢         (     -   1     )     k     ⁢     k   2     ⁢     f       j   +   1     ,   k               =         g   j     ⁡     (     -   W     )       .             (   7.37   )             
 
   It is easy to decouple the even-indexed Chebyshev coefficients from the odd-indexed Chebyshev coefficients in these equations. If we add eqs.(7.36,7.37), then we get an equation that involves only the even-indexed Chebyshev coefficients. 
                 ∑     l   =   0       N   /   2       ⁢       [     1   +       ⅇ     ⅈ   ⁢           ⁢     π   /   4         ⁢         q       -   1     /   2       ⁡     (     2   ⁢   l     )       2         ]     ⁢     f       j   +   1     ,     2   ⁢   l             =         1   2     ⁡     [         g   j     ⁡     (   W   )       +       g   j     ⁡     (   W   )         ]       .             (   7.38   )             
 
If we subtract eq.(7.37) from eq.(7.36), then we get an equation that involves only the odd-indexed Chebyshev coefficients. 
                 ∑     l   =   1       N   /   2       ⁢       [     1   +       ⅇ     ⅈ   ⁢           ⁢     π   /   4         ⁢         q       -   1     /   2       ⁡     (       2   ⁢   l     -   1     )       2         ]     ⁢     f       j   +   1     ,       2   ⁢   l     -   1             =         1   2     ⁡     [         g   j     ⁡     (   W   )       -       g   j     ⁡     (     -   W     )         ]       .             (   7.39   )             
 
   It is useful to express eqs.(7.14-7.19) and eqs.(7.38,7.39) as a pair of independent matrix equations of the form 
                 [         A       u             v   T         γ         ]     ⁡     [         x           α         ]       =       [         y           β         ]     .             (   7.40   )             
 
First, define 
               τ   k     =     {         0             for   ⁢           ⁢   k     =     N   -   2       ,     N   -   3     ,                 1   4     ⁢     1       (     k   +   1     )     ⁢   k       ⁢   iq               for   ⁢           ⁢   k     =     N   -   4       ,   …   ⁢           ,   2.                     (   7.41   )                 s   k     =     {         1             for   ⁢           ⁢   k     =   N     ,     N   -   1     ,               1   -       1   2     ⁢     1       k   2     -   1       ⁢   iq                 for   ⁢           ⁢   k     =     N   -   2       ,   …   ⁢           ,   2.                     (   7.42   )                 t   k     =     {             1   4     ⁢     1     k   ⁡     (     k   -   1     )         ⁢   iq               for   ⁢           ⁢   k     =   N     ,   …   ⁢           ,   3   ,                 1   4     ⁢   iq             for   ⁢           ⁢   k     =   2.                     (   7.43   )             
 b k =1+e ix/4 q −1/2 k 2  for k=N, . . . ,0.  (7.44)
 
The even-indexed Chebyshev coefficients ƒ j+1,2l  are determined by the matrix equation 
                   [           A   0           u   0               v   0   T           γ   0           ]     ⁢           [           x   0               α   0           ]     =     [           y   0               β   0           ]       ,           (   7.45   )             
 
where 
                 A   0     =     [           s   N           t   N                                                                           T     N   -   2             s     N   -   2             t     N   -   2                                                                               T     N   -   4       ⁢                   s     N   -   4             t     N   -   4                                                                           ⋰       ⋰       ⋰                                                                       ⋰       ⋰       ⋰                                                                         r   4           s   4           t   4                                                                           r   2           s   2           ]       ,           ⁢       u   0     =     [         0           0           0           ⋮           ⋮           0             t   2           ]       ,           (   7.46   )                   v   0     =     [           b   N               b     N   -   2                 b     N   -   4               ⋮           ⋮             b   4               b   2           ]       ,           ⁢       x   0     =     [           f       j   +   1     ,   N                 f       j   +   1     ,     N   -   2                   f       j   +   1     ,     N   -   4                 ⋮           ⋮             f       j   +   1     ,   4                 f       j   +   1     ,   2             ]       ,           ⁢       y   0     =     [           h     j   ,   N                 h     j   ,     N   -   2                   h     j   ,     N   -   4                 ⋮           ⋮             h     j   ,   4                 h     j   ,   2             ]       ,           (   7.47   )             
 
The rows of eq.(7.45) are eqs.(7.14,7.16), eq.(7.18) for k=2l (l=N/2−2, . . . ,2), eq.(7.19) and eq.(7.38).
 
   The odd-indexed Chebyshev coefficients ƒ j+1,2l−1  are determined by the matrix equation 
                   [           A   1           u   1               v   1   T           γ   1           ]     ⁢           [           x   1               α   1           ]     =     [           y   1               β   1           ]       ,           (   7.49   )             
 
where 
                 A   1     =     [           s     N   -   1             t     N   -   1                                                                             τ     N   -   3             s     N   -   3             t     N   -   3                                                                             τ     N   -   5             s     N   -   5             t     N   -   5                                                                           ⋰       ⋰       ⋰                                                                       ⋰       ⋰       ⋰                                                                         τ   5           s   5           t   5                                                                           τ   3           s   3           ]       ,       u   1     =     [         0           0           0           ⋮           ⋮           0             t   3           ]       ,           (   7.50   )                   v   1     =     [           b     N   -   1                 b     N   -   3                 b     N   -   5               ⋮           ⋮             b   5               b   3           ]       ,       x   1     ⁡     [           f       j   +   1     ,     N   -   1                   f       j   +   1     ,     N   -   3                   f       j   +   1     ,     N   -   5                 ⋮           ⋮             f       j   +   1     ,   5                 f       j   +   1     ,   3             ]       ,       y   1     =     [           h     j   ,     N   -   1                   h     j   ,     N   -   3                   h     j   ,     N   -   5                 ⋮           ⋮             h     j   ,   5                 h     j   ,   3             ]       ,           (   7.51   )             
 γ 1 =b 1 , a 1 =ƒ j+1,1 , β 1 =½|g j (W)−g j (−W)|.  (7.52)
 
The rows of eq.(7.49) are eqs.(7.15,7.17), eq. (7.18) for k=2l−1(l=N/2−2, . . . ,2) and eq.(7.39).
 
   We can solve eqs.(7.45,7.49) numerically by the Crout form of the block elimination method. The BEC algorithm to carry out this method has five steps. We outline them for eq.(7.40), which includes eqs.(7.45,7.49) as particular cases. Given y and β, we want to solve that equation for x and a. 
   (1) Solve the equation
 
Aξ=u.  (7.53) 
 
   (2) Use the solution vector ξ to compute the divisor
 
δ=γ−v T ·ξ.  (7.54) 
 
   (3) Solve the equation
 
Aη=y.  (7.55) 
 
   (4) Use the solution vector η to compute the scalar
 
a=(β−v T ·η)/δ.  (7.56) 
 
   (5) Compute the vector
 
x=η−aξ.  (7.57) 
 
We can solve eqs.(7.53,7.55) with a method that uses an LU factorization of the matrix A. If A=A 0  or A=A 1 , then the Crout LU factorization can be done numerically without pivoting and is stable.
 
   Since s N =1, we find that
 
A 0 =L 0 U 0 ,  (7.58) 
 
where 
                 L   0     =     [         1                                                                                     σ     N   -   2                                                                             τ     N   -   4             σ     N   -   4                                                                           ⋰       ⋰                                                                         τ   4           σ   4                                                                           τ   2           σ   2           ]       ,     
     ⁢       U   0     =     [         1         t   N                                                                         1         τ     N   -   2                                                                           1         τ     N   -   4                                                                           ⋰       ⋰                                                                       1         τ   4                                                                         1         ]       ,           (   7.59   )             
 
and 
                               σ     N   -   2       =     s     N   -   2         ,               τ     N   -   2       =       t     N   -   2       /     σ     N   -   2           ,                               σ   k     =       s   k     -       τ   k     ⁢     τ     k   +   2             ,             τ   k     =       t   k     /     σ   k                   for   ⁢           ⁢   k     =     2   ⁢           ⁢     l   ⁡     (       l   =       N   /   2     -   2       ,   …   ⁢           ,   2     )           ,                       σ   2     =       s   2     -       τ   2     ⁢       τ   4     .                 }           (   7.60   )             
 
Since s N−1 =1, we find that
 
A 1 =L 1 U 1 ,  (7.61) 
 
where 
                 L   1     =     [         1                                                                                     σ     N   -   3                                                                             τ     N   -   5             σ     N   -   5                                                                           ⋰       ⋰                                                                         τ   5           σ   5                                                                           τ   3           σ   3           ]       ,     
     ⁢       U   1     =     [         1         t     N   -   1                                                                           1         τ     N   -   3                                                                           1         τ     N   -   5                                                                           ⋰       ⋰                                                                       1         τ   5                                                                         1         ]       ,           (   7.62   )             
 
and 
                               σ     N   -   3       =     s     N   -   3         ,               τ     N   -   3       =       t     N   -   3       /     σ     N   -   3           ,                         σ   k     =       s   k     -       τ   k     ⁢     τ     k   +   2             ,       τ   k     =         t   k     /     σ   k       ⁢           ⁢   for                     k   =       2   ⁢           ⁢   l     -     1   ⁢     (       l   =       N   /   2     -   2       ,   …   ⁢           ,   3     )           ,                 σ   3     =       s   3     -       τ   3     ⁢       τ   5     .                 }           (   7.63   )             
 
   If A=A 0 , then eqs.(7.53,7.55) are of the form A 0 z 0 =w 0 . To solve this equation we introduce a vector p 0  and solve the first intermediate equation L 0 p 0 =w 0  for p 0  by forward substitution. Then we solve the second intermediate equation U 0 z 0 =p 0  for z 0  by backward substitution. 
   If A=A 1 , then eqs.(7.53,7.55) are of the form A 1 z 1 =w 1 . This equation is solved the same way. We introduce a vector p 1  and solve the intermediate equation L 1 p 1 =w 1  for p 1  by forward substitution. Then we solve the intermediate equation U 1 z 1 =p 1 for z 1  by backward substitution. 
   Since the LU solution method for eqs.(7.53,7.55) is stable if A=A 0  or A=A 1 , the BEC algorithm can be numerically unstable in finite-precision arithmetic only if the divisor δ is too small. However, if q=ο(N 1/2 ), then δ=ο(N 1/2 ) in each case.
 
δ 0 ˜(2π) 1/2 e −iπ/8 q 1/4 ,  (7.64) 
 
δ 1 ˜(π/2) 1/2 e −iπ/8 q 1/4 .  (7.65) 
 
Hence, if q=ο(N 2 ), then the BEC algorithm is stable. Each computer solution with the BEC algorithm should be followed by one round of iterative refinement. Solutions of eqs.(7.45,7.49) that are computed this way are accurate to machine precision.
 
8. Efficient evaluation of Chebyshev Interpolation polynomials.
 
   Clenshaw&#39;s algorithm is a stable and accurate method to evaluate the Chebyshev interpolation polynomial 
               u   ⁡     (   ξ   )       =       ∑     k   =   0     N     ⁢       u   k     ⁢       T   k     ⁡     (   ξ   )                   (   8.1   )             
 
at a single argument ξ: u(ξ)=υ 0  where
 
   (1) Start: υ N =u N  and υ N−1 =u N−1 +2ξυ N . 
   (2) For k=N−2, . . . , 1: υ k =u k +2ξυ k+1 −υ k+2.    
   (3) End: υ 0 =u 0 +ξυ 1 −υ 2 . 
   The arithmetic operation count for Clenshaw&#39;s algorithm is about 3N per evaluation. If we use this algorithm to evaluate u(ξ) at the N+1 points
 
ξ n =1−2n/N for n=0, . . . ,N,  (8.2) 
 
then the total operation count for these evaluations is ο(3N 2 ). There are several alternative algorithms that also perform these evaluations but have total operation counts like ο(AN log N), where A is a constant. We shall describe an alternative algorithm that is based on Lagrange polynomial interpolation. If N is moderately large, say N&gt;100, then it is more efficient to evaluate u(ξ) at all N+1 points ξ n  with this algorithm than with Clenshaw&#39;s algorithm.
 
   If υ(θ)=u(cos θ), then N 
               v   ⁡     (   θ   )       =         ∑     l   =   0     N     ⁢       u   k     ⁢     cos   ⁡     (     k   ⁢           ⁢   θ     )       ⁢           ⁢   for   ⁢           ⁢   0       ≤   θ   ≤     π   .               (   8.3   )             
 
Since ξ n =cos φ n  where
 
φ n =2 arcsin{(n/N) 1/2 } for n=0, . . . ,N,  (8.4) 
 
evaluating the polynomial u(ξ) at the regularly spaced points ξ n  is equivalent to evaluating υ(θ) at the irregularly spaced points φ n , which is what the alternative algorithm really does.
 
   The first step of this algorithm is to compute the values of υ(θ) at the regularly spaced points θ j =jπ/3N for j=0, . . . ,3N. A very efficient way to do these computations is to pad the spectral coefficients u n  in eq.(8.3) with 2N trailing zeros and then to take a cosine FFT of length 3N. The total operation count in this approach is ο(15N log N). 
   Next, let K be a positive integer. We shall discuss some choices for K after the algorithm has been described. The second step of the algorithm is to compute the forward differences of orders up to 2K that are associated with the points θ j . Let
 
Δ 0 υ(θ j )=υ(θ j ) for j=0, . . . ,3N.  (8.5) 
 
For k=1, . . . ,2K the forward differences of order k are defined as follows.
 
Δ k υ(θ j )=Δ k−1 υ(θ j+1 )−Δ k−1 υ(θ j ) for j=0, . . . ,3N−k.  (8.6) 
 
The operation count for computing all of these differences is ο(6KN). Note, however, that the subtractions in eq.(8.6) can be done in parallel for each k. Therefore, this step could be expedited on an array processor or on a computer that has a hardware implementation of vector arithmetic.
 
   The third step of the algorithm is based on the Gaussian (forward) form of the Lagrange interpolation polynomials of even order. Let h=π/3N, and let Q 2K (θ j +th) be the following polynomial in t of degree 2K. 
                 Q     2   ⁢           ⁢   K       ⁡     (       θ   j     +   th     )       =       v   ⁡     (     θ   j     )       +     t   ⁢           ⁢   Δ   ⁢           ⁢     v   ⁡     (     θ   j     )         +       1   2     ⁢     t   ⁡     (     t   -   1     )       ⁢     Δ   2     ⁢     v   ⁡     (     θ     j   -   1       )         +       1   6     ⁢     t   ⁡     (       t   2     -   1     )       ⁢     Δ   3     ⁢     v   ⁡     (     θ     j   -   1       )         +       ∑     k   =   2       K   -   1       ⁢       ∏     l   =   1       k   -   1       ⁢       (       t   2     -     l   2       )     ⁡     [           t   ⁡     (     t   -   k     )           (     2   ⁢           ⁢   k     )     !       ⁢     Δ     2   ⁢           ⁢   k       ⁢     v   ⁡     (     θ     j   -   k       )         +         t   ⁡     (       t   2     -     k   2       )           (       2   ⁢           ⁢   k     +   1     )     !       ⁢     Δ       2   ⁢           ⁢   k     +   1       ⁢     v   ⁡     (     θ     j   -   k       )           ]           +       ∏     l   =   1       K   -   1       ⁢       (       t   2     -     l   2       )     ⁢       t   ⁡     (     t   -   K     )           (     2   ⁢           ⁢   K     )     !       ⁢     Δ     2   ⁢           ⁢   K       ⁢       v   ⁡     (     θ     j   -   K       )       .                   (   8.7   )             
 
Then
 
υ(θ j +th)=Q 2K (θ j +th)+R 2K (θ j +th),  (8.8) 
 
where the remainder R 2K (θ j +th) can written in the form 
                 R     2   ⁢           ⁢   K       ⁡     (       θ   j     +   th     )       =     t   ⁢       ∏     l   =   1     K     ⁢       (       t   2     -     l   2       )     ⁢       h       2   ⁢           ⁢   K     +   1           (       2   ⁢           ⁢   K     +   1     )     !       ⁢       ∂       2   ⁢           ⁢   K     +   1         ∂     θ       2   ⁢           ⁢   K     +   1           ⁢     v   ⁡     (   ϑ   )                     (   8.9   )             
 
for some σ depending on t such that θ j−K &lt;σ&lt;θ j+K . The polynomial Q 2K (θ j +th) can be evaluated with about 4K operations as follows. 
             (   1   )             Start   ⁢     :     ⁢           ⁢     q   K       =         t   -   K       2   ⁢           ⁢   K       ⁢     Δ     2   ⁢   K       ⁢       v   ⁡     (     θ     j   -   K       )       .                   (   2   )               For   ⁢           ⁢   k     =     K   -   1       ,   …   ⁢           ,     
     ⁢       1   :     q   k       =           t   -   k       2   ⁢           ⁢   k       ⁡     [         Δ     2   ⁢           ⁢   k       ⁢     v   ⁡     (     θ     j   -   k       )         +         t   +   k         2   ⁢           ⁢   k     +   1       ⁢     (         Δ       2   ⁢           ⁢   k     +   1       ⁢     v   ⁡     (     θ     j   -   k       )         +     q     k   +   1         )         ]       .                 
 
   (3) End: Q 2K (θ j +th)=υ(θ j )+t(Δυ(θ j )+q 1 ). 
   Now for n=0, . . . , N let j n  be the index j for which the absolute difference |φ n −θ j | is smallest, and let t n =(φ n −θ j     n   )/h. Clearly 
              t   n          ≤       1   2     .         
 
Since φ 0 =0=θ 0  and φ n =π=θ 3N , υ(φ 0 )=υ(θ 0 ) and υ(φ n )=υ(θ 3N ). If 0&lt;n&lt;N, then we use eq.(8.8) to approximate υ(φ n ) as follows.
 
υ(φ n )=υ(θ j     n   +t n h)≈Q 2K (θ j     n   +t n h) for n=1, . . . ,N−1.  (8.10) 
 
The total operation count for computing these approximations is about 4KN. Note that in order for thereto be enough sample points θ j , the increment h=π/3N must be small enough that φ 1 /h&gt;K. A sufficient condition for there to be enough of these sample points is that N&gt;(π 2 /36)K 2 .
 
   Let us estimate the interpolation errors
 
R 2K (θ j     n   t n h)=υ(φ n )−Q 2K (θ j     n   +t n h)  (8.11) 
 
in the approximations (8.10). First, we consider the simple cases in which there is just one non-vanishing term in eq.(8.3). Let 0≦k≦N and assume that υ(θ)=cos(kθ). It can be shown that 
                        R     2   ⁢   K       ⁡     (       θ     j   n       +       t   n     ⁢   h       )            ≤       1   2     ⁢       (     π   K     )       1   2       ⁢         (       k   N     ⁢     π   6       )         2   ⁢           ⁢   K     +   1       ⁡     [     1   +     O   ⁡     (     K     -   1       )         ]       ⁢           ⁢   if       ⁢          ⁢       v   ⁡     (   θ   )       =       cos   ⁡     (     k   ⁢           ⁢   θ     )       .               (   8.12   )             
 
For fixed k this bound is an exponentially decreasing function of increasing K because 0≦k/N≦1. For fixed K this bound is a rapidly decreasing function of decreasing k. Thus, the bound on the interpolation error is greatest if k=N. If K=25, however, eq.(8.12) implies double-precision accuracy even in this case because ½(π/25) 1/2 (π/6) 51 =0.8×10 −15 .
 
   In the general case we find that 
                      R     2   ⁢           ⁢   K       ⁡     (       θ     j   n       +       t   n     ⁢   h       )            ≤       1   2     ⁢       (     π   K     )       1   2       ⁢       ∑     k   =   1     N     ⁢            u   k          ⁢           (       k   N     ⁢     π   6       )         2   ⁢           ⁢   K     +   1       ⁡     [     1   +     O   ⁡     (     K     -   1       )         ]       .                   (   8.13   )             
 
In practice, the Chebyshev coefficients u k  eventually decrease exponentially in magnitude as the index k increases. Therefore, eq.(8.13) implies that the approximations in eq.(8.10) may have double-precision accuracy if the integer K is slightly less than 25, e.g., K≈20.
 
   The total operation count for computing all the values u(ξ n )=υ(φ n ) with the alternative algorithm is ο(15N log N+10KN). This should be compared with the total operation count ο( 3 N 2 ) for computing these values directly with Clenshaw&#39;s algorithm. 
   9. The decomposition method for transversely unbounded waveguides. 
   We describe a newly created numerical method for integrating a system of coupled parabolic equations that models one-way propagation of sound in a transversely unbounded waveguide. Referring now to  FIG. 3 , we assume that sound in the waveguide propagates mainly in the direction of increasing range  30 , which is assigned the coordinate x. The waveguide is unbounded in the direction of cross range  32 , which is assigned the coordinate y. This discussion is for slab geometry, but the development for spherical geometry is almost identical. 
   The horizontal variations in the properties of the fluid in this waveguide are confined between vertical planes at y=±W p , which are parallel to the direction of propagation. We assume that the fluid in the waveguide was originally in a horizontally stratified state but has been perturbed between these planes. The unperturbed fluid can have an internal horizontal interface that also has been perturbed between the planes. The perturbations are smooth and vanish at the planes. The waveguide remains horizontally stratified in the two semi-infinite sections where |y|&gt;W p . 
   To make the transverse extent of our problem finite, we put transparent boundaries in the vertical planes at y=±W, where W&gt;W p . These boundaries have no physical reality. Sound in the waveguide must penetrate the transparent boundaries without reflection because they are in the semi-infinite horizontally stratified sections of the waveguide. 
   For simplicity we assume that the reference wavenumbers k m  are positive constants. Thus, the acoustic pressure field is represented by modal amplitudes φ m  for m=1, . . . ,M that satisfy the following system of coupled parabolic equations for x&gt;0. 
                 2   ⁢           ⁢   i   ⁢           ⁢     k   m     ⁢       ∂     φ   m         ∂   x         +         ∂   2     ⁢     φ   m         ∂     y   2           =         (       k   m   2     -     ξ   m   2       )     ⁢     φ   m       +       ∑     n   =   1     M     ⁢       [         (       α     m   ⁢           ⁢   n       +     i   ⁢           ⁢     β     m   ⁢           ⁢   n       ⁢     k   n         )     ⁢     φ   n       +       β     m   ⁢           ⁢   n       ⁢       ∂     φ   n         ∂   x         +       γ     m   ⁢           ⁢   n       ⁢       ∂     φ   n         ∂   y           ]     ⁢     e       -     i   ⁡     (       k   m     -     k   n       )         ⁢   x                     (   9.1   )             
 
for m=1, . . . ,M. These equations are obtained from eq.(2.8) with σ m =k m x for each m.
 
   Each modal wavenumber ξ m (x, y) is constant in the half-planes  34  where |y|&gt;W p , which correspond to the unperturbed horizontally stratified sections of the waveguide. We denote this constant by ξ m   (0) . Note that it is the same for y&lt;−W p  as for y&gt;W p . We expect the reference wavenumber k m  to be close to ξ m   (0) . In addition, the coupling coefficients a mn (x,y),β mn (x,y) and γ mn (x,y) vanish for |y|&gt;W p . 
   Sound enters the waveguide through an aperture  36  in the vertical plane at x=0. The aperture is confined between the lines at y=±W a  in this plane, where W a &lt;W p . The acoustic pressure field in this plane is represented by initial conditions on the modal amplitudes of the following form. 
                 φ   m     ⁡     (     0   ,   y     )       =     {             S   m     ⁡     (   y   )                 if   ⁢           ⁢        y          &lt;     W   s       ,             0           if   ⁢           ⁢        y          &gt;       W   s     .                       (   9.2   )             
 
   Now let us formulate transparent boundary conditions for the modal amplitudes φ m  at the transparent boundaries  38  where y=±W. First, the coupled parabolic equations (9.1) decouple in the half-planes |y|&gt;W p  as follows. For m=1, . . . ,M, 
                 2   ⁢     ik   m     ⁢       ∂     φ   m         ∂   x         +         ∂   2     ⁢     φ   m         ∂     y   2           =         (       k   m   2     -       {     ξ   m     (   0   )       }     2       )     ⁢     φ   m     ⁢           ⁢   if   ⁢           ⁢        y          &gt;       W   p     .               (   9.3   )             
 
We bring each of these equations into the form of a standard parabolic equation by defining the following phase-shifted modal amplitudes.
 
{overscore (φ)} m =φ m e −1x     m     x  for m=1, . . . ,M,  (9.4) 
 
where each κ m  is a constant. If we let
 
κ m =({ξ m   (0) } 2 −k m   2 )/2 k m ,  (9.5) 
 
then {overscore (φ)} m  satisfies the standard parabolic equation 
                 2   ⁢     ik   m     ⁢       ∂       φ   ~     m         ∂   x         +         ∂   2     ⁢       φ   ~     m         ∂     y   2           =       0   ⁢           ⁢   if   ⁢           ⁢        y          &gt;       W   p     .               (   9.6   )             
 
Since W&gt;W p , outgoing transparent boundary conditions can be imposed on the phase-shifted modal amplitudes at the transparent boundaries y=±W. For example, we get the following form of these conditions from eqs.(6.2,6.3). For m=1, . . . ,M, 
                   φ   ~     m     ⁡     (     x   ,     ±   W       )       =         ∓         e     i   ⁢           ⁢     π   /   4         ⁡     (     2   ⁢           ⁢   π   ⁢           ⁢     k   m       )           -   1     /   2         ⁢       ∫   0   x     ⁢       ∂     ∂   y       ⁢         φ   ~     m     ⁡     (     u   ,     ±   W       )       ⁢       (     x   -   u     )         -   1     /   2       ⁢     ⅆ   u     ⁢           ⁢   for   ⁢           ⁢   x         &gt;   0.             (   9.7   )             
 
   Let us reformulate eq.(9.1) in terms of the phase-shifted modal amplitudes. 
                   2   ⁢     ik   m     ⁢       ∂       φ   ~     m         ∂   x         +         ∂   2     ⁢       φ   ~     m         ∂     y   2           =         (         {     ξ   m     (   0   )       }     2     -     ξ   m   2       )     ⁢       φ   ~     m       +       ∑     n   =   1     M     ⁢       [         (       a   mn     +     i   ⁢           ⁢     β   mn     ⁢       k   ~     n         )     ⁢       φ   ~     n       +       β   mn     ⁢       ∂       φ   ~     n         ∂   x         +       γ   mn     ⁢       ∂       φ   ~     n         ∂   y           ]     ⁢     e       -     i   ⁡     (         k   ~     m     -       k   ~     n       )         ⁢   x               ,           (   9.8   )             
 
where the shifted reference wavenumbers {overscore (k)} m  are defined as
 
{overscore (k)} m =k m +κ m =ξ m   (0) +(k m −ξ m   (0) ) 2 /2k m  for m=1, . . . ,M.  (9.9) 
 
Note that {overscore (k)} m  is stationary at k m =ξ m   (0) , where it equals ξ m   (0)  as well.
 
   A sufficient condition that the system of coupled parabolic equations (9.1) have a unique solution for a given set of initial conditions (9.2) is that the reference matrix be positive definite and that the transparent boundary conditions (9.7) be satisfied. 
   Each {overscore (φ)} m  is the sum of a primary component {overscore (φ)} m   (1)  and a secondary component {overscore (φ)} m   (2)  that are the solutions of the following two problems.
 
I. Primary component. Let {overscore (φ)} m   (1)  be a solution of the inhomogeneous parabolic equation 
                 2   ⁢     ik   m     ⁢       ∂       φ   ~     m     (   1   )           ∂   x         +         ∂   2     ⁢       φ   ~     m     (   1   )           ∂     y   2           =         (         {     ξ   m     (   0   )       }     2     -     ξ   m   2       )     ⁢       φ   ~     m       +       ∑     n   =   1     M     ⁢       [         (       α   mn     +     i   ⁢           ⁢     β   mn     ⁢       k   ~     n         )     ⁢       φ   ~     n       +       β   mn     ⁢       ∂       φ   ~     n         ∂   x         +       γ   mn     ⁢       ∂       φ   ~     n         ∂   y           ]     ⁢     ⅇ       -     ⅈ   ⁡     (         k   ~     m     -       k   ~     n       )         ⁢   x                     (   9.10   )             
 
for x&gt;0 and |y|≦W that satisfies the same initial conditions that {overscore (φ)} m , satisfies, 
                   φ   ~     m     (   1   )       ⁡     (     0   ,   y     )       =           φ   ~     m     ⁡     (     0   ,   y     )       =     {             S   m     ⁡     (   y   )                 if   ⁢           ⁢        y          &lt;     W   3       ,             0             if   ⁢           ⁢        y          &gt;     W   3       ,                       (   9.11   )             
 
and that is subject to the following Neumann boundary conditions. 
                   ∂               ∂   y       ⁢         φ   ~     m     (   1   )       ⁡     (     x   ,     ±   W       )         =       0   ⁢           ⁢   for   ⁢           ⁢   x     &gt;   0.             (   9.12   )             
 
Conditions (9.11,9.12) are compatible, and the solution {overscore (φ)} m   (1)  of this problem is unique.
 
II. Secondary component. Let {overscore (k)} m   (2)  be a solution of the standard parabolic equation 
                 2   ⁢     ik   m     ⁢       ∂       φ   ~     m     (   2   )           ∂   x         +         ∂   2     ⁢       φ   ~     m     (   2   )           ∂     y   2           =       0   ⁢           ⁢   for   ⁢           ⁢   x     &gt;     0   ⁢           ⁢   and   ⁢           ⁢        y          ≤   W             (   9.13   )             
 
that satisfies the homogeneous initial condition
 
{overscore (φ)} m   (2) (0,y)=0 for |y|≦W,  (9.14) 
 
and that is subject to the following inhomogeneous Neumann boundary conditions. 
                   ∂               ∂   y       ⁢         φ   ~     m     (   2   )       ⁡     (     x   ,     ±   W       )         =           ∂               ∂   y       ⁢         φ   ~     m     ⁡     (     x   ,     ±   W       )       ⁢           ⁢   for   ⁢           ⁢   x     &gt;   0.             (   9.15   )             
 
Conditions (9.14,9.15) are compatible, and the solution {overscore (φ)} m   (2)  of this problem is unique.
 
   It can be shown that
 
{overscore (φ)} m ={overscore (φ)} m   (1) +{overscore (φ)} m   (2)  for m=1, . . . ,M.  (9.16) 
 
This identity suggests the following basic decomposition method for finding the phase-shifted modal amplitudes {overscore (φ)} m . We regard each {overscore (φ)} m  as an unknown in problems I and II and solve these problems simultaneously for the triples {{overscore (φ)} m ,{overscore (φ)} m   (1) ,{overscore (φ)} m   (2) )} subject to the constraints that {overscore (φ)} m ={overscore (φ)} m   (1) +{overscore (φ)} m   (2)  and that {overscore (φ)} m  satisfy the transparent boundary conditions (9.7) for m=1, . . . ,M. This approach gives us the freedom to use different techniques for solving problems I and II numerically.
 
   In practice, we solve discrete forms of problems I and II numerically. First, we develop a discrete form of problem I from an exact integral of eq.(9.10). Fix a range increment Δx and let x j =jΔx for j=0,1, . . . . Define 
                   ψ   ~     m     (   1   )       =         k   m     ⁢       ψ   ~     m     (   1   )         +       (     i   /   2     )     ⁢       ∑     n   =   1     M     ⁢           ⁢       θ   mn     ⁢       ψ   ~     n     ⁢     ⅇ       -     ⅈ   ⁡     (       k   m     -     k   n       )         ⁢   x                 ,           (   9.17   )             
 
and 
                 q   _     m     =         (         {     ξ   m     (   0   )       }     2     -     ξ   m   2       )     ⁢       φ   ~     m       +       ∑     n   =   1     M     ⁢       [         (       α   mn     +     i   ⁢       k   _     m     ⁢     β   mn       -       ∂     β   mn         ∂   x         )     ⁢       φ   ~     n       +       γ   mn     ⁢       ∂       φ   ~     n         ∂   y         +       i     2   ⁢     k   m         ⁢       ∂   2       ∂     y   2         ⁢     (       β   mn     ⁢       φ   ~     n       )         ]     ⁢       e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢   x       .                   (   9.18   )             
 
We can verify that eq.(9.10) is equivalent to the following equation. 
                 2   ⁢     ik   m     ⁢       ∂       ψ   ~     m     (   1   )           ∂   x         +         ∂   2     ⁢       ψ   ~     m     (   1   )           ∂     y   2           =       k   m     ⁢         q   ~     m     .               (   9.19   )             
 
Note that
 
{overscore (ψ)} m   (1) (x,y)=k m {overscore (φ)} m   (1) (x,y) and {overscore (q)} m (x,y)=0 if |y|&gt;W p   (9.20) 
 
because ξ m =ξ m   (0)  and the coupling coefficients vanish for |y|&gt;W p . Therefore, {overscore (ψ)} m   (1) (x, y) satisfies Neumann boundary conditions like eq.(9.12) at y=±W, and we can use the Fourier cosine transform ℑ to integrate eq.(9.19) from y=−W to y=W. Next, we use the integrating factor exp(iη 2 x/2k m ), where η is the spectral parameter, to integrate the transformed equation from x=x j  to x=x j+1 . Applying ℑ −1  to the result, we get 
                   ψ   ~     m     (   1   )       ⁡     (     x     j   +   1       )       =         ℱ     -   1       ⁡     (     exp   ⁢     {       -   i     ⁢           ⁢     η   2     ⁢   Δ   ⁢           ⁢       x   /   2     /     k   m         }     ⁢   ℱ   ⁢         ψ   ~     m     (   1   )       ⁡     (     x   j     )         )       -       (     i   /   2     )     ⁢       ∫     x   j       x     j   +   1         ⁢         -   1       ⁢     (       -   exp     ⁢     {     i   ⁢           ⁢         η   2     ⁡     (       x     j   +   1       -   x     )       /   2     ⁢     k   m       }     ⁢   ℱ   ⁢           ⁢         q   _     m     ⁡     (   x   )         )     ⁢       ⅆ   x     .                     (   9.21   )             
 
   Now we introduce three numerical approximations in eq.(9.21) that lead to a system of linear equations for the values of {overscore (φ)} m  at discrete mesh points in the x, y plane. First, we use the trapezoidal rule to approximate the integral on the RHS of this equation. Dropping the ο((Δx) 3 ) remainder in this approximation, we get
 
{overscore (ψ)} m   (1) (x j+1 )+(iΔx/4) {overscore (q)} m (x j+1 ) =ℑ −1 (exp{−iη 2 Δx/2k m }ℑ[{overscore (ψ)} m   (1) (x j )−(iΔx/4{overscore (q)} m (x j )]).  (9.22) 
 
This equation is inconsistent with eq.(9.19). However, it is consistent with a useful approximation to eq. (9.19), which is to regard {overscore (ψ)} m   (1)  as a solution of the standard parabolic equation 
                 2   ⁢     ik   m     ⁢       ∂       ψ   ~     m     (   1   )           ∂   x         +         ∂   2     ⁢       ψ   ~     m     (   1   )           ∂     y   2           =       0   ⁢           ⁢   for   ⁢           ⁢        y          ≤   W             (   9.23   )             
 for |y|≦W  (9.23)
 
over the disjoint open intervals x j &lt;x&lt;x j+1  for j=0,1, . . . that is subject to the following limit conditions at the range steps x j .
 
{overscore (ψ)} m   (1) (x j ±0,y)={overscore (ψ)} m   (1) (x j ,y)∓(iΔx/4){overscore (q)} m (x j ,y) for |y|&lt;W.  (9.24) 
 
We can use these conditions to write eq.(9.22) as
 
{overscore (ψ)} m   (1) (x j+1 −0)=ℑ −1 (exp{−iη 2 Δx/2k m }ℑ[{overscore (ψ)} m   (1) (x j +0)]),  (9.25) 
 
and we can combine them to get the following identities. 
                     ψ   ~     m     (   1   )       ⁡     (       x   j     ,   y     )       =       1   2     ⁡     [           ψ   ~     m     (   1   )       ⁡     (         x   j     +   0     ,   y     )       +         ψ   ~     m     (   1   )       ⁡     (         x   j     -   0     ,   y     )         ]         ,           (   9.26   )                       ψ   ~     m     (   1   )       ⁡     (         x   j     +   0     ,   y     )       -         ψ   ~     m     (   1   )       ⁡     (         x   j     -   0     ,   y     )         =       (       -   i     ⁢           ⁢   Δ   ⁢           ⁢     x   /   2       )     ⁢           q   _     m     ⁡     (       x   j     ,   y     )       .               (   9.27   )             
 
   Next, we fix a positive integer N, let Δy=2W/N and define the points y k =kΔy−W for k=0,1, . . . ,N. Our second approximation is to replace the operators ℑ and ℑ −1  in eq. (9.22) by the operators ℑ N  and ℑ N   1 . Since these operators transform sequences to sequences, we must replace the functions {overscore (ψ)} m   (1) ±(iΔx/4){overscore (q)} m  in this equation with the sequences of their values at the points y k . Let us write these sequences as {{overscore (ψ)} m   (1) (y k )±(iΔx/4){overscore (q)} m (y k )} k=0   N . Under these approximations, eq.(9.22) becomes
 
{{overscore (ψ)} m   (1) (x j+1 ,y k )+(iΔx/4){overscore (q)} m (x j+1 ,y k )} k=0   N =ℑ −1 (exp{−iη 2 Δx/2k m }ℑ N {overscore (ψ)} m   (1) (x j ,y l )−(iΔx/4){overscore (q)} m (x j ,y l )} l=0   N   (9.28) 
 
The spectral parameter η in this equation is bounded because η=nπ/2W for n=0,1, . . . ,N. Hence, 0≦η≦π/Δy.
 
   We must approximate the partial derivatives σ{overscore (φ)} n /σy and σ 2 (β mn {overscore (φ)} n )/σ 2 y that appear in eq. (9.18) for {overscore (φ)} m  before we can evaluate the sequences in eq.(9.28). Since the variable y in this equation is restricted to the uniformly spaced points y k , our third approximation is to replace the partial derivatives at these points with central differences. At the lowest order of approximation, we use the three-point formulas 
                   ∂     ∂   y       ⁢         φ   ~     n     ⁡     (     y   k     )         =         1     2   ⁢           ⁢   Δ   ⁢           ⁢   y       ⁡     [           φ   ~     n     ⁡     (     y     k   +   1       )       -         φ   ~     n     ⁡     (     y     k   -   1       )         ]       +     O   ⁡     (       (     Δ   ⁢           ⁢   y     )     2     )           ,           (   9.29   )                     ∂   2       ∂     y   2         ⁢     (       β   mn     ⁢       φ   ~     n       )     ⁢     (     y   k     )       =         1       (     Δ   ⁢           ⁢   y     )     2       ⁡     [         (       β   mn     ⁢       φ   ~     n       )     ⁢     (     y     k   +   1       )       -     2   ⁢     (       β   mn     ⁢       φ   ~     n       )     ⁢     (     y   k     )       +       (       β   mn     ⁢       φ   ~     n       )     ⁢     (     y     k   -   1       )         ]       +       O   ⁡     (       (     Δ   ⁢           ⁢   y     )     2     )       .               (   9.30   )             
 
   Now we combine these approximations with eq.(9.18). Rearranging terms in the result, we find that 
                   q   _     m     ⁡     (     y   k     )       =         (         {     ξ   m     (   0   )       }     2     -       ξ   m   2     ⁡     (     y   k     )         )     ⁢         φ   ~     m     ⁡     (     y   k     )         +       ∑     n   =   1     M     ⁢       [         α   mn     ⁡     (     y   k     )       +       (       i   ⁢       k   _     m       -     i         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         )     ⁢       β   mn     ⁡     (     y   k     )         -       ∂     ∂   x       ⁢       β   mn     ⁡     (     y   k     )           ]     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢   x       ⁢         φ   ~     n     ⁡     (     y   k     )           +       ∑     n   =   1     M     ⁢       (         1     2   ⁢   Δ   ⁢           ⁢   y       ⁢       γ   mn     ⁡     (     y   k     )         +       i     2   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (     y     k   +   1       )           )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢   x       ⁢         φ   ~     n     ⁡     (     y     k   +   1       )           +       ∑     n   =   1     M     ⁢       (         -     1     2   ⁢   Δ   ⁢           ⁢   y         ⁢       γ   mn     ⁡     (     y   k     )         +       i     2   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (     y     k   -   1       )           )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢   x       ⁢         φ   ~     n     ⁡     (     y     k   -   1       )           +       O   ⁡     (       (     Δ   ⁢           ⁢   y     )     2     )       .               (   9.31   )             
 
Therefore, neglecting terms of ο(Δx(Δy) 2 ), we find that 
                     ψ   ~     m     (   1   )       ⁡     (       x     j   +   1       ,     y   k       )       +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢         q   _     m     ⁡     (       x     j   +   1       ,     y   k       )           =           k   _     m     ⁢         φ   ~     m     (   1   )       ⁡     (       x     j   +   1       ,     y   k       )         +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢     (         {     ξ   m     (   0   )       }     2     -       ξ   m   2     ⁡     (       x     j   +   1       ,     y   k       )         )     ⁢         φ   ~     m     ⁡     (       x     j   +   1       ,     y   k       )         +       (       i   2     -           k   _     m     ⁢   Δ   ⁢           ⁢   x     4     +       Δ   ⁢           ⁢   x       4   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2           )     ⁢       ∑     n   =   1     M     ⁢         β   mn     ⁡     (       x     j   +   1       ,     y   k       )       ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x     j   +   1           ⁢         φ   ~     n     ⁡     (       x     j   +   1       ,     y   k       )             +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢       ∑     n   =   1     M     ⁢       (         α   mn     ⁡     (       x     j   +   1       ,     y   k       )       -       ∂     ∂   x       ⁢       β   mn     ⁡     (       x     j   +   1       ,     y   k       )           )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x     j   +   1           ⁢         φ   ~     n     ⁡     (       x     j   +   1       ,     y   k       )             +       ∑     n   =   1     M     ⁢       (           i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢           ⁢   Δ   ⁢           ⁢   y       ⁢       γ   mn     ⁡     (       x     j   +   1       ,     y   k       )         -         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (       x     j   +   1       ,     y     k   +   1         )           )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x     j   +   1           ⁢         φ   ~     n     ⁡     (       x     j   +   1       ,     y     k   +   1         )           +       ∑     n   =   1     M     ⁢       (         -       i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢           ⁢   Δ   ⁢           ⁢   y         ⁢       γ   mn     ⁡     (       x     j   +   1       ,     y   k       )         -         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (       x     j   +   1       ,     y     k   -   1         )           )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x     j   +   1           ⁢           φ   ~     n     ⁡     (       x     j   +   1       ,     y     k   -   1         )       .                   (   9.32   )             
 
   Next, we use the constraint that {overscore (φ)} m ={overscore (φ)} m   (1) +{overscore (φ)} m   (2)  to replace the term k m {overscore (φ)} m   (1) (x +1 ,y k ) on the RHS of eq.(9.32) with k m {overscore (φ)} m (x +1 ,y k ) by adding k m {overscore (φ)} m   (2) (x +1 ,y k ) to each side of this equation. Thus,
 
{overscore (ψ)} m   (1) (x j+1 ,y k )+(iΔx/4){overscore (q)} m (x j+1 ,y k )={overscore (ψ)} m (x j+1 −0,y k )−k m {overscore (ψ)} m   (2) (x j+1 ,y k )  (9.33) 
 
where 
                   Ψ   ~     m     ⁡     (         x     j   +   1       -   0     ,     y   k       )       =         [       k   m     +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢     (         {     ξ   m     (   0   )       }     2     -       ξ   m   2     ⁡     (       x     j   +   1       ,     y   k       )         )         ]     ⁢         φ   ~     m     ⁡     (       x     j   +   1       ,     y   k       )         +       (       i   2     -           k   _     m     ⁢   Δ   ⁢           ⁢   x     4     +       Δ   ⁢           ⁢   x       4   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2           )     ⁢       ∑     n   =   1     M     ⁢       β   mn     ⁢     (       x     j   +   1       ,     y   k       )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x     j   +   1           ⁢         φ   ~     n     ⁡     (       x     j   +   1       ,     y   k       )             +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢       ∑     n   =   1     M     ⁢       (         α   mn     ⁡     (       x     j   +   1       ,     y   k       )       -       ∂     ∂   x       ⁢     β   mn     ⁢     (       x     j   +   1       ,     y   k       )         )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x     j   +   1           ⁢         φ   ~     n     ⁡     (       x     j   +   1       ,     y   k       )             +       ∑     n   =   1     M     ⁢       (           i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢           ⁢   Δ   ⁢           ⁢   y       ⁢       γ   mn     ⁡     (       x     j   +   1       ,     y   k       )         -         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (       x     j   +   1       ,     y     k   +   1         )           )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x     j   +   1           ⁢         φ   ~     n     ⁡     (       x     j   +   1       ,     y     k   +   1         )           +       ∑     n   =   1     M     ⁢       (         -       i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢           ⁢   Δ   ⁢           ⁢   y         ⁢       γ   mn     ⁡     (       x     j   +   1       ,     y   k       )         -         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (       x     j   +   1       ,     y     k   -   1         )           )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x     j   +   1           ⁢           φ   ~     n     ⁡     (       x     j   +   1       ,     y     k   -   1         )       .                   (   9.34   )             
 
Similarly, and to the same order of approximation,
 
{overscore (ψ)} m   (1) (x j ,y l )−(iΔx/4){overscore (q)} m (x j ,y l )={overscore (ψ)} m (x j +0,y l )−k m {overscore (ψ)} m   (2) (x j ,y l )  (9.35) 
 
where 
                   Ψ   ~     m     ⁡     (         x   j     -   0     ,     y   l       )       =         [       k   m     -       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢     (         {     ξ   m     (   0   )       }     2     -       ξ   m   2     ⁡     (       x   j     ,     y   l       )         )         ]     ⁢         φ   ~     m     ⁡     (       x   j     ,     y   l       )         +       (       i   2     +           k   _     m     ⁢   Δ   ⁢           ⁢   x     4     -       Δ   ⁢           ⁢   x       4   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2           )     ⁢       ∑     n   =   1     M     ⁢       β   mn     ⁢     (       x   j     ,     y   l       )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x   j         ⁢         φ   ~     n     ⁡     (       x   j     ,     y   l       )             +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢       ∑     n   =   1     M     ⁢       (         α   mn     ⁡     (       x   j     ,     y   l       )       -       ∂     ∂   x       ⁢       β   mn     ⁡     (     x   ,     y   l       )           )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x     j   +   1           ⁢         φ   ~     n     ⁡     (       x   j     ,     y   l       )             +       ∑     n   =   1     M     ⁢       (         -       i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢           ⁢   Δ   ⁢           ⁢   y         ⁢       γ   mn     ⁡     (       x   j     ,     y   l       )         +         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (       x   j     ,     y     l   +   1         )           )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x   j         ⁢         φ   ~     n     ⁡     (       x   j     ,     y     l   +   1         )           +       ∑     n   =   1     M     ⁢       (           i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢           ⁢   Δ   ⁢           ⁢   y       ⁢       γ   mn     ⁡     (       x   j     ,     y   l       )         -         Δ   ⁢           ⁢   x       8   ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β   mn     ⁡     (       x   j     ,     y     l   -   1         )           )     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x   j         ⁢           φ   ~     n     ⁡     (       x   j     ,     y     l   -   1         )       .                   (   9.36   )             
 
   We shall assume that Δy&lt;W−W p . Under this assumption the coupling coefficients a mn (x,y k ), β mn (x,y k ) and γ mn (x,y k ) vanish for x&gt;0 and k=1, N−1. This eliminates the boundary values {overscore (σ)} m (x j+1 ,±W) from the RHS of eq.(9.34), which makes it possible for us to find the internal values {overscore (σ)} m (x j+1 ,y k ) (k=1, . . . ,N−1) from this equation. Let us combine eqs.(9.28,9.33,9.35) as follows.
 
{{overscore (Ψ)} m (x j+1 −0,y k )−k m {overscore (ψ)} m   (2) (x j+1 ,y k )} k=0   N 
 
=ℑ N   −1 (exp{−iη 2 Δx/2k m }ℑ N {{overscore (ψ)} m (x j +0,y l )−k m {overscore (ψ)} m   (2) (x j ,y l )} l=0   N ).  (9.37)
 
If the values {overscore (ψ)} m   (2) (x j+1 ,y k )(k=1, . . . ,N−1) at the range step x j+1  are known, then eqs.(9.34,9.37) determine a system of linear equations for the values {overscore (ψ)} m (x j+1 ,y k )(k=1, . . . ,N−1) at the range step x j+1 . Let us put these equations in matrix form. First, define the column vectors
 
{overscore (ψ)} k =[{overscore (ψ)} 1 (x j+1 ,y k ), . . . ,{overscore (ψ)} M (x j+1 ,y k )] T   (9.38) 
 
{overscore (Ψ)} k   − =[{overscore (Ψ)} 1 (x j+1 −0,y k ), . . . ,{overscore (Ψ)} M (x j+1 −0,y k ] T ,  (9.39) 
 
for k=1, . . . ,N−1. Next, define the M×M matrices A kk  and A kk+1  as follows. 
                 (     A     k   ⁢           ⁢   k       )       m   ⁢           ⁢   n       =     {               k   n     +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢     (         {     ξ   n     (   0   )       }     2     -       ξ   n   2     ⁡     (       x     j   +   1       ,     y   k       )       +       α     n   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )         )         ,           (       if   ⁢           ⁢   m     =   n     )               [         (       i   2     -           k   _     m     ⁢   Δ   ⁢           ⁢   x     4     +       Δ   ⁢           ⁢   x       4   ⁢           ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2           )     ⁢       β     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )         +                                     (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢     (         α     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )       -       ∂     ∂   x       ⁢       β     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )           )       ]     ⁢     e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x     j   +   1             ,           (       if   ⁢           ⁢   m     ≠   n     )                     (   9.40   )                   (     A       k   ⁢           ⁢   k     ±   1       )       m   ⁢           ⁢   n       =       (         ±       i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢   Δ   ⁢           ⁢   y         ⁢       γ     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )         -         Δ   ⁢           ⁢   x       8   ⁢           ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y     k   ±   1         )           )     ⁢       e       -     i   ⁡     (         k   _     m     -       k   _     n       )         ⁢     x     j   +   1           .               (   9.41   )             
 
   The assumption that Δy&lt;W−W p  eliminates the boundary values {overscore (ψ)} m (x j+1 , ±W) from the RHS of eq.(9.34). Thus, it follows from this equation that 
                 [           A   11           A   12                                                   A   21           A   22           A   23                                                 ⋰       ⋰       ⋰                                                 A     N   -     2   ⁢   N     -   3             A     N   -     2   ⁢   N     -   2             A     N   -     2   ⁢   N     -   1                                                     A     N   -     1   ⁢   N     -   2             A     N   -     1   ⁢   N     -   1             ]     ·     [             φ   ~     1                 φ   ~     2             ⋮               φ   ~       N   -   2                   φ   ~       N   -   1             ]       =       [             Ψ   ~       1   _                   Ψ   ~       2   _               ⋮               Ψ   ~         N   _     -   2                   Ψ   ~         N   _     -   1             ]     .             (   9.42   )             
 
If we use eq.(9.37) and the values {overscore (ψ)} m   (2) (x j+1 ,y k )(k=1, . . . ,N−1) to evaluate the vectors {overscore (Ψ)} k   −  on the RHS of this equation, then we can solve it for the internal solution vectors {overscore (ψ)} k .
 
   We can solve eq.(9.42) numerically by Gaussian elimination with partial pivoting (GEPP) followed by one round of iterative refinement. The block tri-diagonal system matrix on the LHS of eq.(9.42) is also a band matrix of order M(N−1) that has 2M−1 lower diagonals and 2M−1 upper diagonals. Since M&lt;&lt;N in practice, it is most efficient to do the GEPP computations with subroutines that take advantage of band storage. 
   The system matrix on the LHS of eq. (9.42) is nonsingular if it is strictly diagonally dominant by columns. We see from eqs.(9.40,9.41) that this condition requires the following strict inequalities to be satisfied for the block indices k=1, . . . ,N−1: 
     |k n +(iΔx/4)({ξ n   (0) } 2 −ξ n   2 (x j+1 y k )+a nn (x j+1 ,y k ))|                            k   n     +       (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢     (         {     ξ   n     (   0   )       }     2     -       ξ   n   2     ⁡     (       x     j   +   1       ,     y   k       )       +       α     n   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )         )              &gt;                 ∑       m   =   1       (     m   ≠   n     )       M     ⁢              (       i   2     -           k   _     m     ⁢   Δ   ⁢           ⁢   x     4     +       Δ   ⁢           ⁢   x       4   ⁢           ⁢         k   m     ⁡     (   Δy   )       2           )     ⁢       β     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )         +                         (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4       )     ⁢     (         α     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )       -       ∂     ∂   x       ⁢       β     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )           )            +                   ∑     m   =   1     M     ⁢                i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢           ⁢   Δ   ⁢           ⁢   y       ⁢       γ     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y     k   -   1         )         -         Δ   ⁢           ⁢   x       8   ⁢           ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )                  +                 ∑     m   =   1     M     ⁢              -       i   ⁢           ⁢   Δ   ⁢           ⁢   x       8   ⁢           ⁢   Δ   ⁢           ⁢   y         ⁢       γ     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y     k   +   1         )         -         Δ   ⁢           ⁢   x       8   ⁢           ⁢         k   m     ⁡     (     Δ   ⁢           ⁢   y     )       2         ⁢       β     m   ⁢           ⁢   n       ⁡     (       x     j   +   1       ,     y   k       )                                for   ⁢           ⁢   n     =   1     ,   …   ⁢           ,     M   .                   (   9.43   )             
 for n= 1 , . . . ,M.  ( 9.43 ) 
   If this condition is satisfied, then Gaussian elimination requires no pivoting and is stable. 
   Note the factors {overscore (k)} m Δx, Δx/Δy and Δx/k m (Δy) 2  that appear in this equation. We recommend choosing the increments Δx and Δy so that k m Δx≦ο(1) and k m Δy≧ο(1). These restrictions imply that the ratios Δx/Δy and Δx/k m (Δy) 2  are at most ο(1). 
   If Δy&lt;W−W p , then we can find the boundary values {overscore (ψ)} m   (1) (x j+1 ,±W) of the primary component from eq.(9.37) without the values {overscore (ψ)} m   (2) (x j+1 ,y k ). In this case, it follows from eq.(9.32) that {overscore (ψ)} m   (1) (x j+1 ,±W)+(iΔx/4){overscore (q)} m (x j+1 ,±W)=k m {overscore (ψ)} m   (1) (x j+1 ,±W). Thus, eqs.(9.33,9.37) imply that 
                 k   m     ⁢         ψ   ~     m     (   1   )       ⁡     (       x     j   +   1       ,     ±   W       )         =         N     -   1         (     vh   ±   w     )       ⁢       (     exp   ⁢     {       -   i     ⁢           ⁢     η   2     ⁢   Δ   ⁢           ⁢     x   /   2     ⁢     k   m       }     ⁢     N     ⁢       {           Ψ   ~     m     ⁡     (         x   j     +   0     ,   yl     )       -       k   m     ⁢         ψ   ~     m     (   2   )       ⁡     (       x   j     ,   yn     )           }       t   =   0     N       )     .               (   9.44   )             
 
Note that it is not necessary to compute the whole inverse DFT to evaluate the RHS of this equation because we need only the values from the inverse DFT at the endpoints y k =±W.
 
   This completes our development of the discrete form of problem I and the method for solving it. Next, we develop a discrete form of problem II. We begin by integrating eq.(9.13) from x=x j  to x=x j+1  with the trapezoidal rule. Dropping the ο((Δx) 3 ) remainder in this approximation, we get the semi-discrete parabolic equation 
                     φ   ~     m     (   2   )       ⁡     (       x     j   +   1       ,   y     )       -         φ   m     (   2   )       ~     ⁡     (       x   j     ,   y     )         =           (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4     ⁢           ⁢     k   m       )     ⁡     [           ∂   2       ∂     y   2         ⁢         φ   ~     m     (   2   )       ⁡     (       x   j     ,   y     )         +         ∂   2       ∂     y   2         ⁢         φ   ~     m     (   2   )       ⁡     (       x     j   +   1       ,   y     )           ]       ⁢           ⁢   for   ⁢           ⁢        y          ≤     W   .               (   9.45   )             
 
We adopt conditions (9.14,9.15) without change. In particular, we require that the inhomogeneous Neumann boundary conditions (9.15) be satisfied at each range step. 
                 ∂     ∂   y       ⁢         φ   ~     m     (   2   )       ⁡     (       x     j   +   1       ,     ±   W       )         =       ∂     ∂   y       ⁢           φ   ~     m     ⁡     (       x     j   +   1       ,     ±   W       )       .               (   9.46   )             
 
   In order to formulate transparent boundary conditions for the discretized amplitude {overscore (ψ)} m , we assume that it satisfies the semi-discrete parabolic equation that corresponds to eq.(9.6). 
                     φ   ~     m     ⁡     (       x     j   +   1       ,   y     )       -         φ   ~     m     ⁡     (       x   j     ,   y     )         =           (     i   ⁢           ⁢   Δ   ⁢           ⁢     x   /   4     ⁢           ⁢     k   m       )     ⁡     [         ∂     ∂     y   2         ⁢         φ   ~     m     ⁡     (       x   j     ,   y     )         +         ∂   2       ∂     y   2         ⁢         φ   ~     m     ⁡     (       x     j   +   1       ,   y     )           ]       ⁢           ⁢   for   ⁢           ⁢        y          &gt;       W   p     .               (   9.47   )             
 
Under this assumption, the discretized amplitude {overscore (ψ)} m  satisfies discrete transparent boundary conditions like eq.(6.31). We can adapt the algorithm described by eqs.(6.26-6.30) to formulate these conditions as follows. 
                       φ   ~     m     ⁡     (       x     j   +   1       ,     ±   W       )       ±       1     p   m       ⁢     ∂     ∂   y       ⁢         φ   ~     m     ⁡     (       x     j   +   1       ,     ±   W       )           =       2   ⁢           ⁢     p   m     ⁢       ∑     k   =   0     j     ⁢         (     -   1     )       j   -   k       ⁢       a       k   +   1     ,   m     k     ⁡     (     ±   W     )       ⁢           ⁢   for   ⁢           ⁢   j         =   0       ,     1   ⁢   …     ⁢           ,           (   9.48   )             
 
where
 
p m =e −iπ/4 2(k m /Δx) 1/2   (9.49) 
 
and the coefficients a k+1,m   k (±W) are computed recursively. First, a 1,m   0 (±W)=0. Second, for j=1,2, . . . , 
                   a     k   ,   m     j     ⁡     (     ±   W     )       =         -     1   2       ⁢       a     k   ,   m       j   -   1       ⁡     (     ±   W     )         -       1   2     ⁢       a       k   -   2     ,   m       j   -   1       ⁡     (     ±   W     )             ,     for   ⁢           ⁢       (       k   =   1     ,   …   ⁢           ,     j   -   1       )       ︸       Skip   ⁢           ⁢   this   ⁢           ⁢   when   ⁢           ⁢   j     =   1.           ,           (   9.50   )                     a     j   ,   m     j     ⁡     (     ±   W     )       =         ±     1   2       ⁢       p   m     -   2       ⁡     [         ∂     ∂   y       ⁢         φ   ~     m     ⁡     (       x     j   -   1       ,     ±   W       )         +       ∂     ∂   y       ⁢         φ   ~     m     ⁡     (       x   j     ,     ±   W       )           ]         -       a     j   ,   m       j   -   1       ⁡     (     ±   W     )       -       1   2     ⁢       a       j   -   2     ,   m       j   -   1       ⁡     (     ±   W     )             ,           (   9.51   )                     a       j   +   1     ,   m     j     ⁡     (     ±   W     )       =         ∓     1   2       ⁢       p   m     -   2       ⁡     [         ∂     ∂   y       ⁢         φ   ~     m     ⁡     (       x     j   -   1       ,     ±   W       )         +       ∂     ∂   y       ⁢         φ   ~     m     ⁡     (       x   j     ,     ±   W       )           ]         -       1   2     ⁢       a       j   -   1     ,   m       j   -   1       ⁡     (     ±   W     )             ,           (   9.52   )             
 
where
 
a −1,m   j (±W)=0, a 0,m   j (±W)=0 for j=0,1, . . .  (9.53) 
 
   Now we combine the constraints {overscore (ψ)} m (x j+1 ,±W)={overscore (ψ)} m   (1) (x j+1 ,±W)+{overscore (ψ)} m   (2) (x j+1 ,±W) and eq.(9.49) with eq.(9.48) to develop new boundary conditions for {overscore (ψ)} m   (2) . Note that the boundary values of the primary component {overscore (ψ)} m   (1)  appear as source terms in these conditions. 
                       φ   ~     m     (   2   )       ⁡     (       x     j   +   1       ,     ±   W       )       ±       1     p   m       ⁢     ∂     ∂   y       ⁢         φ   ~     m     (   2   )       ⁡     (       x     j   +   1       ,     ±   W       )           =         -         φ   ~     m     (   1   )       ⁡     (       x     j   +   1       ,     ±   W       )         +     2   ⁢           ⁢     p   m     ⁢       ∑     k   =   0     j     ⁢         (     -   1     )       j   -   k       ⁢       a       k   +   1     ,   m     k     ⁡     (     ±   W     )       ⁢           ⁢   for   ⁢           ⁢   j           =   0       ,   1   ,     …   ⁢           .             (   9.54   )             
 
We also combine condition (9.46) with recursion formulas (9.51,9.52) to reformulate them in terms of {overscore ( 104  )} m   (2) . 
                   a     j   ,   m     j     ⁡     (     ±   W     )       =         ±     1   2       ⁢       p   m     -   2       ⁡     [         ∂     ∂   y       ⁢         φ   ~     m     (   2   )       ⁡     (       x     j   -   1       ,     ±   W       )         +       ∂     ∂   y       ⁢         φ   ~     m     (   2   )       ⁡     (       x   j     ,     ±   W       )           ]         -       a     j   ,   m       j   -   1       ⁡     (     ±   W     )       -       1   2     ⁢       a       j   -   2     ,   m       j   -   1       ⁡     (     ±   W     )             ,           (   9.55   )                   a       j   +   1     ,   m     j     ⁡     (     ±   W     )       =         ∓     1   2       ⁢       p   m     -   2       ⁡     [         ∂     ∂   y       ⁢         φ   ^     m     (   2   )       ⁡     (       x     j   -   1       ,     ±   W       )         +       ∂     ∂   y       ⁢         φ   ~     m     (   2   )       ⁡     (       x   j     ,     ±   W       )           ]         -       1   2     ⁢         a       j   -   1     ,   m       j   -   1       ⁡     (     ±   W     )       .                 (   9.56   )             
 
Therefore, in principle the discrete form of problem II consists of
         (1) the semi-discrete parabolic equation (9.45), subject to   (2) the initial condition (9.14), and   (3) the boundary conditions (9.54), with   (4) the recursion formulas (9.50,9.53,9.55,9.56).
 
However, in practice we would use the fixed-length truncated algorithm described by eqs.(6.72-6.76) to evaluate the boundary conditions for this problem. Let L J  denote the length of the algorithm, which could be chosen according to criterion. (6.82) for calculations in double-precision arithmetic. First, we use eq.(6.78) to modify the boundary conditions (9.54) as follows. 
                       φ   ~     m     (   2   )       ⁡     (       x     j   +   1       ,     ±   W       )       ±       1     p   m       ⁢     ∂     ∂   y       ⁢         φ   ~     m     (   2   )       ⁡     (       x     j   +   1       ,     ±   W       )           =         -         φ   ~     m     (   1   )       ⁡     (       x     j   +   1       ,     ±   W       )         +     2   ⁢           ⁢     p   m     ⁢       ∑     k   =   0     j     ⁢         (     -   1     )       j   -   k       ⁢       α       L   J     ,   m     k     ⁡     (     ±   W     )       ⁢           ⁢   for   ⁢           ⁢   j           =   0       ,   1   ,     …   ⁢           .             (   9.57   )             
 
Second, we use eqs.(6.72-6.76) to replace eqs.(9.50,9.53,9.55,9.56) with the following recursion equations for the coefficients a L     J     m   k (±W). 
                   α     l   ,   m     j     ⁡     (     ±   W     )       =         -     1   2       ⁢       α       l   +   1     ,   m       j   -   1       ⁡     (     ±   W     )         -       1   2     ⁢       α       l   -   1     ,   m       j   -   1       ⁡     (     ±   W     )             ⁢     
     ⁢         for   ⁢           ⁢   l     =   1     ,   …   ⁢           ,       L   J     -   2     ,             (   9.58   )                     α         L   J     -   1     ,   m     j     ⁡     (     ±   W     )       =         ±     1   2       ⁢       p   m     -   2       ⁡     [         ∂     ∂   y       ⁢         φ   ~     m     (   2   )       ⁡     (       x     j   -   1       ,     ±   W       )         +       ∂     ∂   y       ⁢         φ   ~     m     (   2   )       ⁡     (       x   j     ,     ±   W       )           ]         -       α       L   J     ,   m       -   1       ⁡     (     ±   W     )       -       1   2     ⁢       α         L   J     -   2     ,   m       j   -   1       ⁡     (     ±   W     )             ,           (   9.59   )                     α       L   J     ,   m     j     ⁡     (     ±   W     )       =         ∓     1   2       ⁢       p   m     -   2       ⁡     [         ∂     ∂   y       ⁢         φ   ~     m     (   2   )       ⁡     (       x     j   -   1       ,     ±   W       )         +       ∂     ∂   y       ⁢         φ   ~     m     (   2   )       ⁡     (       x   j     ,     ±   W       )           ]         -       1   2     ⁢       α         L   J     -   1     ,   m       j   -   1       ⁡     (     ±   W     )             ,           (   9.60   )             
 
where
 
a l,m   0 (±W)=0 for l=1, . . . ,L J , a 0,m   j (±W)=0 for j=0,1, . . .  (9.61) 
       

   This completes the development of the discrete form of problem II. A sufficient condition for it to be stable is that 
                   Δ   ⁢           ⁢   x           k   0     ⁡     (     Δ   ⁢           ⁢   y     )       2       &lt;     2   π       ,           (   9.62   )             
 
where k 0 =min k m (m=1, . . . ,M). Numerical integration of eq.(9.45) can introduce numerical artifacts in the solutions of the discrete form of problem II. However, these artifacts will be suppressed if the spatial frequency limit of the integration method does not exceed the following bound.
 
λ c =3 −1/4 2(k 0 /Δx) 1/2 .  (9.63) 
 
This criterion imposes a simple constraint on the N-th order Chebyshev tau integration method. Since the spatial frequency limit of the N-th order Chebyshev tau integration method is 2/Δy, this limit is less than λ c  if 
                 Δ   ⁢           ⁢   x           k   0     ⁡     (     Δ   ⁢           ⁢   y     )       2       &lt;       3       -   1     /   2       .             (   9.64   )             
 
This condition is only slightly more restrictive than condition (9.62).
 
   Referring now to  FIG. 4 , we outline the operation of the preferred embodiment. We begin the process in block  40  by computing the modal wavenumbers and the coupling coefficients at the mesh points (x j ,y k )(k=0,1, . . . ,N) from environmental data. They can be computed for all j=0,1, . . . and stored for subsequent use or computed as needed, beginning at j=0. In block  42  we select the initial values {overscore (ψ)} m (x 0 ,y k )(j=0 and k=0,1, . . . ,N) for the source field in the waveguide at the plane x=x 0 . Now we go to block  44 . This block is the start of the general cycle of steps for solving problems I and II simultaneously. We must compute the values {overscore (ψ)} m (x j+1 ,y k ) and {overscore (ψ)} m   (2) (x j+1 ,y k ) for m=1, . . . ,M and k=0,1, . . . ,N at an arbitrary range step x j +1 from the values {overscore (ψ)} m (x j ,y k ) and {overscore (ψ)} m   (2) (x j ,y l ) for m=1, . . . ,M and l=0,1, . . . ,N at the preceding range step x j . 
   (1) In block  46 , for each m we compute the boundary values {overscore (ψ)} m   (1) (x j+1 ,±W) of the primary component from eq.(9.44). These values are used in the next step. 
   (2) In block  48 , for each m we integrate the semi-discrete parabolic equation (9.45) subject to the boundary conditions (9.57), which we evaluate with eqs.(9.58-9.61). We use the Chebyshev tau method to do the integration numerically. This method transforms the integration to an algebraic problem for the Chebyshev coefficients of {overscore (ψ)} m   (2) (x j+1 ,y) that we solve with the BEC algorithm. It is then necessary to compute the values {overscore (ψ)} m   (2) (x j+1 ,y k ) for k=0,1, . . . ,N from these coefficients. Depending on the size of N, we use either Clenshaw&#39;s algorithm or an alternative algorithm to compute them. We use the internal values {overscore (ψ)} m   (2) (x j+1 ,y k )(k=1, . . . ,N−1) in the next step. We compute the differential boundary values σ{overscore (ψ)} m   (1) (x j+1 ,±W)/σy and save them for use by eqs.(9.59,9.60) at the following two range steps, x j+2  and x j+3 . Having found {overscore (ψ)} m   (1) (x j+1 ,±W) in the previous step and {overscore (ψ)} m   (2) (x j+1 ,±W) in the present step, we compute the boundary values {overscore (ψ)} m (x j+1 ,±W) from the constraints that {overscore (ψ)} m (x j+1 ,±W)={overscore (ψ)} m   (1) (x j+1 ,±W)+{overscore (ψ)} m   (2) (x j+1 ,±W). 
   (3) In block  50 , we solve eq.(9.42) numerically by GEPP for the internal values {overscore (ψ)} m (x j+1 ,y k ) (k=1, . . . ,N−1). This completes the general cycle of steps for solving problems I and II. 
   In block  52  we use the modal amplitudes to evaluate the eigenfunction expansion of the complex pressure at desired depths in the waveguide at range x j+1  and cross ranges y k (k=0,1, . . . ,N). We predict the physical acoustic pressure from these computations. Finally, we increase j to j+1 in block  54  and return to block  44  for the next cycle of computations. This completes the outline. 
   The discretized solutions {overscore (ψ)} m  turn out to be inaccurate at very large range from the source due to small inaccuracies in the discretized solutions {overscore (ψ)} m   (2) . However, an extrapolation method can be applied to the discretized solutions {overscore (ψ)} m   (2)  to compensate for the inaccuracies in them under conditions of practical interest. Before describing this method, let us modify the notation for these solutions to indicate their dependence on the increment Δx.
 
{overscore (ψ)} m   (2) (x j ,y)→{overscore (ψ)} m   (2) (x j ,y i Δx)  (9.65) 
 
A simple way to improve the accuracy of these approximations is to reduce the step size Δx. Since they are even functions of Δx, it is efficient to combine this approach with Richardson extrapolation.
 
   We shall discuss a procedure that uses Richardson extrapolation to estimate {overscore (ψ)} m   (2) (x j ,y) from the computed values of {overscore (ψ)} m   (2) (x j ,y;Δx/υ) for υ=1,2,3,4 if Δx=ο(π/k 0 ). First, we introduce some new notation. Define the interpolated range steps x μ   (υ)  as follows.
 
x μ   (υ) =μΔx/y for μ=1,2, . . . and υ=1,2, . . .  (9.66) 
 
Note that x υj   (υ)=x   j  for all j and υ.
 
   The procedure requires that the discrete form of problem II be solved four times over each interval x j ≦x≦x j+1  using fractional steps Δx/υ for υ=1,2,3,4. For each υ we must replace Δx in eqs.(9.45,9.49) with Δx/υ and solve this problem at the interpolated range steps x υj+μ   (υ)  for μ=1, . . . ,υ. If we use the Chebyshev tau method to do the cross range integration in eq.(9.45), then we can do these calculations in terms of the Chebyshev expansion coefficients of the solutions {overscore (ψ)} m   (2) (x υj+μ   (υ) ,y;Δx/υ) without having to evaluate the solutions at any cross range y. These solutions require ten times as much work to compute as {overscore (ψ)} m   (2) (x j+1 ,y;Δx) alone, but this computational expense can be justified by the great increase in range over which the extrapolated solutions are accurate. 
   We denote these estimates by  Δx   (4) ∘{overscore (ψ)} m   (2) (x j+1 ,y). They are computed as follows.
 
 Δx   (4) ∘{overscore (ψ)} m   (2) (x j+1 ,y)=[8195{overscore (ψ)} m   (2) (x j+1 ,y;Δx/4)−6561{overscore (ψ)} m   (2) (x j+1 ,y;Δx/3)
 
+896{overscore (ψ)} m   (2) (x j+1 ,y;Δx/2)−7{overscore (ψ)} m   (2) (x j+1 ,y;Δx)]/2520  (9.67)
 
In the modified solution method {overscore (ψ)} m   (2) is replaced by  Δx   (4) ∘{overscore (ψ)} m   (2) . Thus, we use the modified solution method to find the values {overscore (ψ)} m (x j+1 ,y k ) and the extrapolates  Δx   (4) ∘{overscore (ψ)} m   (2) (x j+1 ,y k ) for m=1, . . . ,M and k=0,1, . . . ,N from the values {overscore (ψ)} m (x j ,y l ) and from the extrapolates  Δx   (4) ∘{overscore (ψ)} m   (2) (x j ,y l ) for m=1, . . . ,M and l=0,1, . . . ,N. These quantities are related by eq. (9.37) after we have replaced {overscore (ψ)} m   (2)  with  Δx   (4) ∘{overscore (ψ)} m   (2)  in it.
 
{{overscore (Ψ)} m (x j+1 −0,y k )−k m   Δx   (4) ∘{overscore (ψ)} m   (2) (x j+1 ,y k )} k=0   N 
 
=ℑ N   −1 (exp{−iη 2 Δx/2k m }ℑ N {{overscore (Ψ)} m (x j +0,y l )−k m   Δx   (4) ∘{overscore (ψ)} m   (2) (x j ,y l )} l=0   N ).  (9.68)
 
   We also need the boundary values {overscore (ψ)} m   (1) (x υj+μ   (υ) ,±W) at the interpolated range steps to evaluate the boundary conditions that are satisfied by the solutions {overscore (ψ)} m   (2) (x υj+μ   (υ) ,y;Δx/υ). We compute them by replacing Δx with μΔx/υ in eq.(9.44). Thus, for μ=1, . . . ,υ, 
                 k   m     ⁢         φ   ~     m     (   1   )       ⁡     (       x     vj   +   μ       (   v   )       ,     ±   W       )         =         N     -   1         (       y   k     =     ±   W       )       ⁢       (     exp   ⁢     {       -     in   2       ⁢   μΔ   ⁢           ⁢     x   /     v2k   m         }     ⁢     N     ⁢       {           Ψ   ~     m     ⁡     (         x   j     +   0     ,     y   l       )       -       k   m     ⁢         Δ   ⁢           ⁢   x       (   4   )       ∘         φ   ~     m     (   2   )       ⁡     (       x   j     ,     y   l       )             }       l   =   0     N       )     .               (   9.69   )             
 
Note that it is not necessary to compute the whole inverse DFT to evaluate the RHS of this equation because we need only the values from the inverse DFT at the endpoints y k =±W.
 
   In medium range applications, the lower order Richardson extrapolates
 
 Δx   (2) ∘{overscore (ψ)} m   (2) (x j+1 ,y)=[4{overscore (ψ)} m   (2) (x j+1 ,y;Δx/2)−{overscore (ψ)} m   (2) (x j+1 ,y;Δx)]/3,  (9.70) 
 
 Δx   (3) ∘{overscore (ψ)} m   (2) (x j+1 ,y) =[243{overscore (ψ)} m   (2) (x j+1 ,y;Δx/ 3 )−128{overscore (ψ)} m   (2) (x j+1 ,y;Δx/2)+5{overscore (ψ)} m   (2) (x j+1 ,y;Δx)]/120  (9.71) 
 
can be used in place of  Δx   (4) ∘{overscore (ψ)} m   (2) (x j+1 ,y) to reduce the computational cost of the procedure.
 
10. Ramifications.
 
   Many modifications of the invention are possible in light of the above teachings. For example, it can be adapted to waveguides with multiple internal interfaces. In addition, loss can be incorporated in the model. One way to do this is to modify the modal wavenumbers by a perturbative method that is well-known in the art. Furthermore, a variable step size and variable reference wavenumbers can be used in place of the fixed step size and fixed reference wavenumbers in the preferred embodiment. 
   Additional variations of the invention also are possible in light of the above teachings. Thus, it can be adapted for predicting propagation of seismic waves in the earth&#39;s crust, propagation of radio waves in the earth &#39;s atmosphere, and propagation of light waves in dielectric material. Furthermore, computer implementation of the invention is most efficient on parallel computers because many of the operations in the decomposition method can be done in parallel or can be individually vectorized. It is therefore to be understood that the invention may be practiced otherwise than has been specifically described.