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

Re: integration using PSQL algorithm

• To: mathgroup at smc.vnet.net
• Subject: [mg52368] Re: integration using PSQL algorithm
• From: Paul Abbott <paul at physics.uwa.edu.au>
• Date: Thu, 25 Nov 2004 05:49:38 -0500 (EST)
• Organization: The University of Western Australia
• References: <co1dl9$sb8$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

In article <co1dl9$sb8$1 at smc.vnet.net>, Arturas Acus <acus at itpa.lt> 
wrote:

> I have some reasons to suspect that some integral may have closed form
> in terms of powers of Pi: a+b*Pi+c*Pi^2+d*Pi^3+e*Pi^4 (a,b,c,d,e being
> algebraic numbers, probably with a=d=0). 

If the coefficients were positive integers less than the base Pi you 
could just use RealDigits. For example,

  num = 1 + Pi + 2 Pi^2 + Pi^3 + 2 Pi^4

  coefs = RealDigits[N[num, 25], Pi]

  num == Expand[FromDigits[coefs, Pi]]

> Assuming I can calculate the
> integral numerically to desired precision, is it possible to verify this
> hypothesis. As far as I know, the PSLQ algorithm, which helps to find
> algebraic relations between real numbers, can be useful here. The only
> (Mathematica related) reference I know is
> http://library.wolfram.com/infocenter/MathSource/4263/
> , but the example here is too simple. 

In The Mathematica Journal 6(2): 29-30, you will find code for 
TranscendentalRecognize[] that rationalizes the given transcendental 
basis and then uses rational arithmetic and LatticeReduce to find the 
"simplest" (rational) representation for the floating point number in 
that basis:

TranscendentalRecognize[n_, basis_] := 
  Module[{c, d, digs, e, id, lat, powerten, r, s, vals},    
    {d, e} = RealDigits[n]; 
    s = Sign[n]; 
    c = FromDigits[d]; 
    powerten = 10^(Length[d] - e);     
    digs = (RealDigits[N[#1, -e + Length[d] + 5]] & ) /@ basis;
    r = (FromDigits[Take[First[#1], -e + Last[#1] + 
      Length[d]]] & ) /@ digs;
   lat = Transpose[Append[IdentityMatrix[Length[basis] + 2],
      Flatten[{powerten, r, c}]]];
   vals = Take[First[LatticeReduce[lat]], Length[basis] + 2];     
   Expand[-((s*(Take[vals, {2, -2}] . basis + First[vals]))/Last[vals])]
 ]

This code works fine on the example above

  TranscendentalRecognize[N[num], Pi^Range[4]]

but it also works with more general examples such as

  num = 1/4 - (2/3) Pi + (40/11) Pi^2 + Pi^4;

  TranscendentalRecognize[N[num, 30], Pi^Range[4]]

> As for details, the integral in consideration is
> 
> integral = (-6*(-1 + a)*(4*Sqrt[7 + a] + Sqrt[2]*(5 + a)*Pi))/((5 +
> a)*(7 + a)^(3/2)) - 
>  (24*(-1 + a)*(6 + a)*Log[6 + a - Sqrt[5 + a]*Sqrt[7 + a]])/((5 + a)*(7
> + a))^(3/2)
> 
> with a= Cos[4* phi], integrated from 0 to Pi/2 (or to 2 Pi, no matter)
> The problems arises from second term.

Mathematica can integrate the first term in closed form. The result 
involves EllipticE[4/3] and EllipticE[Pi/3, 4/3]. The first of these can 
be written as Hypergeometric2F1[-1/2, 1/2, 1, 4/3] times Pi/2 but there 
is no simple form for the second. Unless there is cancellation between 
the two terms, I don't see how algebraic coefficients can result.

Working to high precision,

  Clear[num]

  num[n_] := num[n] = NIntegrate[Evaluate[integral /. a -> Cos[4 phi]], 
    {phi, 0, Pi/2}, WorkingPrecision -> n]

I find that

  tr = TranscendentalRecognize[num[150], Pi^Range[4]]

yields

  -(402550287036105137464991/240420095219532969056072) - 
   (398516724123798207595259*Pi)/240420095219532969056072 - 
   (177542572056568175522549*Pi^2)/240420095219532969056072 -   
   (502586551534098939521059*Pi^3)/240420095219532969056072 + 
   (196649748041515569157439*Pi^4)/240420095219532969056072

and this agrees with the numerical value to better than 1 part in 
10^(142):

  num[150] - tr

The rational coefficients in this basis are

  coefs = CoefficientList[tr, Pi]

However, changing the precision causes the coefficients to vary in such 
a way that makes me quite confident your suspicion is wrong.

Alternatively, you can use

  << NumberTheory`Recognize`

(which itself uses LatticeReduce) to attempt to find simple algebraic 
numbers for each of the coefficients. The resulting low order 
polynomials have large coefficients indicating that, again, your 
suspicion is likely unfounded.

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)         
35 Stirling Highway
Crawley WA 6009                      mailto:paul at physics.uwa.edu.au 
AUSTRALIA                            http://physics.uwa.edu.au/~paul