Raymond Toy <toy.raymond <at> gmail.com> writes:
>
> On 12/6/12 9:07 PM, Vickram Ratnam wrote:
> > Hi,
> >
> >
> >
> > I?m new to Maxima. I am using Maxima for windows!
> >
> >
> >
> > I?m trying to compute a definite integral using Woollett?s tanh-sinh
> > maxima code and the maxima built-in integration methods but with no success!
> >
> >
> >
> > My function is the following:
> >
> >
> >
> > 1/sqrt((.0000812*(1+x)^4)+(.2719188*(1+x)^3)+(0*(1+x)^2)+.728)
>
> Just for fun, I thought I'd try to get an "exact" result for this
> elliptic integral.
>
> First, get a copy of ellint.mac from
> https://github.com/rtoy/maxima-hacks/tree/master/specfun and load it.
>
> Let's start with the quartic:
>
> y2 : (.0000812*(1+x)^4)+(.2719188*(1+x)^3)+(0*(1+x)^2)+.728;
>
> We want to compute integrate(1/sqrt(y2), x, 0, 1e13)
>
> Factor it into two quadratics:
>
> qf: expand(float(quartic_factor(y2,x))) ->
> [8.120000000000002e-5*x^2+.2721939664867234*x+.6497400204178281,
> 1.0*x^2+.6112501635044124*x+1.539077121112826]
>
> Try to convert these quadratics into the form a*x^2+b:
>
> subs:expand(float(ell_subst(qf[1],qf[2],x,s) ->
>
> [s = x/(x-.01670058983603612)+4.790748539009956/(x-.01670058983603612),
> .001898096671116776,1.0*s^2-.9971317440669707,1.0*s^2+13.91488129120323,
> -.2080105214217932]
>
> This means that the substitution s = x/(x-0.016)+4.9/(x-0.016) will
> convert integrate(1/sqrt(y2),x) to
>
> integrate(subs[5]/sqrt(subs[2])/sqrt(subs[3]*subs[4]), s)
>
> Or
>
> integrate(-4.774480664563152/(sqrt(1.0*s^2-.9971317440669707)
> *sqrt(1.0*s^2+13.91488129120323)), s)
>
> Let's convert that to elliptic functions:
>
> int: expand(float(-4.774480664563152*ellint1(subs[3],subs[4],s))) ->
>
> [- 1.236397167842715 inverse_jacobi_nc(1.001437220443617 s,
> .9331323181043008), 1.236397167842715
> inverse_jacobi_ds(.2589595088360969 s,
> .9331323181043008)]
>
> Two equivalent results are returned. (The routine doesn't know which is
> "better", so it returns both.) Let's use first one. This means
>
> integrate(1/sqrt(y2),x) =
> -1.236397167842715*inverse_jacobi_nc(1.001437220443617*s,.9331323181043008)
>
> where s = s =
> x/(x-.01670058983603612)+4.790748539009956/(x-.01670058983603612).
>
> The original limit of x from 1 to 1e13 become s from -286.86 to 1. So
> the final value is 3.373846103759024. (Whew! I wasn't sure I'd get
> anywhere near the answer you gave in different message!)
>
> You can use bigfloats to do the final evaluation, or maybe use bigfloats
> throughout the calculation above.
>
> But still, unless the original coefficients are truly exact values, I
> don't see the point of computing the integral to 300+ digits. Seems
> like it would be just random noise.
>
> Also, I think a little bit of analysis as done here goes a long way to
> simplify the problem. You now have some insight into what the function
> looks like for various limits of integration.
>
> Ray
>
I am trying to do so for two main reasons:
1) This definite integral is a computation for the Hubble parameter. I use the
Hubble parameter in a python code to calculate cosmological distances from a
given redshift. In order to calculate the conditions when the universe was
several seconds old, the redshift must be in the order of 1e+13 and above, up to
1e+43! However, even using such large redshift the conditions become negative or
meaningless, if higher precision is not ordered. For me to get a meaningful or
positive figure for such large redshifts, I need a minimum of 200 digits
accuracy.
2) My main program to perform these calculations is python with mpmath. Python
mpmath has always been very accurate and allow for the computation of larger
digits accuracy. However, I?m not so sure with regard to its quadrature module.
And so I checked the answer against mathematica (online version), and the result
diverged at about the 186th digit. So I needed a third source to know which one
is correct!!!
However yesterday i downloaded mathematica fully-functional trial version and
its result concur with python mpmath. The error lies with mathematica online
version which uses a lower recursive degree.
The coefficients in the hubble parameter which are .0000812, .2719188, 0 and
.728 are for the omega radiation, omega matter, omega curvature and omega lambda
respectively. The omega curvature in the present scenario is 0 for a flat
universe. Otherwise these figures can change depending on the models of the
universe being used.