home *** CD-ROM | disk | FTP | other *** search
- Newsgroups: sci.math.symbolic
- Path: sparky!uunet!utcsri!torn!watserv2.uwaterloo.ca!watdragon.uwaterloo.ca!daisy.uwaterloo.ca!gjfee
- From: gjfee@daisy.uwaterloo.ca (Greg Fee)
- Subject: Re: Log integral
- Message-ID: <Bx8BMI.4CC@watdragon.uwaterloo.ca>
- Sender: news@watdragon.uwaterloo.ca (USENET News System)
- Organization: University of Waterloo
- Date: Thu, 5 Nov 1992 06:06:17 GMT
- Lines: 260
-
-
- > 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
-