home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #26 / NN_1992_26.iso / spool / sci / math / symbolic / 2910 next >
Encoding:
Text File  |  1992-11-04  |  9.7 KB  |  271 lines

  1. Newsgroups: sci.math.symbolic
  2. Path: sparky!uunet!utcsri!torn!watserv2.uwaterloo.ca!watdragon.uwaterloo.ca!daisy.uwaterloo.ca!gjfee
  3. From: gjfee@daisy.uwaterloo.ca (Greg Fee)
  4. Subject: Re: Log integral
  5. Message-ID: <Bx8BMI.4CC@watdragon.uwaterloo.ca>
  6. Sender: news@watdragon.uwaterloo.ca (USENET News System)
  7. Organization: University of Waterloo
  8. Date: Thu, 5 Nov 1992 06:06:17 GMT
  9. Lines: 260
  10.  
  11.  
  12. >                Integrate[Log[(1 - (bx/(a+cx^2))]/x,{x,0,Infinity}]
  13.  
  14. In Maple form, this integral is:
  15.                                              b x
  16.                           infinity ln(1 - --------)
  17.                              /                   2
  18.                             |             a + c x
  19.                             |      ---------------- dx
  20.                             |              x
  21.                            /
  22.                            0
  23.  
  24. We can make Maple generate the solution.   In fact, we can even get Maple
  25. to generate the form of the indefinite integral as well.
  26.  
  27.  
  28.     |\^/|      MAPLE V
  29. ._|\|   |/|_.  Copyright (c) 1981-1990 by the University of Waterloo.
  30.  \  MAPLE  /   All rights reserved.  MAPLE is a registered trademark of
  31.  <____ ____>   Waterloo Maple Software.
  32.       |        Type ? for help.
  33.  
  34. --> we want to solve:
  35. > Int(log(1 - b*x/(a+c*x^2))/x, x = 0..infinity);
  36.                                              b x
  37.                           infinity ln(1 - --------)
  38.                              /                   2
  39.                             |             a + c x
  40.                             |      ---------------- dx
  41.                             |              x
  42.                            /
  43.                            0
  44.  
  45. --> First, grab the integrand.
  46. --> Take the part inside the `log', put under a common denominator and expand:
  47.  
  48. > 1 - b*x/(a+c*x^2):
  49. > normal("):
  50. > arg:=expand(log(")/x);
  51.                                       2              2
  52.                             ln(a + c x )   ln(a + c x  - b x)
  53.                    arg := - ------------ + ------------------
  54.                                   x                 x
  55.  
  56. --> Ok, we now have a canonical form for the integrand.
  57. --> Next we find the roots for both polynomials inside the logs:
  58.  
  59. > res1:=solve(a+c*x^2,x);
  60.                                        1/2           1/2
  61.                                 (- c a)       (- c a)
  62.                         res1 := ----------, - ----------
  63.                                      c             c
  64.  
  65. --> here, we 2 roots of equal magnitude and opposite sign
  66. --> and the other polynomial has 2 distinct roots:
  67.  
  68. > res2:=solve(a+c*x^2-b*x,x):
  69.  
  70. --> So `arg' can be put into the form:
  71.  
  72. > arg2:= -(ln(c^2*(x-r1)*(x+r1))/x) + (ln(c^2*(x-r2)*(x-r3))/x);
  73.                        2                          2
  74.                    ln(c  (x - r1) (x + r1))   ln(c  (x - r2) (x - r3))
  75.          arg2 := - ------------------------ + ------------------------
  76.                                x                          x
  77.  
  78. > arg2:=expand(arg2);
  79.                     ln(x - r1)   ln(x + r1)   ln(x - r2)   ln(x - r3)
  80.           arg2 := - ---------- - ---------- + ---------- + ----------
  81.                          x            x            x            x
  82.  
  83. --> and the terms in c^2 have cancelled out.
  84. --> Next, we get the indefinite integral term by term:
  85.  
  86. > ans:=map(int,arg2,x);
  87.                    x                     x                 x           x + r1
  88.   ans := - dilog(----) - ln(x - r1) ln(----) - ln(r1) ln(----) + dilog(------)
  89.                   r1                    r1                r1             r1
  90.  
  91.                  x                     x             x                     x
  92.        + dilog(----) + ln(x - r2) ln(----) + dilog(----) + ln(x - r3) ln(----)
  93.                 r2                    r2            r3                    r3
  94.  
  95. --> Next, plug in the roots:
  96.  
  97. > ans:=subs(r1=res1[1],r2=res2[1],r3=res2[2],ans);
  98.  
  99.                                               1/2
  100.                       x c              (- c a)            x c
  101.    ans := - dilog(----------) - ln(x - ----------) ln(----------)
  102.                          1/2                c                1/2
  103.                   (- c a)                             (- c a)
  104.  
  105.                                                 /           1/2\
  106.                                                 |    (- c a)   |
  107.                     1/2                         |x + ----------| c
  108.              (- c a)            x c             \         c    /
  109.         - ln(----------) ln(----------) + dilog(------------------)
  110.                   c                1/2                     1/2
  111.                             (- c a)                 (- c a)
  112.  
  113.                     x c                b + %1         x c               x c
  114.         + dilog(2 ------) + ln(x - 1/2 ------) ln(2 ------) + dilog(2 ------)
  115.                   b + %1                  c         b + %1            b - %1
  116.  
  117.                      b - %1         x c
  118.         + ln(x - 1/2 ------) ln(2 ------)
  119.                         c         b - %1
  120.  
  121.                                  2         1/2
  122. %1 :=                          (b  - 4 c a)
  123.  
  124. --> Note that the above is a solution to the INDEFINITE integral.
  125.  
  126. --> Next, plug in the endpoints using limits:
  127. --> First, limit at x->0 from the right:
  128.  
  129. > l1:=op(1,series(series(ans,x),x,1));
  130.  
  131.                                   1/2
  132.                     2      (- c a)
  133.         l1 := 1/6 Pi  - ln(----------) (ln(x) + 1/2 ln(c) - 1/2 ln(- a))
  134.                                 c
  135.  
  136.                         b - %1
  137.              + ln(- 1/2 ------) (ln(2) + ln(x) + ln(c) - ln(b - %1))
  138.                            c
  139.  
  140.                            1/2
  141.                     (- c a)
  142.              - ln(- ----------) (ln(x) + 1/2 ln(c) - 1/2 ln(- a))
  143.                          c
  144.  
  145.                         b + %1
  146.              + ln(- 1/2 ------) (ln(2) + ln(x) + ln(c) - ln(b + %1))
  147.                            c
  148.  
  149.                                  2         1/2
  150. %1 :=                          (b  - 4 c a)
  151.  
  152.  
  153. --> That term has ln(x) behavior which is cause for concern
  154. --> but computing the coefficient for the ln(x) term in the form ln(  ):
  155.  
  156. > culprit:=combine(coeff(collect(l1,ln(x)),ln(x),1),ln);
  157.                                   2         1/2         2         1/2
  158.                            (b - (b  - 4 c a)   ) (b + (b  - 4 c a)   )
  159.          culprit := ln(1/4 -------------------------------------------)
  160.                                                c a
  161.  
  162. --> we see that:
  163.  
  164. > expand(op(1,culprit));
  165.                                        1
  166.  
  167. --> and ln(1)=0, so the coefficient for the ln(x) term is zero.
  168.  
  169. --> Good that means the limit at x->0 is finite
  170.  
  171. > l1:=subs(ln(x)=0,l1):
  172.  
  173. --> Next, limit for `ans' as x-> infinity
  174.  
  175. > l2:=subs(O=0,convert(asympt(asympt(ans,x),x,1),polynom)):
  176. > l2:=expand(l2);
  177.  
  178.                 2            2                                  2        2
  179.   l2 := - 1/3 Pi  - 3/4 ln(c)  - 1/2 ln(c) ln(- a) + 1/4 ln(- a)  - ln(2)
  180.  
  181.                                                                              2
  182.        - 2 ln(2) ln(c) + ln(2) ln(b + %1) + ln(c) ln(b + %1) - 1/2 ln(b + %1)
  183.  
  184.                                                              2
  185.        + ln(2) ln(b - %1) + ln(c) ln(b - %1) - 1/2 ln(b - %1)
  186.  
  187.                                  2         1/2
  188. %1 :=                          (b  - 4 c a)
  189.  
  190. --> Good so that result is finite.
  191. --> So our solution is:
  192.  
  193. > sol:=expand(l2-l1);
  194.                  2                   1/2                          1/2
  195.   sol := - 1/2 Pi  + 1/2 ln(- (- c a)   ) ln(c) - 1/2 ln(- (- c a)   ) ln(- a)
  196.  
  197.        - ln(- b - %1) ln(2) - ln(- b - %1) ln(c) + ln(- b - %1) ln(b + %1)
  198.  
  199.        - ln(- b + %1) ln(2) - ln(- b + %1) ln(c) + ln(- b + %1) ln(b - %1)
  200.  
  201.                   2                            2
  202.        + 1/2 ln(c)  + 1/2 ln(c) ln(- a) + ln(2)  + 2 ln(2) ln(c)
  203.  
  204.                        2                 2
  205.        - 1/2 ln(b + %1)  - 1/2 ln(b - %1)
  206.  
  207.                                  2         1/2
  208. %1 :=                          (b  - 4 c a)
  209.  
  210. --> Below is a number of manipulations to whittle down `sol':
  211.  
  212. > sol:=map(simplify,sol);
  213. > sol:=expand(sol):
  214. > sol:=subs(I=0,sol):
  215. > collect(",ln(c)):
  216. > map(combine,",ln):
  217. > sol:=map(normal,",expanded):
  218. > collect(",ln(2)):
  219. > map(combine,",ln):
  220. > sol:=map(normal,",expanded);
  221.  
  222.             2                                                      2
  223. sol := ln(2)  - ln(2) ln(4) - ln(2) ln(c) - ln(2) ln(a) - 1/4 ln(c)
  224.  
  225.                                2
  226.      - 1/2 ln(c) ln(a) - 3/4 Pi  + ln(- b + %1) ln(b - %1)
  227.  
  228.                                           2                 2                 2
  229.      + ln(- b - %1) ln(b + %1) - 1/4 ln(a)  - 1/2 ln(b + %1)  - 1/2 ln(b - %1)
  230.  
  231.                                  2         1/2
  232. %1 :=                          (b  - 4 c a)
  233.  
  234. --> So this is our solution!  Now if we want to plug in values like
  235. --> say a=b=c=1 then we have to consider the case b^2-4*a*c < 0
  236. --> so let's change the form of our solution:
  237.  
  238. > subs(sqrt(b^2-4*c*a)=I*sqrt(4*c*a-b^2),sol):
  239. > evalc("):
  240. > sol:=subs(I=0,"):
  241.  
  242. --> Finally our the solution:
  243.  
  244. > eq:= Int(log(1 - b*x/(a+c*x^2))/x, x = 0..infinity)=sol;
  245.  
  246.                             b x
  247.          infinity ln(1 - --------)
  248.             /                   2
  249.            |             a + c x             2
  250.   eq :=    |      ---------------- dx = ln(2)  - ln(2) ln(4) - ln(2) ln(c)
  251.            |              x
  252.           /
  253.           0
  254.  
  255.                                 2                           2                2
  256.        - ln(2) ln(a) - 1/4 ln(c)  - 1/2 ln(c) ln(a) - 3/4 Pi  + 1/4 ln(4 c a)
  257.  
  258.        - arctan(%1, - b) arctan(- %1, b) - arctan(- %1, - b) arctan(%1, b)
  259.  
  260.                   2                    2                      2
  261.        - 1/4 ln(a)  + 1/2 arctan(%1, b)  + 1/2 arctan(- %1, b)
  262.  
  263.                                          2 1/2
  264. %1 :=                          (4 c a - b )
  265.  
  266. --> And the solution is in closed form!
  267. # The above solution can be simplified to
  268. > answer := -arcsin(1/2*b/(a*c)^(1/2))*(arcsin(1/2*b/(a*c)^(1/2))+Pi):
  269.  
  270.    Tony Scott  and  Greg Fee
  271.