Elliptic integrals



Michel Talon wrote:
> Richard Fateman wrote:
>
>   
>> Unfortunately it doesn't.  The "steps" indicated are not explained
>> algorithmically, and in fact there may be
>> extraordinary difficulty in making them algorithmic.  That is why the
>> programs to do them have probably not
>> been written in the last several decades.  
>>     
>
> I don't know if a "general" program has been written or not to do elliptic
> integrals, (i thought that Trager had something like that in Axiom, but i
> have just checked Axiom doesn't do elliptic integrals) but i have a
> tendency to think that special cases are much more useful. Anyways Maple
> has no problem to do:
>
> rose% maple
>     |\^/|     Maple 10 (X86 64 LINUX)
>
>   
>> lprint(int(1/sqrt(x^4+1),x));
>>     
> 1/(1/2*2^(1/2)+1/2*I*2^(1/2))*(1-I*x^2)^(1/2)*(1+I*x^2)^(1/2)/(x^4+1)^(1/2)*
> EllipticF(x*(1/2*2^(1/2)+1/2*I*2^(1/2)),I)
>   
Here is what my current elliptic integral code does
(http://common-lisp.net/~rtoy/ellint3.mac).  It's really incomplete
because you have to do a lot of stuff by hand, but it produces something
reasonable.  Unfortunately, it gets the sign wrong sometimes.

(%i2) load(ellint3);
Evaluation took 0.0200 seconds (0.0100 elapsed) using 283.867 KB.
(%o2) "/Users/toy/.maxima/ellint3.mac"
(%i3) quartic_factor(x^4+1,x);
Evaluation took 0.0500 seconds (0.0500 elapsed) using 186.031 KB.
(%o3) [x^2+sqrt(2)*x+1,x^2-sqrt(2)*x+1]

/* ellint1 expects two quadratic factors or a linear and quadratic factor
   and computes integrate(1/sqrt(s1*s2),x)
*/
(%i4) ellint1(%o3[1],%o3[2],x);
Evaluation took 0.0100 seconds (0.0000 elapsed) using 187.586 KB.
(%o4) [-sqrt(2)*?%inverse_jacobi_sc((x+1)/(sqrt(3-2*sqrt(2))*(x-1)),
                                    4*sqrt(2)/(2*sqrt(2)+3))
        /sqrt(2*sqrt(2)+3),
       sqrt(2)*?%inverse_jacobi_cs((x+1)/(sqrt(2*sqrt(2)+3)*(x-1)),
                                   4*sqrt(2)/(2*sqrt(2)+3))
        /sqrt(2*sqrt(2)+3)]

The current code is not smart enough to know which inverse function to
return so it returns both.  Don't know if %o4 is really right, but

(%i5) integrate(1/sqrt(1+t^4),t,0,inf);
Evaluation took 0.0100 seconds (0.0300 elapsed) using 331.008 KB.
(%o5) beta(1/4,1/4)/4

(%i10) bfloat(make_elliptic_f(limit(%o4,x,inf)-limit(%o4,x,0)));
Evaluation took 0.0400 seconds (0.0400 elapsed) using 2.822 MB.
(%o10) [-1.854074677301372b0,1.854074677301372b0]
(%i11) bfloat(%o5);
Evaluation took 0.0200 seconds (0.0200 elapsed) using 3.760 MB.
(%o11) 1.854074677301375b0

So it's not totally wrong. :-)  Also, I didn't show it, but the result
of make_elliptic_f, which converts the inverse functions into an
elliptic_f function is quite a bit different in appearance from what you
give.  I don't know if they're the same or not.

>> quit
>>     
> bytes used=2397208, alloc=2227816, time=0.06
>
> It does also the case x^3+1 but the result is obviously more complicated
> since exp(2i pi/3) appears.
>   
(%i13) ellint1(x+1,x^2+x+1,x);
Evaluation took 0.0000 seconds (0.0000 elapsed) using 58.430 KB.
(%o13) [-?%inverse_jacobi_nc((x+2)/x,3/4),
        ?%inverse_jacobi_ds((x+2)/(2*x),3/4)]

Don't know if that is right or not.

With a bit of work, we can even derive that elliptic_kc(sqrt(3)/4+1/2) =
3^(1/4)/4/sqrt(%pi)*gamma(1/6)*gamma(1/3) if we start with the elliptic
integral

integrate(1/sqrt(1-t^3),t,minf,1)

(Maxima can't do this integral, but if you break the integral into minf
to 0 and 0 to 1, Maxima can do both.)
>
> By the way if we multiply the abelian differential by a reasonable rational
> function, maple is able to express the result on elliptic integrals of
> first and third kind:
>
>   
>> lprint(int((x^3+5*x+7)/((x^2+6*x+8)*sqrt(x^4+1)),x));
>>     
I've made some attempt at integrals of the second kind, but nothing for
the third kind.  I didn't try it out and I'm not sure if it even works.


Ray