Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Mathematica can't do this double integral

  • To: mathgroup at smc.vnet.net
  • Subject: [mg48906] Re: Mathematica can't do this double integral
  • From: Paul Abbott <paul at physics.uwa.edu.au>
  • Date: Wed, 23 Jun 2004 02:50:59 -0400 (EDT)
  • Organization: The University of Western Australia
  • References: <cb0ufp$r4a$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

In article <cb0ufp$r4a$1 at smc.vnet.net>,
 Enrique Aguado <l.e.aguado at leeds.ac.uk> wrote:

> It looks like this:
> 
> Int[Int[E^(a Cos[x]+ b Cos[y]+ k Cos[x - y]) {y, -Pi, Pi}],{x, -Pi, Pi}]
> 
> Any suggestions anyone?

Although Mathematica fails to compute it, the inner integral (either 
with respect to x or y, if you change the order of integration) can be 
computed in closed form.

Using an integral representation of BesselI -- see 

  http://functions.wolfram.com/BesselAiryStruveFunctions/BesselI/07/01/01/

or Abramowitz and Stegen 9.6.16 at 

  http://jove.prohosting.com/~skripty/page_376.htm, 

-- we have

 Integrate[Exp[a Cos[t]], {t, 0, Pi}] == Pi BesselI[0, a]

(the Mathematica integrator gets this wrong in 5.0!). From the symmetry 
of the integrand, we can extend this to an integral over [-Pi,Pi] as

 Integrate[Exp[a Cos[t]], {t,-Pi, Pi}] == 2 Pi BesselI[0, a]       (1)

Hence, if k=0, the double integral is separable and reduces to

 Integrate[Exp[a Cos[x] + b Cos[y]],{y, -Pi, Pi},{x, -Pi, Pi}] ==
  (2 Pi)^2 BesselI[0,a] BesselI[0,b]

Moreover, if a=0, one obtains

 Integrate[Exp[b Cos[y] + k Cos[x - y]],{y, -Pi, Pi},{x, -Pi, Pi}] ==
  (2 Pi)^2 BesselI[0,b] BesselI[0,k]

with an analogous result for b=0.

Since the integral in (1) is over one period of the integrand, we see 
that

 Integrate[Exp[a Cos[t+d]], {t,-Pi, Pi}] == 2 Pi BesselI[0, a]

for arbitrary real d. Moreover, we also have that

 Integrate[Exp[a Sin[t+d]], {t,-Pi, Pi}] == 2 Pi BesselI[0, a]

since Sin and Cos are related by a phase shift of Pi/2. 

Using TrigExpand on the argument of the exp in the given integral, 
a Cos[x] + k Cos[x - y] + b Cos[y], and collecting terms, leads to

 a Cos[x] + (b + k Cos[x]) Cos[y] + k Sin[x] Sin[y]

which can be expressed as

 a Cos[x] + Sqrt[b^2 + 2 b k Cos[x] + k^2] Sin[d + y]

where d = ArcTan[b + k Cos[x], k Sin[x]]. Hence 

 Integrate[E^(a Cos[x]+ b Cos[y]+ k Cos[x - y]), {y, -Pi, Pi}] ==
  E^(a Cos[x]) 2 Pi BesselI[0, Sqrt[b^2 + 2 b k Cos[x] + k^2]]

(Note that d does not appear in this expression). If one does the 
integral over x first, an analogous expression arises with a and b 
interchanged. Hence we have, at least, reduced the double integral to a 
single integral.

One approach for computing the remaining integral is to expand the 
BesselI into a series, leading to integer powers of 

  (b^2 + 2 b k Cos[x] + k^2), 

and, since

  BesselI[n, a] == (1/Pi) Integrate[E^(a Cos[x]) Cos[n x], {x, 0, Pi}]

one obtains an infinite sum over BesselI[n, a] with coefficients that 
are polynomials in b and k. However, this is unsymmetrical. A second 
approach is to expand the argument of the integral into a Maclaurin 
series in k and compute each integral in closed form. In this way, I 
have managed to show that the double integral has the following nice 
symmetrical expression:

 (2Pi)^2 (BesselI[0, a] BesselI[0, b] BesselI[0, k] +       
  2 Sum[BesselI[n, a] BesselI[n, b] BesselI[n, k], {n, 1, Infinity}])

Clearly this reduces to the correct answer when one of a, b, or k is 0. 
Limited numerical testing shows that it looks ok. It appears to be a 
generalization of the identity 

 functions.wolfram.com/BesselAiryStruveFunctions/BesselI/23/01/0008/

Cheers,
Paul

-- 
Paul Abbott                                   Phone: +61 8 9380 2734
School of Physics, M013                         Fax: +61 8 9380 1014
The University of Western Australia      (CRICOS Provider No 00126G)         
35 Stirling Highway
Crawley WA 6009                      mailto:paul at physics.uwa.edu.au 
AUSTRALIA                            http://physics.uwa.edu.au/~paul


  • Prev by Date: Re: A module to write a module
  • Next by Date: multigraphics
  • Previous by thread: Re: Mathematica can't do this double integral
  • Next by thread: Re: Mathematica can't do this double integral