A problem that arose in a recent calculus class involved practical methods
for finding a numerical estimate of the arc length of the curve y = x^3
from (0,0) to (2,8). So, one might ask, could one use SAGE to perform
this computation, say, to 10000 decimal places? (By the way, I don't think
this is such an unreasonable request. On a computer with 8 GB of RAM, one
can easily compute the arc length of a circle to one billion decimal
places, for example, using Hanhong Xue's implementation of the Chudnovsky
algorithm in GMP.)
Here is how not to do it:
sage: %time n(1/sqrt(3)*integrate(sqrt(1+x^4), x, 0, 2*sqrt(3)),
> digits=10000)
> CPU times: user 1.53 s, sys: 0.07 s, total: 1.60 s
> Wall time: 1.60 s
> 8.630329222227694
>
Well, at least it didn't take very long. Pretend you didn't know how to
antidifferentiate the integrand using special functions. We could proceed
as follows:
sage: from sympy import integrate, sqrt, var
> sage: x=var('x')
> sage: f=sqrt(1+x**4)
> sage: f.integrate(x)
> x*gamma(1/4)*hyper((-1/2, 1/4), (5/4,),
> x**4*exp_polar(I*pi))/(4*gamma(5/4))
>
Performing the obvious simplifications in our heads, we now use this
antiderivative to evaluate our definite integral:
sage: from mpmath import mp, hyp2f1
> sage: mp.dps = 10001
> sage: %time arclength=2*hyp2f1(-1/2,1/4,5/4,-144)
> CPU times: user 144.20 s, sys: 0.20 s, total: 144.40 s
> Wall time: 144.55 s
> sage: place = file('/tmp/arclength.txt','w')
> sage: place.write(str(arclength))
> sage: place.close()
>
The file /tmp/arclength.txt now contains the value of the definite integral
to 10000 decimal places. (I haven't tested it, but I'm pretty sure the
alternatives
sage: mpmath.quadts(lambda x: 1/sqrt(3)*sqrt(1+x^4), (0, 2*sqrt(3)))
>
and
sage: new_prec=gp.set_precision(10000)
sage: gp('intnum(x = 0, 2*sqrt(3), 1/sqrt(3)*sqrt(1+x^4))')
>
would take much too long.)
The good news is that this problem can be solved in SAGE. There was a time
when one of the few areas in which Maple and Mathematica could claim
superiority to SAGE was in this sort of integration problem. The procedure
given above at least provides evidence, if not proof, that this is no
longer the case.
The bad news is that the commands I used were highly unintuitive! I had to
know to use sympy. I had to know which of its routines to import. I had
to know the syntax sympy would accept in defining the function (e.g. for
exponentiation). Then I had to know to call mpmath. I had to know the
name mpmath uses for the hypergeometric function, which is different from
the name sympy uses.
It would be great if all of this were seamlessly integrated together in
SAGE, so that if a newcomer simply types in the obvious thing (such as "how
not to do it" above), he or she gets the desired result. Short of that, is
there anything simpler than the somewhat tortuous path that I followed?
--
You received this message because you are subscribed to the Google Groups
"sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.
