Here is how I would implement it:
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.
