From: gjfee@daisy.uwaterloo.ca (Greg Fee)
Subject: Re: Log integral
Organization: University of Waterloo
Date: Thu, 5 Nov 1992 06:06:17 GMT


>                Integrate[Log[(1 - (bx/(a+cx^2))]/x,{x,0,Infinity}]

In Maple form, this integral is:
                                             b x
                          infinity ln(1 - --------)
                             /                   2
                            |             a + c x
                            |      ---------------- dx
                            |              x
                           /
                           0

We can make Maple generate the solution.   In fact, we can even get Maple
to generate the form of the indefinite integral as well.


    |\^/|      MAPLE V
._|\|   |/|_.  Copyright (c) 1981-1990 by the University of Waterloo.
 \  MAPLE  /   All rights reserved.  MAPLE is a registered trademark of
 <____ ____>   Waterloo Maple Software.
      |        Type ? for help.

--> we want to solve:
> Int(log(1 - b*x/(a+c*x^2))/x, x = 0..infinity);
                                             b x
                          infinity ln(1 - --------)
                             /                   2
                            |             a + c x
                            |      ---------------- dx
                            |              x
                           /
                           0

--> First, grab the integrand.
--> Take the part inside the `log', put under a common denominator and expand:

> 1 - b*x/(a+c*x^2):
> normal("):
> arg:=expand(log(")/x);
                                      2              2
                            ln(a + c x )   ln(a + c x  - b x)
                   arg := - ------------ + ------------------
                                  x                 x

--> Ok, we now have a canonical form for the integrand.
--> Next we find the roots for both polynomials inside the logs:

> res1:=solve(a+c*x^2,x);
                                       1/2           1/2
                                (- c a)       (- c a)
                        res1 := ----------, - ----------
                                     c             c

--> here, we 2 roots of equal magnitude and opposite sign
--> and the other polynomial has 2 distinct roots:

> res2:=solve(a+c*x^2-b*x,x):

--> So `arg' can be put into the form:

> arg2:= -(ln(c^2*(x-r1)*(x+r1))/x) + (ln(c^2*(x-r2)*(x-r3))/x);
                       2                          2
                   ln(c  (x - r1) (x + r1))   ln(c  (x - r2) (x - r3))
         arg2 := - ------------------------ + ------------------------
                               x                          x

> arg2:=expand(arg2);
                    ln(x - r1)   ln(x + r1)   ln(x - r2)   ln(x - r3)
          arg2 := - ---------- - ---------- + ---------- + ----------
                         x            x            x            x

--> and the terms in c^2 have cancelled out.
--> Next, we get the indefinite integral term by term:

> ans:=map(int,arg2,x);
                   x                     x                 x           x + r1
  ans := - dilog(----) - ln(x - r1) ln(----) - ln(r1) ln(----) + dilog(------)
                  r1                    r1                r1             r1

                 x                     x             x                     x
       + dilog(----) + ln(x - r2) ln(----) + dilog(----) + ln(x - r3) ln(----)
                r2                    r2            r3                    r3

--> Next, plug in the roots:

> ans:=subs(r1=res1[1],r2=res2[1],r3=res2[2],ans);

                                              1/2
                      x c              (- c a)            x c
   ans := - dilog(----------) - ln(x - ----------) ln(----------)
                         1/2                c                1/2
                  (- c a)                             (- c a)

                                                /           1/2\
                                                |    (- c a)   |
                    1/2                         |x + ----------| c
             (- c a)            x c             \         c    /
        - ln(----------) ln(----------) + dilog(------------------)
                  c                1/2                     1/2
                            (- c a)                 (- c a)

                    x c                b + %1         x c               x c
        + dilog(2 ------) + ln(x - 1/2 ------) ln(2 ------) + dilog(2 ------)
                  b + %1                  c         b + %1            b - %1

                     b - %1         x c
        + ln(x - 1/2 ------) ln(2 ------)
                        c         b - %1

                                 2         1/2
%1 :=                          (b  - 4 c a)

--> Note that the above is a solution to the INDEFINITE integral.

--> Next, plug in the endpoints using limits:
--> First, limit at x->0 from the right:

> l1:=op(1,series(series(ans,x),x,1));

                                  1/2
                    2      (- c a)
        l1 := 1/6 Pi  - ln(----------) (ln(x) + 1/2 ln(c) - 1/2 ln(- a))
                                c

                        b - %1
             + ln(- 1/2 ------) (ln(2) + ln(x) + ln(c) - ln(b - %1))
                           c

                           1/2
                    (- c a)
             - ln(- ----------) (ln(x) + 1/2 ln(c) - 1/2 ln(- a))
                         c

                        b + %1
             + ln(- 1/2 ------) (ln(2) + ln(x) + ln(c) - ln(b + %1))
                           c

                                 2         1/2
%1 :=                          (b  - 4 c a)


--> That term has ln(x) behavior which is cause for concern
--> but computing the coefficient for the ln(x) term in the form ln(  ):

> culprit:=combine(coeff(collect(l1,ln(x)),ln(x),1),ln);
                                  2         1/2         2         1/2
                           (b - (b  - 4 c a)   ) (b + (b  - 4 c a)   )
         culprit := ln(1/4 -------------------------------------------)
                                               c a

--> we see that:

> expand(op(1,culprit));
                                       1

--> and ln(1)=0, so the coefficient for the ln(x) term is zero.

--> Good that means the limit at x->0 is finite

> l1:=subs(ln(x)=0,l1):

--> Next, limit for `ans' as x-> infinity

> l2:=subs(O=0,convert(asympt(asympt(ans,x),x,1),polynom)):
> l2:=expand(l2);

                2            2                                  2        2
  l2 := - 1/3 Pi  - 3/4 ln(c)  - 1/2 ln(c) ln(- a) + 1/4 ln(- a)  - ln(2)

                                                                             2
       - 2 ln(2) ln(c) + ln(2) ln(b + %1) + ln(c) ln(b + %1) - 1/2 ln(b + %1)

                                                             2
       + ln(2) ln(b - %1) + ln(c) ln(b - %1) - 1/2 ln(b - %1)

                                 2         1/2
%1 :=                          (b  - 4 c a)

--> Good so that result is finite.
--> So our solution is:

> sol:=expand(l2-l1);
                 2                   1/2                          1/2
  sol := - 1/2 Pi  + 1/2 ln(- (- c a)   ) ln(c) - 1/2 ln(- (- c a)   ) ln(- a)

       - ln(- b - %1) ln(2) - ln(- b - %1) ln(c) + ln(- b - %1) ln(b + %1)

       - ln(- b + %1) ln(2) - ln(- b + %1) ln(c) + ln(- b + %1) ln(b - %1)

                  2                            2
       + 1/2 ln(c)  + 1/2 ln(c) ln(- a) + ln(2)  + 2 ln(2) ln(c)

                       2                 2
       - 1/2 ln(b + %1)  - 1/2 ln(b - %1)

                                 2         1/2
%1 :=                          (b  - 4 c a)

--> Below is a number of manipulations to whittle down `sol':

> sol:=map(simplify,sol);
> sol:=expand(sol):
> sol:=subs(I=0,sol):
> collect(",ln(c)):
> map(combine,",ln):
> sol:=map(normal,",expanded):
> collect(",ln(2)):
> map(combine,",ln):
> sol:=map(normal,",expanded);

            2                                                      2
sol := ln(2)  - ln(2) ln(4) - ln(2) ln(c) - ln(2) ln(a) - 1/4 ln(c)

                               2
     - 1/2 ln(c) ln(a) - 3/4 Pi  + ln(- b + %1) ln(b - %1)

                                          2                 2                 2
     + ln(- b - %1) ln(b + %1) - 1/4 ln(a)  - 1/2 ln(b + %1)  - 1/2 ln(b - %1)

                                 2         1/2
%1 :=                          (b  - 4 c a)

--> So this is our solution!  Now if we want to plug in values like
--> say a=b=c=1 then we have to consider the case b^2-4*a*c < 0
--> so let's change the form of our solution:

> subs(sqrt(b^2-4*c*a)=I*sqrt(4*c*a-b^2),sol):
> evalc("):
> sol:=subs(I=0,"):

--> Finally our the solution:

> eq:= Int(log(1 - b*x/(a+c*x^2))/x, x = 0..infinity)=sol;

                            b x
         infinity ln(1 - --------)
            /                   2
           |             a + c x             2
  eq :=    |      ---------------- dx = ln(2)  - ln(2) ln(4) - ln(2) ln(c)
           |              x
          /
          0

                                2                           2                2
       - ln(2) ln(a) - 1/4 ln(c)  - 1/2 ln(c) ln(a) - 3/4 Pi  + 1/4 ln(4 c a)

       - arctan(%1, - b) arctan(- %1, b) - arctan(- %1, - b) arctan(%1, b)

                  2                    2                      2
       - 1/4 ln(a)  + 1/2 arctan(%1, b)  + 1/2 arctan(- %1, b)

                                         2 1/2
%1 :=                          (4 c a - b )

--> And the solution is in closed form!
# The above solution can be simplified to
> answer := -arcsin(1/2*b/(a*c)^(1/2))*(arcsin(1/2*b/(a*c)^(1/2))+Pi):

   Tony Scott  and  Greg Fee