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