On Fri, Apr 11, 2014 at 1:41 PM, Greg Marks <[email protected]> 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)
Why not
integrate(sqrt(1+9*x^4), x, 0, 2)
?
If you look up "arc length" in wikipedia (my preference over a
calculus textbook:-)
you get
integral(sqrt(1+(diff(y,x)^2),x,a,b)
for the arc length of the curve (x,y) from x=a to x=b.
>> 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.
--
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.
