From: rusin@vesuvius.math.niu.edu (Dave Rusin)
Newsgroups: sci.math
Subject: Re: Does this have a solution?
Date: 31 Mar 1999 22:14:35 GMT
Keywords: Yet another difficult definite integral

Andrew Sendonaris  <bla@rice.edu> wrote:
> Let f(p) = V*cos(p)
> where V>0
>
>I want to find the following integral:
>
>  Integral of { f(p) * exp( f(p)^2 / 2) * (1-0.5*erfc( f(p)/sqrt(2))) }
>       where p=0..Pi

"Most" integrals admit no closed-form expression, so it's probably not
worth even looking at these, but since no one has responded I took at a
look and found it's easily integrable.

First note that since  erfc(x) = 1- erf(x)  we really have
>  Integral of { f(p) * exp( f(p)^2 / 2) * (0.5+0.5*erf( f(p)/sqrt(2))) }
>       where p=0..Pi
The first part,
>  Integral of { f(p) * exp( f(p)^2 / 2) * (0.5) }
>       where p=0..Pi
is zero by symmetry: the integrals from 0..Pi/2 and Pi/2..Pi are equal in
magnitude but opposite in sign. So we need only deal with the second part.
This time, the symmetry shows it suffices to integral from 0..Pi/2  and double.

Since  erf  is itself defined as an integral, this is now the integral of 
> f(p) * exp( f(p)^2 / 2) * (1/sqrt(Pi)*int((exp(-t^2), t=0..f(p)/sqrt(2))))
>       where p=0..Pi/2

This may be written as a double integral, of
 f * exp( f^2 / 2) * (1/sqrt(Pi)) * exp( -t^2 )  dt dp
where p=0..Pi/2, t=0..f/sqrt(2), and  f  is an abbreviation for  V*cos(p).

We can let  g = V/sqrt(2) sin(p)  and replace  p  in favor of  g  here; 
note that the shorthand  f = V cos(p)  must now be written  
f = sqrt(V^2 - 2 g^2). There is no ambiguity about signs since we have 
restricted to  p < Pi/2. Fortunately,  dg = f/sqrt(2) dp  so in the integrand 
itself we don't need the square roots; it's just part of the limits of 
integration. So we are integrating
  sqrt(2/Pi) * exp( V^2/2 - g^2 - t^2 )  dt dg
over the region with  t>0, g>0, and t^2 + g^2 < V^2/2  -- a quarter circle.

Now we can switch to polar coordinates: if
g = r cos(u)  and  t = r sin(u)  then our integrand becomes
  exp(V^2/2) * sqrt(2/Pi) * exp( -r^2) * r dr du
The integral over  u=0..Pi/2  is then the integral of
  exp(V^2/2) * sqrt(Pi/2) * exp( -r^2) * r dr
but of course  int(exp(-r^2) r ,r=0..R) = (1-exp(-R^2))/2;
with limits  r = 0..V/sqrt(2), this leaves us with 
  exp(V^2/2) * sqrt(Pi/2) * (1 - exp(-V^2/2))/2.

Doubling as indicated before gives the explicit formula

  ( exp(V^2/2) - 1 ) * sqrt(Pi/2)

for the original integral. Sample check:

> I1:=Int(subs(f=V*cos(p),f*exp(f^2/2)*(1-(1/2)*erfc(f/sqrt(2)))),p=0..Pi);
          Pi
         /
        |                      2       2                              1/2
 I1 :=  |    V cos(p) exp(1/2 V  cos(p) ) (1 - 1/2 erfc(1/2 V cos(p) 2   )) dp
        |
       /
         0

> I2:=( exp(V^2/2) - 1 ) * sqrt(Pi/2);                                     
                                        2        1/2   1/2
                    I2 := 1/2 (exp(1/2 V ) - 1) 2    Pi

> evalf(subs(V=.1234,I1));
                                 .009578877572

> evalf(subs(V=.1234,I2));
                                 .009578878170

dave