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