From: gjfee@daisy.uwaterloo.ca (Greg Fee)
Subject: definite integration
Organization: University of Waterloo
Date: Mon, 9 Nov 1992 16:26:08 GMT

Re: the definite integral

integrate(ln(1 - b*x/(a+c*x^2))/x, x=0..infinity)

Here is another derivation of the closed form based on differentiation
with respect to a parameter.

    |\^/|      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.
> f := ln(1 - b*x/(a+c*x^2))/x;
                                            b x
                                  ln(1 - --------)
                                                2
                                         a + c x
                             f := ----------------
                                          x

> M := a^(1/2)/c^(1/2);
                                         1/2
                                        a
                                   M := ----
                                         1/2
                                        c

# Invert integrand about M to transform [M,infinity] to [0,M] .
> f1 := -subs(x=M^2/x,f)*diff(M^2/x,x);
                                             b a
                                ln(1 - --------------)
                                           /      2 \
                                           |     a  |
                                       c x |a + ----|
                                           |       2|
                                           \    c x /
                          f1 := ----------------------
                                           x

> g := normal(f+f1);
                                              2
                                     - a - c x  + b x
                                ln(- ----------------)
                                                2
                                         a + c x
                         g := 2 ----------------------
                                           x

# New upper limit of integration is M, but transform it to 1
> g := normal(subs(x=M*x,g)*diff(M*x,x));
                                   1/2      2  1/2      1/2
                              - a c    - a x  c    + b a    x
                         ln(- -------------------------------)
                                       1/2         2
                                      c    a (1 + x )
                  g := 2 -------------------------------------
                                           x

> g1 := op(1,x/2*g);
                                 1/2      2  1/2      1/2
                            - a c    - a x  c    + b a    x
                    g1 := - -------------------------------
                                     1/2         2
                                    c    a (1 + x )

> g1 := convert(g1,parfrac,x);
                                            b x
                          g1 := 1 - ------------------
                                     1/2  1/2       2
                                    a    c    (1 + x )

# Let  b/a^(1/2)/c^(1/2) = 2*sin(k)
> ss := b=2*sin(k)*a^(1/2)*c^(1/2);
                                              1/2  1/2
                          ss := b = 2 sin(k) a    c

> h1 := subs(ss,g1);
                                          sin(k) x
                              h1 := 1 - 2 --------
                                                2
                                           1 + x

> h := 2/x*ln(h1);
                                           sin(k) x
                                  ln(1 - 2 --------)
                                                 2
                                            1 + x
                           h := 2 ------------------
                                           x

> hk := normal(diff(h,k));
                                          cos(k)
                         hk := - 4 -------------------
                                        2
                                   1 + x  - 2 sin(k) x

> hk := -4*cos(k)/(cos(k)^2+sin(k)^2+x^2-2*sin(k)*x);
                                          cos(k)
                 hk := - 4 -----------------------------------
                                 2         2    2
                           cos(k)  + sin(k)  + x  - 2 sin(k) x

# integrate with resect to x between 0 and 1
> ih := int(hk,x);
                                           2 x - 2 sin(k)
                      ih := - 4 arctan(1/2 --------------)
                                               cos(k)

> nih := normal(ih);
                                           x - sin(k)
                         nih := - 4 arctan(----------)
                                             cos(k)

> ihx := eval(subs(x=1,nih)-subs(x=0,nih));
                                  1 - sin(k)             sin(k)
                ihx := - 4 arctan(----------) - 4 arctan(------)
                                    cos(k)               cos(k)

> a1 := eval(subs(sin(k)/cos(k)=tan(k),ihx));
                                        1 - sin(k)
                       a1 := - 4 arctan(----------) - 4 k
                                          cos(k)

> a2 := subs(arctan((1-sin(k))/cos(k))=Pi/4-k/2,a1); 
                                a2 := - Pi - 2 k

# Now integrate k back
> a3 := int(a2,k);
                                               2
                               a3 := - Pi k - k

# Find the missing constant
> int(eval(subs(k=0,h)),x=0..1);
                                       0

> subs(k=0,a3);
                                       0

# Missing constant in a3 is 0 so a3 is correct
# Put back original parameters
> sk := solve(ss,{k});
                                                 b
                       sk := {k = arcsin(1/2 ---------)}
                                              1/2  1/2
                                             a    c

> answer := subs(sk,a3);
                                       b                       b     2
         answer := - Pi arcsin(1/2 ---------) - arcsin(1/2 ---------)
                                    1/2  1/2                1/2  1/2
                                   a    c                  a    c

> answer := factor(answer);
                                    b      /                    b     \
         answer := - arcsin(1/2 ---------) |Pi + arcsin(1/2 ---------)|
                                 1/2  1/2  |                 1/2  1/2 |
                                a    c     \                a    c    /