MathGroup Archive 2006

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Singularities at end point in integrations...

In article <drq0q6$muq$1 at>, ashesh <ashesh.cb at> 

>  I Need to perform an integration with poles and zeros in the integrand. 
>  Please let me know if there a way in Mathematican that can be used to 
>  perform the definite integral 
>     sqrt((x-a)*(x-b)/((x-c)*(x-d))) 
> between the limits (c,d), (a,d), (a,b) or (b,c).

Each definite integral can be expressed as a combination of elliptic 
integrals. Note that Mathematica can compute the indefinite integral:

  Integrate[Sqrt[(x-a) (x-b)/((x-c) (x-d))], x]

and from this you can obtain closed-form expressions for each case of 

From the symmetry of your integrand, there are really only three cases 
to consider: 

  (c,d) (both poles)
  (a,b) (both zeros)
  (a,d) of (b,c) (one zero and one pole)

To give a concrete example, suppose that we are interested in computing 
the integral over (a,b) and that d > c > b > x > a. Note that the 
definite integral will be pure imaginary since (x-a) (x-b) < 0, so we 
modify the integrand before computing the indefinite integral:

  Collect[Integrate[Sqrt[(x-a) (b-x)/((c-x) (d-x))], x], 
  {EllipticPi[a_, b_, c_], EllipticF[d_, e_], EllipticE[f_, g_]},   
    Simplify[#, d > c > b > x > a] & ]

Simplifying this expression,

  Simplify[%, d > c > b > x > a]

and then collecting terms again,

 intab[a_, b_][c_, d_][x_] = 
   Collect[%, {EllipticPi[a_,b_,c_], EllipticF[d_,e_], EllipticE[f_,g_]},
    Simplify[#, d > c > b > x > a] & ]

we end up with a moderately complicated exact expression for the 
indefinite integral. Now, since there are no zeros or poles for a < x < 
b, we can compute the definite integral over (a,x) (where x < b, because 
there is a singularity at x = b) by taking the difference of the 
indefinite integral evaluated at the limits a and x, still in terms of 
the symbolic parameters a, b, c, and d.

 defint[a_, b_][c_, d_][x_] = intab[a, b][c, d][x] - intab[a, b][c, d][a]

As a check, try {a -> 1, b -> 2, c -> 3, d -> 4}. The plot is reasonable.

  Plot[Evaluate@Re[defint[1, 2][3, 4][x]], {x, 0, 2}, PlotRange -> All]



(evaluated with high precision upper limit very near to but just less 
than 2) with

  NIntegrate[Sqrt[(x-1) (2-x)/((3-x) (4-x))], {x, 1, 2}]

shows excellent agreement.

> I have read about the routine in quadpack called "dqawse.f" which can perform 
> "integration of functions having algebraico-logarithmic end point 
> singularities".

No need to use quadpack here. NIntegrate can handle such integrals -- 
and many more besides.


Paul Abbott                                      Phone:  61 8 6488 2734
School of Physics, M013                            Fax: +61 8 6488 1014
The University of Western Australia         (CRICOS Provider No 00126G)    

  • Prev by Date: Re: Re: String Input and Variable Setting
  • Next by Date: Re: summing a series in mathematica
  • Previous by thread: Re: Singularities at end point in integrations...
  • Next by thread: Re: How to get range of InterpolatingFunction?