On 11 April 2014 20:19, John Cremona <john.crem...@gmail.com> wrote:
> Here is how I would implement it:
Apologies, wrong thread.
>
> sage: v = [385, 231, 165, 105]
> sage: _,_,U = Matrix(v).smith_form()
> sage: U.column(0)
> (-38, 76, -19, 2)
> sage: U.column(0) * vector(v)
> 1
>
>
>
> On 11 April 2014 19:42, kcrisman <kcris...@gmail.com> wrote:
>>
>>
>> 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 sage-support+unsubscr...@googlegroups.com.
>> To post to this group, send email to sage-support@googlegroups.com.
>> 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 sage-support+unsubscr...@googlegroups.com.
To post to this group, send email to sage-support@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.
