On Friday, April 11, 2014 1:41:55 PM UTC-4, Greg Marks wrote:
>
> 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?
>
This is a great post, thanks! The problem is that Maxima knows some
antiderivatives for such things, sympy knows others, and I don't know where
the intersection is.
n(1/sqrt(3)*integrate(sqrt(1+x^4), x, 0, 2*sqrt(3),algorithm='sympy'),
digits=100)
is another thing I would have gone with, but
AttributeError: 'gamma' object has no attribute '_sage_'
It would probably not be a good idea to have integration try every single
subsystem for every integral though.
There are a lot of tickets which would do exactly what you are talking
about. Gluing sympy (which is really much improved from its early days,
and very useful) in better, as we see didn't happen in my example, is one
way to help; http://trac.sagemath.org/ticket/2516 might be another...
see http://trac.sagemath.org/ticket/14446 for something similar. Aargh.
I thought one could get arbitrary numerical integration precision from the
GSL interface (numerical_integral) but I'm having trouble interpreting it
to do so.
Are you interested in helping us get these things organized better? It
would be greatly appreciated.
- kcrisman
--
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.
