DDiscrete Bessel and Mathieu functions
Kenan Uriostegui and Kurt Bernardo WolfInstituto de Ciencias F´ısicasUniversidad Nacional Aut´onoma de M´exicoAv. Universidad s/n, Cuernavaca, Morelos 62251, M´exico Abstract
The two-dimensional Helmholtz equation separates in elliptic co-ordinates based on two distinct foci, a limit case of which includespolar coordinate systems when the two foci coalesce. This equationis invariant under the Euclidean group of translations and orthogonaltransformations; we replace the latter by the discrete dihedral groupof N discrete rotations and reflections. The separation of variablesin polar and elliptic coordinates is then used to define discrete Besseland Mathieu functions, as approximants to the well-known continuous
Bessel and Mathieu functions, as N -point Fourier transforms approx-imate the Fourier transform over the circle, with integrals replacedby finite sums. We find that these ‘discrete’ functions approximatethe numerical values of their continuous counterparts very closely andpreserve some key special function relations. The role of the Euclidean group of translations, reflections and rotationsin the determination of the coordinate systems that separate the solutionsof the two-dimensional Helmholtz equation is well known from the work byWillard Miller Jr. [1, Ch. 1]. This symmetry accounts for their separabilityin four coordinate systems: Cartesian, polar, parabolic and elliptic. Only theelliptic system is generic; when the two foci coalesce, this system becomesthe polar one with angular and radial coordinates; when one focus departsto infinity the system becomes parabolic; and when both foci do, it becomesCartesian. Posgrado en Ciencias F´ısicas, Universidad Nacional Aut´onoma de M´exico a r X i v : . [ m a t h - ph ] F e b he polar decomposition was used by Biagetti et al. [2] to first introducea discrete version of Bessel functions based on an expansion of plane wavesinto a finite number of polar components —that was not quite complete.This was properly completed in Ref. [3], defining discrete Bessel functions B N n ( ρ ), which approximate the usual continuous Bessel functions J n ( ρ ) byreplacing Fourier series over a circle S by the finite Fourier transform on N equidistant points on that circle, θ m = 2 πm/N, m ∈ { , , . . . , N − } =: S ( N ) , (1)where m is counted modulo N . It was found that these discrete functionsapproximated very closely (of the order 10 − ) the corresponding continuousones over a region, roughly 0 ≤ n + ρ < N .Several authors have introduced functions that approximate the well-known continuous Bessel functions J n ( ρ ), for the purpose of reducing com-putation time, or to provide new classes of solutions to difference equationsthat will share some of their salient properties [4, 5, 6]. Our approach followsthe well known approximation afforded by the N -point finite Fourier trans-form to the integral Fourier transform over the circle. This is done for polarand elliptic coordinates, and introduces both ‘discrete’ Bessel and Mathieufunctions. These functions, we should emphasize, differ from those proposedin the works cited above, which are also distinct in definition and purposeamong themselves. By construction it will follow that under N → ∞ , thesediscrete functions become the continuous ones, although this limit requiresfurther mathematical precision, as it may involve Gibbs-type oscillation phe-nomena that we cannot address here.In Sect. 2 we present this discretization method and a resum´e of theresults in Ref. [3] for Bessel functions, to note that the discrete functionsthus defined approximate the continuous ones remarkably well. In Sect. 3 ofthe present paper we apply the strategy of replacing harmonic analysis on S by S ( N ) to define discrete approximants to the Mathieu functions of firstand second kind in the elliptic coordinate system. All relations are backedby numerical verification. In the concluding Sect. 4 we provide some furtherconnections and preliminary conclusions.2 Continuous and discrete Bessel functions
The Helmholtz equation for wavefields f ( x, y ) of (fixed) real wavenumber κ ∈ R , is ( ∂ x + ∂ y + κ ) f ( x, y ) = 0 , (2)with ∂ z ≡ ∂/∂ z and ( x, y ) ∈ R . In this section we follow the well knowncase of polar coordinates, x = r cos θ, y = r sin θ, r ∈ [0 , ∞ ) , θ ∈ ( − π, π ] = S . (3)A key assumption is a Hilbert space structure for the solutions f ( x, y ) bywhich one can write them as the two-dimensional Fourier transform, f ( x, y ) = 12 π (cid:90)(cid:90) R d κ x d κ y exp i( xκ x + yκ y ) (cid:101) f ( κ x , κ y ) . (4)The Helmholtz equation (2) is then correspondingly transformed to a conju-gate space ( κ x , κ y ) ∈ R where it reads ( κ − κ x − κ y ) (cid:101) f ( κ x , κ y ) = 0, whichwe can also refer to polar coordinates κ x = κ cos φ , κ y = κ sin φ , with thesurface element d κ x d κ y = κ d κ d φ . The solutions to the Fourier-transformedHelmholtz equation are thus reduced by a Dirac δ -distributions in the ra-dius [1, Ch. 1], as (cid:101) f ( κ x , κ y ) = √ πκ − δ ( κ − (cid:101) κ ) f ◦ ( φ ), with f ◦ ( φ ) a functionon the φ -circle S of radius (cid:101) κ , that we write again κ , understanding thatit is the fixed wavenumber. The Helmholtz solutions (4) thus acquire thesingle-integral form f ( x, y ) = 1 √ π (cid:90) S d φ exp i κ ( x cos φ + y sin φ ) f ◦ ( φ ) , (5)with the Hilbert space structure based on the inner product of functions f (1) ◦ ( φ ) and f (2) ◦ ( φ ) on the circle,( f (1) ◦ , f (2) ◦ ) ◦ := (cid:90) S d φ f (1) ◦ ( φ ) ∗ f (2) ◦ ( φ ) . (6)It is here that we reduce the continuous circle Fourier transform to the N -point discrete Fourier transform, from S to S ( N ) , replacing integrals bysummations and the continuous variable φ ∈ S with φ m ∈ S ( N ) , as (cid:90) S d φ F ◦ ( φ ) ↔ (cid:88) m ∈S N ) F ( φ m ) , π ↔ N,φ m = 2 πm/N, (7)3or m ∈ { , , . . . , N − } counted modulo N ; the set of N discrete angles φ m are thus equidistant by 2 π/N . The functions f ( φ m ) ≡ f m can be interpretedas sample points of a continuous function, or as the index for the list ofcomponents of an N -cyclic vector. In either case, the inner product of twodiscrete functions f (1) n and f (2) n is naturally( f (1) , f (2) ) ( N ) := N − (cid:88) n =0 f (1) ∗ n f (2) n , (8)and it is clear that the N → ∞ limit will lead back from the discrete to thecontinuum, with the approximations and limits familiar from Fourier theory.The Helmholtz equation (2) in polar coordinates, multiplied by r ,( r ∂ r + r∂ r + ∂ φ + κ ) f ( r, φ ) = 0 , (9)shows that solutions can be factored into a function of the radius times afunction of the angle as f ( r, φ ) = R ( r ) Φ( φ ), while (5) implies that solutionsΦ( φ ) for the angular factor will determine a corresponding radial factor R ( r ).An orthonormal and complete set of eigenfunctions of ∂ φ over the circle φ ∈ S is the set of phases Φ n ( φ ) := (2 π ) − / exp(i nφ ), with integer n ∈{ , ± , . . . } , and inner products (Φ n , Φ n (cid:48) ) ◦ = δ n,n (cid:48) . When the domain ofthese functions is restricted from φ ∈ S to φ m ∈ S ( N ) as in (1), we retainthe subset of N functions on the N points in S ( N ) , given byΦ ( N ) n ( φ m ) := 1 √ N exp(i nφ m ) = 1 √ N exp (cid:18) π i mnN (cid:19) = Φ ( N ) n ± N ( φ m ) , (10)labeled by the cyclic subset n ∈ { , , . . . , N − } , that are also orthonor-mal under the common inner product (8) for discrete functions on S ( N ) , andcomplete:(Φ ( N ) n , Φ ( N ) n (cid:48) ) ( N ) = δ n,n (cid:48) , N − (cid:88) n =0 Φ ( N ) n ( φ m ) ∗ Φ ( N ) n ( φ m (cid:48) ) = δ m,m (cid:48) . (11)Returning to (5) with ( x, y ) in the polar coordinates ( r, θ ) of (3), and takingfor f ◦ ( φ m ) the basis functions (10) on the discrete points of S ( N ) , we writethe N solutions to the discretized Helmholtz equation, labeled by cyclical4 ∈ { , , . . . , N − } , as f n ( r, θ k ) = 1 √ N (cid:88) m ∈S N ) exp[i κr (cos θ k cos φ m + sin θ k sin φ m )] Φ ( N ) n ( φ m )= 1 N (cid:88) m ∈S N ) exp[i κr cos( θ k − φ m )] exp(i nφ m )= e i n ( θ k + π/ N (cid:88) m ∈S N ) exp(i κr sin ϕ m ) exp( − i nϕ m ) , (12)having replaced ϕ m := θ k − φ m + π in the summation over the N discretepoints on the circle.Following Miller [1, p. 29], the phase in front of (12), e i nθ n ( θ k + π/ =i n e π i nk/N = i n √ N Φ ( N ) n ( θ k ), is extracted to write the functions as f n ( r, θ k ) = i n √ N B ( N ) n ( κr ) Φ ( N ) n ( θ k ) , (13)where the radial factor B ( N ) n ( ρ ), ρ := κr , are the discrete Bessel functions .From (12) these functions are seen to be real and their parities, using co-efficients { c n , s n } := { , } for n even or { , } for n odd, can be writtenas B ( N ) n ( ρ ) = 1 N (cid:88) m ∈S N exp(i ρ sin ϕ m ) [ c n cos( nϕ m ) − i s n sin( nϕ m )]= 1 N (cid:88) m ∈S N exp(i ρ sin ϕ m ) × (cid:40) cos nϕ m , n even , − i sin nϕ m , n odd . (14)The distinction between even and odd cases of n , as done in [3], is subtle butimportant to obtain the correct result for all n ’s ( cf . [2, Eq. (9)]). It resultsin the parity and cyclicity properties B ( N ) n ( ρ ) = B ( N ) n ± N ( ρ ) = ( − n B ( N ) − n ( ρ ) = ( − n B ( N ) n ( − ρ ) , B ( N ) n (0) = δ n, , (15)which also hold for the continuous Bessel functions J n ( ρ ) of integer order [8].A plane wave of wavenumber κ along the y -axis in a Helmholtz mediumthat allows only N equidistant directions of propagation on the circle, canbe obtained from (14) using the completeness relation (11) to expand the5iddle term and writeexp(i ρ sin ϕ m ) = B ( N ) ( ρ ) + 2 N − (cid:88) n =1 B ( N ) n ( ρ ) cos(2 nϕ m )+ 2i N − (cid:88) n =0 B ( N ) n +1 ( ρ ) sin((2 n +1) ϕ m ) , (16)showing how the discrete Bessel functions can take the place of the continuousones, cf. [7, Eq. KU120(13)].In Ref. [3] we proved analytically, and verified numerically, that the fol-lowing expressions for the discrete Bessel functions are exact analogues ofthose valid for continuous Bessel functions. Corresponding to [7, WA44] forodd N =: 2 j + 1, in Ref. [3] we proved the linear relations involving the evenand odd- n discrete Bessel functions, B ( ρ ) + j (cid:88) n =1 B n ( ρ ) cos(2 nϕ m ) = cos( ρ sin ϕ m ) , (17) j (cid:88) n =0 B n +1 ( ρ ) sin((2 n +1) ϕ m ) = sin( ρ sin ϕ m ) . (18)The quadratic formulas [9, § D functions, under contraction from the rotation to theEuclidean group. These relations, of group-theoretical origin, retain theirvalidity under the discretization of the rotation subgroup, and lead to j (cid:88) n = − j B n ( ρ ) B n (cid:48) − n ( ρ (cid:48) ) = B n (cid:48) ( ρ + ρ (cid:48) ) , (19)keeping in mind the parity property (15) for the negative n -indices in thesum for odd N , addressing the vector rather than spin representations of therotation group.In Fig. 1 we essentially repeat the figure in Ref. [3] where we comparedthe discrete and continuous Bessel functions, B ( N ) n ( ρ ) and J n ( ρ ), to supportthe claim that the approximation is indeed remarkable within an intervalthat is roughly 0 ≤ n + ρ < N . A similar set of figures is presented belowfor Mathiew functions.Now, having N basis functions B ( N ) n ( ρ ), numbered by cyclic n modulo N , it is natural to inquire whether the argument ρ can or should be also6igure 1: The ‘discrete’ Bessel functions B ( N ) n ( ρ ) on continuous intervals 0 ≤ ρ ≤ (2 N − vs . the ‘continuous’ Bessel functions J n ( ρ ) (thin black lines), for orders n ∈{ , , , } and point numbers N ∈ { , , } . Heavy black lines replace bothwhere the ‘discrete’ and the ‘continuous’ Bessel functions differ by less than 10 − . N integer values ρ k = k ∈ { , , . . . , N − } . This wasdone in Ref. [2] while in [3] the plot in Fig. 1 marked these points and usedthem to define a kernel B ( N ) n ( ρ k ) for a ‘discrete Bessel transform’ betweentwo N -vectors of components f n and (cid:101) f k . The fact is that while the angle ϕ isdiscretized naturally to N points on the circle, the radial coordinate ρ is not subject to a similarly compelling set of points, but is valid and non-cyclicover the complex ρ -plane. The same discretization process for the angular—but not the radial— coordinate will be applied to the Mathieu case below. Elliptic coordinates on the plane generalize the previous polar coordinates(3); they are defined in terms of Cartesian coordinates through x = cosh (cid:37) cos ψ, y = sinh (cid:37) sin ψ, (cid:37) ∈ [0 , ∞ ) , ψ ∈ ( − π, π ] = S . (20)where ( (cid:37), ψ ) are analogues of the previous polar coordinates ( r, φ ) for whichwe retain the names as ‘radial’ and ‘angular’ variables. For fixed (cid:37) or forfixed ψ , the locus of points ( x, y ) ∈ R that satisfy x / cosh (cid:37) + y / sinh (cid:37) = 1 , x / cos ψ − y / sin ψ = 1 , (21)draw families of confocal ellipses or hyperbolas respectively. At (cid:37) = 0, ψ ∈ S draws twice the line between the two foci ( x, y ) = ( ± ,
0) for ψ = (0 , π ). Themajor and minor semi-axes of the ellipses are cosh (cid:37) and sinh (cid:37) respectively,so their eccentricities are 1 / cosh (cid:37) , that tend to circles when (cid:37) → ∞ . Onthe other hand, for fixed ψ ∈ S in each of the four quadrants, since (cid:37) ≥ four parity cases out of the two reflections, across the x and y axes. Comparethis with the case of polar coordinates where r ≥ − n in (14) provides the two Bessel parity cases.The Helmholtz differential equation (2), written in the elliptic coordinates(20), is clearly separable,[( ∂ (cid:37) + κ cosh (cid:37) ) + ( ∂ ψ − κ cos ψ )] f ( (cid:37), ψ ) = 0 , (22)so that solutions can be written in the product form f ( (cid:37), ψ ) ∼ P ( (cid:37) ) Ψ( ψ ).Dividing by f , one obtains two coupled equations in (cid:37) and ψ , the latter is8igure 2: Set of equally-spaced discrete points on ellipses (20) of the ‘angular’ coordinates { ψ m } ∈ S ( N ) for N = 21, and hyperbolas of the ‘radial’ coordinate for (cid:37) ∈ { . , , . } . an eigenvalue equation in the angular coordinate,( ∂ ψ − q cos 2 ψ ) Ψ( ψ, q ) = ν Ψ( ψ, q ) , q := κ , (23)known as the Mathieu differential equation. The angular coordinate ψ isperiodic and a well-known solution method consists in expanding solutionsof (23) in the Fourier basis ∼ exp(i nψ ) over all integer n . This defines theMathieu functions of the first kind ce n ( ψ, q ) and se n ( ψ, q ) with integer n [11],characterized by a parity index p ∈ { , } for even and odd cases [7, Eqs.8.61], and distinct for even and odd indices. In a two-line expression all casescan be written as (cid:20) ce n + p ( ψ, q )se n + p ( ψ, q ) (cid:21) = ∞ (cid:88) s =0 (cid:20) A n + p s + p cos((2 s + p ) ψ ) B n + p s + p sin((2 s + p ) ψ ) (cid:21) . (24)The parities are even ce n ( − ψ, q ) = ce n ( ψ, q ), odd se n ( − ψ, q ) = − se n ( ψ, q ),and se ( ψ, q ) ≡
0. The coefficients A ns , B ns are found introducing this expan-sion into (23) to find recursion relations [7, Eqs. 8.62] that lead to efficientnumerical computation. For use below, we write them using Fourier series as (cid:20) A ns B ns (cid:21) = 1 π (cid:90) S d ψ (cid:20) cos( sψ ) ce n ( ψ, q )sin( sψ ) se n ( ψ, q ) (cid:21) , (25)for n (cid:54) = 0, while A n = (2 π ) − (cid:82) S d ψ ce n ( ψ, q ), and B n ≡
0. The Mathieufunctions (24) are orthogonal under the inner product (6) over the circle,9ce m , ce n ) ◦ = πδ m,n , (se m , se n ) ◦ = πδ m,n for n (cid:54) = 0 —zero otherwise, and(ce m , se n ) ◦ = 0.We now restrict the range of the angular coordinate ψ from S to S ( N ) ,shown for the elliptic coordinates in Fig. 2, in correspondence with the pre-vious discrete phase functions in the Bessel case (10), and thus defining the‘angular’ discrete Mathieu functions of the first type over ψ m ∈ S ( N ) as (cid:20) ce ( N ) n + p ( ψ m , q )se ( N ) n + p ( ψ m , q ) (cid:21) := N − (cid:88) s =0 (cid:20) a n + p s + p cos((2 s + p ) ψ m ) b n + p s + p sin((2 s + p ) ψ m ) (cid:21) , (26)with coefficients a ns , b ns . The finite N -point Fourier transform approximatesthem through the replacement (7) to the functions and coefficients A ns , B ns of the continuous case in (25), as (cid:20) a ns b ns (cid:21) := 1 N N − (cid:88) m =0 (cid:20) cos( sψ m ) ce ( N ) n ( ψ m , q )sin( sψ m ) se ( N ) n ( ψ m , q ) (cid:21) (cid:39) (cid:20) A ns B ns (cid:21) , (27)for s (cid:54) = 0, while a n = A n , b n = 0, and also b s = 0.The last relation in (27) is an approximate equality, the validity of whichis contingent upon the numerical computation and comparison between thelower- and upper-case coefficients within a range of their indices in, say,0 ≤ n, s ≤ N −
1, which is reflected in turn by the discrete and continuousMathieu functions themselves. In Fig. 3 we compare a sample of continuousangular Mathieu functions of the first kind with their discrete approximationsfrom Eq. (26). In favor of the thus defined discrete Mathieu functions, wenote that they satisfy orthogonality relations under the discrete inner product(8), namely(ce ( N ) n , ce ( N ) n (cid:48) ) ( N ) = N δ n,n (cid:48) , (se ( N ) n , se ( N ) n (cid:48) (cid:54) =0 ) ( N ) = N δ n,n (cid:48) , (ce ( N ) n , se ( N ) n (cid:48) ) ( N ) = 0 . (28)By construction, the parities of the discrete Mathieu functions are also evence ( N ) n ( − ψ m , q ) = ce ( N ) n ( ψ m , q ), or odd se ( N ) n ( − ψ m , q ) = − se ( N ) n ( ψ m , q ).At this point it is illuminating to inquire into the manner in which thediscrete functions approximate the continuous ones. Consider for examplehow ce ( N ) ( ψ m , q ), whose definition (26) allows us to compute it for continuous ψ m ∈ S , matches ce ( ψ, q ) in the whole ψ range. In Fig. 4 we do so for small N , noting that where the continued ψ m lines of the former take theirvalues, they intersect the properly continuous line of the latter; although thetwo lines intersect also at other points, the two lines remain notably distinct.10igure 3: Discrete vs. continuous ‘angular’ Mathieu functions for N = 41, q = 2. Thevalues of the discrete functions ce ( N ) n ( ψ m , q ) and se ( N ) n ( ψ m , q ) at ψ m , 0 ≤ m ≤ N −
1, areindicated by circles. The continuous Mathieu functions ce n ( ψ, q ) and se n ( ψ, q ) are markedby black lines in their full range 0 ≤ ψ < π . Their difference is less than 10 − for allpoints ψ m . As the figure shows, the approximation is not valid over presumably smallranges around these intersections, but only at the prescribed ψ m = 2 πm/N points. We intend to elaborate on such and similar limits elsewhere.Proceeding now as we did in (12), but using the discrete Mathieu func-tions of the first kind ce ( N ) n ( ψ m , q ) and se ( N ) n ( ψ m , q ) in place of the plain phasefunctions Φ ( N ) n ( φ m ), we again have Helmholtz solutions f n ( (cid:37), ψ k ) on the planethat are characterized by parities and sub-indices n , whose radial factor willbe the discrete ‘radial’ Mathieu functions of the second kind, to be indicatedcorrespondingly as Ce ( N ) n ( (cid:37), q ) and Se ( N ) n ( (cid:37), q ), (cid:20) f c n + p ( (cid:37), ψ k ) f s n + p +1 ( (cid:37), ψ k ) (cid:21) = 1 N N − (cid:88) m =1 (cid:20) ce ( N ) n ( ψ m , q )se ( N ) n ( ψ m , q ) (cid:21) exp[i κ ( x cos ψ m + y sin ψ m ) (29)=: (cid:20) c n ( q ) Ce ( N ) n ( (cid:37), q ) ce ( N ) n ( ψ k , q ) s n ( q ) Se ( N ) n ( (cid:37), q ) se ( N ) n ( ψ k , q ) (cid:21) , (30)11igure 4: Comparison between the ‘discrete’ Mathieu functions ce ( N ) ( ψ m , q ) whose argu-ments are continued to ψ m ∈ S (gray line) vs . the ‘continuous’ Mathieu function ce n ( ψ, q )(black line), for point numbers N ∈ { , , } and q = 2. The discrete points ψ m ∈ S ( N ) lie at a subset of the intersections marked with circles. where c n ( q ) and s n ( q ) are constants. Using the elliptic coordinates with adiscretized angular part, x ( (cid:37), ψ k ) = cosh (cid:37) cos ψ k and y ( (cid:37), ψ k ) = sinh (cid:37) sin ψ k in (20), the phase exponent is then i κ = 2i √ q times x cos ψ m + y sin ψ m =cosh (cid:37) cos ψ k cos ψ m + sinh (cid:37) sin ψ k sin ψ m . As was done before in (12) and(13), we extract the new discrete ‘radial’ functions using the orthogonality(28) of the previous discrete ‘angular’ Mathieu functions, as (cid:20) Ce ( N ) n + p ( (cid:37), q )Se ( N ) n + p +1 ( (cid:37), q ) (cid:21) = (cid:20) /N c n + p ( q ) ce ( N ) n + p ( ψ k , q )1 /N s n + p +1 ( q ) se ( N ) n + p +1 ( ψ k , q ) (cid:21) N − (cid:88) m =0 (cid:20) ce ( N ) n + p ( ψ m , q )se ( N ) n + p +1 ( ψ m , q ) (cid:21) × exp[2i √ q (cosh (cid:37) cos ψ k cos ψ m + sinh (cid:37) sin ψ k sin ψ m )] . (31)The coefficients in front of the summation will be now determined throughconsidering specific values for the ‘angular’ coordinate ψ ↔ ψ m , comparingthem with expressions of the continuous Mathieu functions of the secondkind obtained from integrals that are tabulated in Ref. [7, § ψ m = 0 or π , althoughthe latter is not in the set S ( N ) if ψ = 0, since N was assumed to be odd.Let us first consider the case of even parity p = 0 and the angle ψ = π in (31), where the previous remark applies. Based on the close approxima-tion between the discrete and continuous Mathieu functions, we may simplyreplace the latter for the former, so that the two lines in that expression read (cid:20) Ce ( N ) n ( (cid:37), q )Se ( N ) n +1 ( (cid:37), q ) (cid:21) = (cid:20) K c n K s n +1 (cid:21) N − (cid:88) m =0 (cid:20) ce ( N ) n ( ψ m , q )se ( N ) n +1 ( ψ m , q ) (cid:21) exp(2i √ q sinh (cid:37) sin ψ m ) . (32)When this summation formula is compared with the integral expressions12abulated in [7, § (cid:20) Ce n ( (cid:37), q )Se n +1 ( (cid:37), q ) (cid:21) = (cid:20) ce n (0 , q ) / π A n − ise (cid:48) n +1 (0 , q ) / π B n +11 √ q (cid:21) × (cid:90) S d ψ (cid:20) ce n ( ψ, q )se n +1 ( ψ, q ) (cid:21) exp(2i √ q sinh (cid:37) sin ψ ) , (33)we conclude that the constants in the summation (32), after identifying 2 π ↔ N , A n = a n and B n +11 (cid:39) b n +11 , are K c n = ce n (0 , q ) a n N , K s n +1 = − i se (cid:48) n +1 (0 , q )2 b n +11 N √ q . (34)where se (cid:48) n (0 , q ) := d se n ( ψ, q ) / d ψ | ψ =0 . In Fig. 5 we compare a sample of thediscrete and continuous ‘radial’ Mathieu functions, noting that the two linesare quite coincident in the range (cid:37) ∈ [0 , π ), but that the discrete approximantoscillates wildly beyond π . Again, here we can only justify this statementnumerically.Next we consider the case of odd parity p = 1 at the value ψ = ψ = 0.The upper line in (32) readsCe ( N ) n +1 ( (cid:37), q ) = K c n +1 N − (cid:88) m =0 ce ( N ) n +1 ( ψ m , q ) exp(2i √ q cosh (cid:37) cos ψ m ) , (35)that we compare with the integral for the continuous Mathieu functions ofthe second kind in [7, § n +1 ( (cid:37), q ) = i ce (cid:48) n +1 ( π, q )2 πA n +11 √ q (cid:90) S d ψ ce n +1 ( ψ, q ) exp(2i √ q cosh (cid:37) cos ψ ) , (36)where ce (cid:48) n +1 ( ψ, q ) is the derivative of the Mathieu function. Again exploitingthe correspondences (7), 2 π ↔ N and A n +11 (cid:39) a n +11 , we conclude theconstant in (38) to be K c n +1 = i ce (cid:48) n +1 ( π, q )2 a n +11 N √ q . (37)The remaining case to be considered is that of odd parity and even index,namely for Se ( N ) n +2 ( (cid:37), q ). This presents a problem though, because the sum-mation (31) is identically zero for both ψ k = 0 and π due to the parities of13igure 5: Discrete vs. continuous ‘radial’ Mathieu functions in the interval 0 ≤ (cid:37) < . N ∈ { , , } and here for q = 2. The ‘discrete’ functions Ce ( N ) n ( (cid:37), q ) and Se ( N ) n ( (cid:37), q )with the (continuous) argument (cid:37) (gray line), is compared with the ‘continuous’ functionsCe n ( (cid:37), q ) and Se n ( (cid:37), q ) (thin black line). As before, where both coincide within 10 − theyare replaced by a thick black line. The radial Mathieu functions, when computed with thecommercial Mathematica algorithm, oscillate wildly after an upper value that decreaseswith increasing values of q . the terms in the sum. It is different from zero for 0 < ψ k < π however, so ifwe choose ψ k = π , where both summands in the exponent appear as 1 / √ ( N ) n +2 ( (cid:37), q ) = K s n +2 N − (cid:88) m =0 se ( N ) n +2 ( ψ m , q ) × exp[i (cid:113) q (cosh (cid:37) cos ψ m + sinh (cid:37) sin ψ m )] . (38)For the corresponding continuous case, we could not find a correspondingintegral in [7, § K s n +2 in (38). The lack of a similar plane-wave integral expression for thecontinuous Mathieu functions Se n +2 ( (cid:37), q ) has been noted also in Ref. [12]without explanation. However, we have checked numerically that the simileof the discrete to continuous functions approximation provided bySe ( N ) n +2 ( (cid:37), q ) (cid:39) − i se n +2 (i (cid:37), q ) = Se n +2 ( (cid:37), q ) , (39)which is an equality for continuous functions, cf . [7, 8.611.4, 8.631.4]. For0 < (cid:37) < − . We should note that gener-ally the discrete ‘radial’ and ‘angular’ Mathieu functions for pure imaginary14rguments are not related to similar equalities of their continuous integralexpressions, because the summation definitions in (26) involve hyperbolicfunctions. In particular, say,Se ( N ) n + p +1 ( (cid:37), q ) (cid:54) = − ise ( N ) n + p +1 (i (cid:37), q ) = N − (cid:88) m =0 b n + p +12 m + p +1 sinh[(2 m + p +1) (cid:37) ] . (40) The expansion of Helmholtz plane waves in series of radial Bessel and angu-lar trigonometric functions has its discrete analogue in Eq. (16), which tellsus that the wavefield due to a finite number N of plane waves at equidistantdirection angles can be expanded in discrete Bessel radial functions and cor-responding trigonometric angular functions. A similar statement will holdwhen the wavefield is expanded in discrete Mathieu functions with the phasesdetermined by the points on an ellipse as depicted in Fig. 2. Conceivablysuch fields can be produced in resonant two-dimensional micro-cavities fedby a number of activation channels.We recognize that the full treatment and exploration of properties for thediscrete Bessel and Mathieu function presented here is not exhaustive, butthat it should be sufficient to indicate that the approximation method con-sisting in the replacement of a continuous closed subgroup of the symmetrygroup of a partial differential equation by a finite discrete group is definitelyof interest. In the present case of two dimensions, the orthogonal group wasreduced to the dihedral group. In three dimensions, the symmetry Euclideansymmetry group could reduce its three-dimensional rotation subgroup by anyof its polyhedral subgroups, whose functions may serve to describe wavefieldswith a corresponding subset of wave propagation directions. Here we havepresented a set of exact relations, others whose approximation closeness wasestimated through numerical computation, and others that have been onlysuggested by that approach, and for which we expect to present further re-sults from ongoing work. Acknowledgments
We thank the support of the Universidad Nacional Aut´onoma de M´exicothrough the PAPIIT-DGAPA project AG100120 ´Optica Matem´atica .15 eferences [1] W. Miller Jr.,
Symmetry and Separation of Variables , Encyclopedia ofMathematics, Vol. 4 (Cambridge Univerity Press, 1984).[2] G. Biagetti, P. Crippa, L. Falaschetti, and C. Turchetti, Discrete Besselfunctions for representing the Class of Finite Duration Decaying Se-quences, European Signal Analysis Conference, pp. 2126–2130 (Bu-dapest, 2016).[3] K. Uriostegui and K.B. Wolf, Discrete Bessel functions and transform,(submitted) arXiv:2005.06076 [math-ph][4] R.H. Boyer, Discrete Bessel functions,
J. Math. Anal. Appl. , 509–524(1961).[5] M. Bohner and T. Cuchta, The Bessel difference equation, Proc. Amer.Math. Soc. , 1567–1580 (2017).[6] A. Slav´ık, Discrete Bessel functions and partial differential equations,
J.Diff. Eqs. Applics.
DOI:10.1080/10236198.2017.141610 (2017).[7] I.S. Gradshteyn and I.M. Ryzhik,
Tables of Integrals, Series and Prod-ucts , A. Jeffrey and D. Zwillinger Eds, (Academic Press, 2007).[8] G.N. Watson,
Theory of Bessel Functions (Cambridge University Press,1922).[9] A. Erd´elyi et al. , Higher Transcendental Functions (Based on notes byH. Bateman) Vol. 2 (McGraw-Hill, New York, 1953).[10] P. Winternitz, K.B. Wolf, G.S. Pogosyan, and A.N. Sissakian, Graf’saddition theorem obtained from SO(3) contraction,
Theor. Mat. Phys. , 1501–1503 (2001).[11] N.W. McLachlan,
Theory and Application of Mathieu Functions (OxfordUniversity Press, 1947).[12] L. Chaos-Cador and E. Ley-Koo, Mathieu functions, matrix evaluation,and generating functions,
Rev. Mex. F´ıs.48