numerical integration comparisons (GSl vs. mpmath vs. mpmath_rel)

gistfile1.py

test_functions = [ (e**(-x**2)*log(x), 17, 42),
                   (sqrt(1-x**2), 0, 1),
                   (sin(sin(x)), 0, 1),
                   (max_symbolic(sin(x), cos(x)), 0, pi),
                   (sin(x)*exp(cos(x)), 0, pi),
                   (exp(-x**100), 0, 1.1) ]


########################
# Burcin

import sage.libs.mpmath.all as mpmath
from mpmath.calculus.quadrature import TanhSinh, GaussLegendre

class RelativeErrorMixin(object):
    def estimate_error(self, results, prec, epsilon):
        mag = abs(results[-1])
        if len(results) == 2:
            return abs(results[0]-results[1]) / mag
        try:
            if results[-1] == results[-2] == results[-3]:
                return self.ctx.zero
            D1 = self.ctx.log(abs(results[-1]-results[-2])/mag, 10)
            D2 = self.ctx.log(abs(results[-1]-results[-3])/mag, 10)
            if D2 == 0:
                D2 = self.ctx.inf
            else:
                D2 = self.ctx.one / D2
        except ValueError:
            return epsilon
        D3 = -prec
        D4 = min(0, max(D1**2 * D2, 2*D1, D3))
        return self.ctx.mpf(10) ** int(D4)

class TanhSinhRel(RelativeErrorMixin, TanhSinh):
    pass

class GaussLegendreRel(RelativeErrorMixin, GaussLegendre):
    pass

########################
# Test code

def num_int_test_v2(f, a, b):
    LJ  = 25
    LJr = 35

    # GSL
    t = walltime()
    r = numerical_integral(f, a, b)[0]
    t = walltime(t)
    print "GSL:".ljust(LJ), str(r).ljust(LJr), " time: ", t

    # mpmath
    t = walltime()
    mp_f = lambda z: f(x = mpmath.mpmath_to_sage(z, 53))
    mpmath.mp.prec = 53
    r = mpmath.call(mpmath.quad, mp_f, [a, b])
    t = walltime(t)
    print "mpmath:".ljust(LJ), str(r).ljust(LJr), " time: ", t
    
    # quad to various precisions 
    for p in [64, 100]:
        mp_f = lambda z: f(x = mpmath.mpmath_to_sage(z, p))
        mpmath.mp.prec = p
        t = walltime()
        r = mpmath.call(mpmath.quad, mp_f, [a, b], method=GaussLegendreRel)
        t = walltime(t)
        s = "mpmath_rel (prec=%d):" % p
        print s.ljust(LJ), str(r).ljust(LJr), " time: ", t
        
def num_int_test_all():
    for f,a,b in test_functions:
        if callable(f):
            print "function: %s  on  [%s, %s]" % (f(x), a, b)
        else:
            print "function: %s  on  [%s, %s]" % (f, a, b)
        num_int_test_v2(f,a,b)
        print ""