Re: Singularities at end point in integrations...
- To: mathgroup at smc.vnet.net
- Subject: [mg64161] Re: Singularities at end point in integrations...
- From: Paul Abbott <paul at physics.uwa.edu.au>
- Date: Fri, 3 Feb 2006 01:03:48 -0500 (EST)
- Organization: The University of Western Australia
- References: <drq0q6$muq$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
In article <drq0q6$muq$1 at smc.vnet.net>, ashesh <ashesh.cb at gmail.com>
wrote:
> 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
interest.
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]
Comparing
defint[1,2][3,4][2-10^(-20`30)]
(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.
Cheers,
Paul
_______________________________________________________________________
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)
AUSTRALIA http://physics.uwa.edu.au/~paul