Dear sage-devel,
I am in the process of translating some lecture notes on "experimental
mathematics" from Mathematica to Sage, and I'm running into some
issues that I'll eventually ask about here after they've been brewing
in my head for some time.
Here's the first installment: suppose we want to compute a definite
integral with high accuracy. In Mathematica, one does something like
NIntegrate[1/Sqrt[1-x^4], {x, 0, 1}, {WorkingPrecision->200, AccuracyGoal->100}]
which returns
1.31102877714605990523241979494555970684137747571581158140841085190039529353520712511514776648071454672306787633589160902780447263249566096273945720229943750924300489083176406490206156616780035633786212
Apart from the fact that one needs to give both the working precision
and the desired accuracy, this is fairly clean and simple syntax.
(Oh, the answer has 200 digits but of course not all of them are
correct; however, more than 100 of them are correct so the accuracy
goal is attained.)
The situation in Sage is far from clean and simple (but maybe I'm
wrong and someone can show me the right way). I can do
sage: f = 1/sqrt(1-x^4)
sage: f.nintegral(x, 0, 1)
(1.3110287771460889, 1.411326611133745e-10, 315, 0)
Trying introspection to find out what the output values mean
gives me "Please see :obj:`sage.calculus.calculus.nintegral` for more
details", which I guess is fair enough. Reading that docstring tells
me that the function uses Maxima, and that it is possible to increase
the default precision, but only up to a point (in particular, not far
enough for the above example). There is a note at the end of the
docstring explaining that Pari can do what I want, and that what I
need to do in my case is this:
sage: gp.set_real_precision(200)
28
sage: gp.eval('intnum(x=0, 1, 1/sqrt(1-x^4))')
'1.3110287771460599052324197949455597068413774757158115814084108519003952935352071251151477664807145467230678729876212'
sage: RealField(332)(gp.eval('intnum(x=0, 1, 1/sqrt(1-x^4))'))
1.31102877714605990523241979494555970684137747571581158140841085190039529353520712511514776648071455
(Note that RealField(332) is 332 binary digits, hence about 100
decimal digits.)
I would like for this to be much more clueless-user-friendly. I'm not
sure what the best way to achieve that would be, though. One
possibility is to have something like
sage: f.nintegral(x, 0, 1, implementation='pari', precision=332)
(precision given in binary digits because that's how it works in other
situations in Sage.)
This would set the gp real precision to some initial guesstimate, run
the intnum command, and look at how many decimals are in the result.
If that's less than 100, increase the real precision and try again.
Note that the output of this would not be a quadruple as it is now (as
a wrapper to Maxima), so that's a bit of a problem.
There's also a function nintegrate which wraps GSL and returns a pair
(value, error). It would be nice to uniformize the whole thing
instead of having somewhat different behaviours depending on which of
the underlying implementation is used. So maybe the default simplest
form could be
sage: f.nintegral(x, 0, 1, precision=my_favorite_number)
and the implementation used is chosen based on what my_favorite_number
is. If it's small enough that GSL would give the right answer, use
GSL (very fast!) If it's bigger than that but still ok for Maxima
(i.e. small enough for quadpack), use Maxima. Otherwise use Pari and
fiddle with the real precision until you get enough digits.
Any other ideas? I think that numerical integration is a common
enough task in many circles, that it would do Sage some good to have a
simple and clean way of doing it.
Best,
Alex
--
Alex Ghitza -- Lecturer in Mathematics -- The University of Melbourne
-- Australia -- http://www.ms.unimelb.edu.au/~aghitza/
--
To post to this group, send an email to sage-devel@googlegroups.com
To unsubscribe from this group, send an email to
sage-devel+unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sage-devel
URL: http://www.sagemath.org
