#8321: numerical integration with arbitrary precision
-------------------------+--------------------------------------------------
Reporter: burcin | Owner: burcin
Type: defect | Status: needs_work
Priority: major | Milestone: sage-4.7.2
Component: symbolics | Keywords: numerics,integration
Work_issues: | Upstream: N/A
Reviewer: | Author: Stefan Reiterer
Merged: | Dependencies:
-------------------------+--------------------------------------------------
Comment(by benjaminfjones):
I wrote a small function to compare errors and timings of the GSL
implementation (of `numerical_integral`) and the mpmath implementation
provided by the patch. Here is the testing code:
{{{
#!python
def num_int_test(f, a, b):
""" Input: `f` should be a function of a single variable, [a, b] is a
domain of integration """
LJ = 25 # left justification
from sage.symbolic.integration.integral import definite_integral
exact_I = definite_integral(f, f.variables()[0], a, b)
print "Exact ".ljust(LJ) + " = %s" % exact_I
print "Exact .n()".ljust(LJ) + " = %s" % exact_I.n()
print "GSL".ljust(LJ) + " = %s" % numerical_integral(f, a,
b)[0]
print "mpmath".ljust(LJ) + " = %s" % definite_integral(f,
f.variables()[0], a, b, hold=True).n()
num_I = definite_integral(f, f.variables()[0], a, b, hold=True)
for p in [53, 64, 100, 200, 500, 1000]:
s = "mpmath (prec=%d)" % p
print s.ljust(LJ) + " = %s" % num_I.n(p)
print "Timings at prec=53:"
print " GSL: ",
timeit('numerical_integral(%s, %s, %s)' % (f, a, b))
print " mpmath: ",
timeit('definite_integral(%s, %s, %s, %s, hold=True).n()' % (f,
f.variables()[0], a, b))
}}}
Now I took 3 examples from pg. 132 of the PDF suggested above at
http://cannelle.lateralis.org/sagebook-1.0.pdf and tested them using
`num_int_test`:
{{{
# applied patch http://trac.sagemath.org/sage_trac/raw-
attachment/ticket/8321/trac_8321_rebase.patch
sage: x = var('x')
sage: from sage.symbolic.integration.integral import definite_integral
# f = e^(-x^2)*log(x) on [17, 42]
sage: num_int_test(e^(-x^2)*log(x), 17, 42)
Exact = integrate(e^(-x^2)*log(x), x, 17, 42)
Exact .n() = 2.59225286296247e-127
GSL = 2.5657285007e-127
mpmath = 2.59225286296247e-127
mpmath (prec=53) = 2.59225286296247e-127
mpmath (prec=64) = 2.59225286296226841e-127
mpmath (prec=100) = 2.5922528629622683570941971829e-127
mpmath (prec=200) =
2.5922528629622683570941971813123497913329093314031958452707e-127
mpmath (prec=500) =
2.56572850056105148291735639613047859001477095540203266250504462960653767523831161390072264481216550025343198689835357618681873598023026225266716067780e-127
mpmath (prec=1000) =
2.56572850056105148291735639613047859001477095540203266250504462960653767360416188079136395575326953119218247602307727367985551096000368640359367812179070686479198046287233104280204937504901221620134046153583613193738177820412122516350777255525035947116513676784199592200655526485894447669230515221776e-127
Timings at prec=53:
GSL: 625 loops, best of 3: 971 µs per loop
mpmath: 25 loops, best of 3: 26.9 ms per loop
# f = sqrt(1-x^2) on [0, 1]
sage: num_int_test(sqrt(1-x^2), 0, 1)
Exact = 1/4*pi
Exact .n() = 0.785398163397448
GSL = 0.785398167726
mpmath = 0.785398163397448
mpmath (prec=53) = 0.785398163397448
mpmath (prec=64) = 0.785398163397448310
mpmath (prec=100) = 0.78539816339744830961566084582
mpmath (prec=200) =
0.78539816339744830961566084581987572104929234984377645524374
mpmath (prec=500) =
0.785398163397448309615660845819875721049292349843776455243736148076954101571552249657008706335529266995537021628320576661773461152387645557931339852032
mpmath (prec=1000) =
0.785398163397448309615660845819875721049292349843776455243736148076954101571552249657008706335529266995537021628320576661773461152387645557931339852032120279362571025675484630276389911155737238732595491107202743916483361532118912058446695791317800477286412141730865087152613581662053348401815062285318
Timings at prec=53:
GSL: 625 loops, best of 3: 652 µs per loop
mpmath: 25 loops, best of 3: 12 ms per loop
# f = sin(sin(x)) on [0, 1]
sage: num_int_test(sin(sin(x)), 0, 1)
Exact = integrate(sin(sin(x)), x, 0, 1)
Exact .n() = 0.430606103120691
GSL = 0.430606103121
mpmath = 0.430606103120691
mpmath (prec=53) = 0.430606103120691
mpmath (prec=64) = 0.430606103120690605
mpmath (prec=100) = 0.43060610312069060491237735525
mpmath (prec=200) =
0.43060610312069060491237735524846578643360804182199746950463
mpmath (prec=500) =
0.430606103120690604912377355248465786433608041821997469504633350750892193636074792502000332212863495547968512886769444385260392350928954849458834511854
mpmath (prec=1000) =
0.430606103120690604912377355248465786433608041821997469504633350750892193636074792502000332212863495547968512886769444385260392350928954849458834511854394326788473583253436780737313870079328121429092122005425057044706514198162061316772646582265252772251628205725432156943890956907988745419355505731945
Timings at prec=53:
GSL: 625 loops, best of 3: 134 µs per loop
mpmath: 25 loops, best of 3: 11.2 ms per loop
}}}
In the second example, GSL's answer compared to the numerically evaluated
exact symbolic answer is significantly off. The mpmath answer at default
precision matches the numerically evaluated exact symbolic answer
perfectly.
As for the timings, as you can see GSL is *much* faster than the patch
code using mpmath. I tried using `fast_callable` in the case where the
domain precision is the same as RDF and that gives a big speedup, though
still not as fast as GSL:
{{{
#!python
def _evalf( ... ):
...
try:
precision = parent.prec()
mpmath.mp.prec = precision
if precision == RDF.precision():
mp_f = fast_callable(f, vars=[x], domain=RDF)
else:
mp_f = lambda z: f(x = mpmath.mpmath_to_sage(z,
precision))
except AttributeError:
mp_f = fast_callable(f, vars=[x], domain=RDF)
return mpmath.call(mpmath.quad, mp_f, [a,b])
}}}
New timings using `fast_callable`:
{{{
# Testing: added fast_callable to _evalf_
# results: definite_integral( ... ) is around 5x faster when called at RDF
precision
sage: timeit('numerical_integral(e^(-x^2)*log(x), 17, 42)') # base case
using GSL at float precision
625 loops, best of 3: 959 µs per loop
sage: timeit('definite_integral(e^(-x^2)*log(x), x, 17, 42,
hold=True).n()') # using fast_callable
125 loops, best of 3: 5.41 ms per loop
sage: timeit('definite_integral(e^(-x^2)*log(x), x, 17, 42,
hold=True).n(53)') # using fast_callable
125 loops, best of 3: 5.41 ms per loop
sage: timeit('definite_integral(e^(-x^2)*log(x), x, 17, 42,
hold=True).n(54)') # *not* using fast_callable
25 loops, best of 3: 28.1 ms per loop
sage: timeit('definite_integral(e^(-x^2)*log(x), x, 17, 42,
hold=True).n(64)') # *not* using fast_callable
25 loops, best of 3: 26.3 ms per loop
sage: timeit('definite_integral(e^(-x^2)*log(x), x, 17, 42,
hold=True).n(100)') # *not* using fast_callable
25 loops, best of 3: 32.2 ms per loop
}}}
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/8321#comment:31>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
