sympy/numerics/quad.py

   1 """
   2 Routines for arbitrary-precision numerical quadrature
   3 """
   4
   5 from sympy import oo
   6
   7 #----------------------------------------------------------------------------#
   8 #                        General classes / utilities                         #
   9 #----------------------------------------------------------------------------#
  10
  11 from float_ import Float
  12 from constants import *
  13 from functions import *
  14 from evalf_ import polyfunc
  15
  16 import math
  17
  18
  19 class Quadrature:
  20     """
  21     This class provides a standard interface for quadrature rules. A
  22     rule is applied through a function call (e.g. to integrate sin(x)
  23     from x=3 to x=8, SomeQuadRule(<options>)(lambda x: sin(x), 3, 8))
  24
  25     Each Quadrature subclass should define an '_eval' method that
  26     returns a tuple containing two values: (y, err) where y is the
  27     calculated value of the integral, and err is an estimate of
  28     the approximation error (or None if the error cannot be estimated
  29     from the rule).
  30
  31     By default, the quadrature rule is assumed to be defined on
  32     [-1, 1] or (-1, 1). __call__ automatically transforms other
  33     (finite) intervals to this default interval. This behavior can
  34     be modified by overriding the 'transform' method.
  35     """
  36
  37     def transform(self, f, a, b):
  38         """Given an integrand function f and boundaries a, b,
  39         return an equivalent integrand defined on [-1, 1]"""
  40         a = Float(a)
  41         b = Float(b)
  42         if (a, b) == (-1, 1):
  43             return f
  44         else:
  45             C = (b-a)/2
  46             D = (b+a)/2
  47             def g(x):
  48                 return C * f(D + C*x)
  49             return g
  50
  51     def __call__(self, f, a, b, extraprec=3, verbose=False):
  52         Float.store()
  53         Float.extraprec(extraprec)
  54         f = self.transform(f, a, b)
  55         try:
  56             s, err = self._eval(f, verbose)
  57         except Exception, e:
  58             Float.revert()
  59             raise e
  60         Float.revert()
  61         return +s, err
  62
  63     def adaptive(self, f, a, b, eps, steps=0, maxsteps=5000, verbose=False):
  64         """Apply rule adaptively (must support error estimation)"""
  65         s, err = self(f, a, b, verbose=verbose)
  66         if err <= eps or steps >= maxsteps:
  67             return s, err, steps+1
  68         if verbose:
  69             print steps, a, b
  70         mid = (a+b)/2
  71         s1, e1, steps = self.adaptive(f, a, mid, eps, steps+1, maxsteps, verbose)
  72         s2, e2, steps = self.adaptive(f, mid, b, eps, steps, maxsteps, verbose)
  73         return s1 + s2, e1 + e2, steps
  74
  75
  76 #----------------------------------------------------------------------------#
  77 #                       Trivial quadrature rules                             #
  78 #----------------------------------------------------------------------------#
  79
  80 # Mainly for demonstration purposes...
  81
  82 class Midpoint(Quadrature):
  83     def _eval(self, f, verbose=False):
  84         return 2*f(Float(0)), None
  85
  86 class Trapezoid(Quadrature):
  87     def _eval(self, f, verbose=False):
  88         return f(Float(-1)) + f(Float(1)), None
  89
  90
  91 #----------------------------------------------------------------------------#
  92 #                          Gaussian quadrature                               #
  93 #----------------------------------------------------------------------------#
  94
  95 def _iterfixpoint(h, x, eps):
  96     while 1:
  97         new = h(x)
  98         if x.ae(new, eps):
  99             return new
 100         x = new
 101
 102 def _gaussnode(pf, k, n):
 103     import math
 104     x = Float(-math.cos(math.pi*(k+1-0.25)/(n+0.5)))
 105     def h(r):
 106         p, q = pf(r)
 107         return r - p/q
 108     eps = Float((1, -Float._prec//2))
 109     x = _iterfixpoint(h, x, eps)
 110     w = 2 / (1-x**2) / (pf(x)[1])**2
 111     return x, w
 112
 113 _gausscache = {}
 114
 115
 116 class GaussLegendre(Quadrature):
 117     """
 118     Gauss-Legendre quadrature is highly efficient for smooth integrands
 119     on finite intervals.
 120     """
 121
 122     def __init__(self, n=None, verbose=False):
 123         self.n = n
 124         prec = dps = Float.getdps()
 125         self.prec = prec
 126
 127         # Reuse old nodes
 128         if n in _gausscache:
 129             cacheddps, xs, ws = _gausscache[n]
 130             if cacheddps >= dps:
 131                 self.x = [x for x in xs]
 132                 self.w = [w for w in ws]
 133                 return
 134
 135         if verbose:
 136             print ("calculating nodes for degree-%i Gauss-Legendre "
 137                 "quadrature..." % n)
 138
 139         Float.store()
 140         Float.setdps(2*prec + 5)
 141
 142         self.x = [None] * n
 143         self.w = [None] * n
 144
 145         from sympy.specfun import legendre
 146         pf = polyfunc(legendre(n, 'x'), True)
 147
 148         for k in xrange(n//2 + 1):
 149             if verbose and k % 4 == 0:
 150                 print "  node", k, "of", n//2
 151             x, w = _gaussnode(pf, k, n)
 152             self.x[k] = x
 153             self.x[n-k-1] = -x
 154             self.w[k] = self.w[n-k-1] = w
 155
 156         _gausscache[n] = (dps, self.x, self.w)
 157
 158         Float.revert()
 159
 160     def _eval(self, f, verbose=False):
 161         s = Float(0)
 162         for i in xrange(self.n):
 163             s += self.w[i] * f(self.x[i])
 164         return s, None
 165
 166
 167 class CompositeGaussLegendre(Quadrature):
 168     """
 169     Gauss-Legendre quadrature with error estimation. Two Gauss-Legendre
 170     rules of different degree are used, and the error is estimated as
 171     the difference between the two results.
 172     """
 173
 174     def __init__(self, deg1, deg2, verbose=False):
 175         self.rule1 = GaussLegendre(deg1, verbose)
 176         self.rule2 = GaussLegendre(deg2, verbose)
 177
 178     def _eval(self, f, verbose=False):
 179         s = self.rule1._eval(f, verbose)[0]
 180         t = self.rule2._eval(f, verbose)[0]
 181         return s, abs(s-t)
 182
 183
 184 #----------------------------------------------------------------------------#
 185 #                            Tanh-sinh quadrature                            #
 186 #----------------------------------------------------------------------------#
 187
 188 def _tsnode(k, h):
 189     # x[k] = tanh(pi/2 * sinh(k*h))
 190     # w[k] = pi/2 * cosh(k*h) / cosh(pi/2 sinh(k*h))**2
 191     t = Float(k) * h
 192     a = exp(t)
 193     pi = pi_float()
 194     ar = Float(1) / a
 195     sinht = (a - ar) / 2
 196     cosht = (a + ar) / 2
 197     b = exp((pi * sinht) / 2)
 198     br = Float(1) / b
 199     sinhb = (b - br) / 2
 200     coshb = (b + br) / 2
 201     return sinhb/coshb, (pi/2)*cosht/coshb**2
 202
 203 _tscache = {}
 204
 205 class TanhSinh(Quadrature):
 206     """
 207     The tanh-sinh quadrature rule (also known as the doubly-exponential
 208     quadrature rule) is well suited for integrals with singularities
 209     and/or extremely high precision levels (tens or hundreds of digits).
 210
 211     The level m is a small integer between, say, 5 for low precision
 212     levels and 10 for extremely high precision (it must be set higher
 213     if the integrand is ill-behaved).
 214
 215     This implementation of the tanh-sinh algorithm is based on the
 216     description given in Borwein, Bailey & Girgensohn, "Experimentation
 217     in Mathematics - Computational Paths to Discovery", A K Peters,
 218     2003, pages 312-313. It also described in various documents
 219     available online.
 220     """
 221
 222     def __init__(self, eps, m, verbose=False):
 223         # Compute abscissas and weights
 224
 225         prec = Float.getdps()
 226
 227         self.prec = prec
 228         self.eps = eps
 229         self.m = m
 230         self.h = h = Float(1) / 2**m
 231         self.x = []
 232         self.w = []
 233
 234         if (eps, m) in _tscache:
 235             self.x, self.w = _tscache[(eps, m)]
 236             return
 237
 238         if verbose:
 239             print ("calculating nodes for tanh-sinh quadrature with "
 240                 "epsilon %s and degree %i" % (eps, m))
 241
 242         Float.store()
 243         Float.setdps(prec + 10)
 244
 245         for k in xrange(20 * 2**m + 1):
 246             x, w = _tsnode(k, h)
 247             self.x.append(x)
 248             self.w.append(w)
 249             diff = abs(self.x[-1] - Float(1))
 250             if diff <= eps:
 251                 break
 252             if verbose and m > 6 and k % 300 == 299:
 253                 Float.store(); Float.setdps(5)
 254                 print "  progress", -log(diff, 10) / prec
 255                 Float.revert()
 256
 257         _tscache[(eps, m)] = self.x, self.w
 258
 259         Float.revert()
 260
 261     def _estimate_error(self, res):
 262         try:
 263             D1 = log(abs(res[-1]-res[-2]), 10)
 264             D2 = log(abs(res[-1]-res[-3]), 10)
 265         except ValueError:
 266             return self.eps
 267         D3 = -self.prec
 268         D4 = min(0, max(D1**2/D2, 2*D1, D3))
 269         return Float('0.1') ** -int(D4)
 270
 271     def _eval(self, f, verbose=False):
 272         S = Float(0)
 273         h = Float(1)
 274         m = self.m
 275         res = []
 276         for k in xrange(1, m+1):
 277             h = h / 2
 278             for i in xrange(0, len(self.x), 2**(m-k)):
 279                 if i % (2**(m-k+1)) != 0 or k == 1:
 280                     if i == 0:
 281                         S = S + self.w[0]*f(Float(0))
 282                     else:
 283                         S = S + (self.w[i])*(f(-self.x[i]) + f(self.x[i]))
 284             res.append(h*S)
 285             if k > 2:
 286                 err = self._estimate_error(res)
 287                 if err <= self.eps:
 288                     break
 289         return +res[-1], self._estimate_error(res)
 290
 291
 292 class AdaptiveTanhSinh(Quadrature):
 293     """
 294     Adaptive tanh-sinh quadrature algorithm.
 295     """
 296
 297     def __init__(self, initial):
 298         self.initial = initial
 299
 300     def adaptive(self, f, a, b, eps, steps=0, maxsteps=5000, verbose=False):
 301         prec = Float.getprec()
 302         for m in xrange(self.initial, 50):
 303             if verbose:
 304                 print "using tanh-sinh rule with m =", m
 305             ts = TanhSinh(eps, m, verbose=verbose)
 306             s, err = ts(f, a, b, verbose=verbose)
 307             steps = 2*len(ts.x)
 308             if err <= eps or steps >= maxsteps:
 309                 return s, err, steps
 310
 311
 312 #----------------------------------------------------------------------------#
 313 #                          User-friendly functions                           #
 314 #----------------------------------------------------------------------------#
 315
 316
 317 def nintegrate(f, a, b, method=0, maxsteps=5000, verbose=False):
 318     """
 319     Basic usage
 320     ===========
 321
 322     nintegrate(f, a, b) numerically evaluates the integral
 323
 324            - b
 325           |    f(x) dx.
 326          -   a
 327
 328
 329     The integrand f should be a callable that accepts a Float as input
 330     and outputs a Float or ComplexFloat; a and b should be Floats,
 331     numbers that can be converted to Floats, or +/- infinity.
 332
 333     A simple example:
 334
 335         >>> Float.setdps(15)
 336         >>> print nintegrate(lambda x: 2*x**2 - 4*x, 2, 3)
 337         2.66666666666667
 338
 339     Calculating the area of a unit circle, with 30 digits:
 340
 341         >>> Float.setdps(30)
 342         >>> print nintegrate(lambda x: 4*sqrt(1-x**2), 0, 1)
 343         3.14159265358979323846264338328
 344
 345     The integration interval can be infinite or semi-infinite:
 346
 347         >>> Float.setdps(15)
 348         >>> print nintegrate(lambda x: exp(-x)*sin(x), 0, oo)
 349         0.5
 350
 351     Integration methods and accuracy
 352     ================================
 353
 354     Nintegrate attempts to obtain a value that is fully accurate within
 355     the current working precision (i.e., correct to 15 decimal places
 356     at the default precision level). If nintegrate fails to reach full
 357     accuracy after a certain number of steps, it prints a warning
 358     message.
 359
 360     This message signifies either that the integral is either divergent
 361     or, if convergent, ill-behaved. It may still be possible to
 362     evaluate an ill-behaved integral by increasing the 'maxsteps'
 363     setting, changing the integration method, and/or manually
 364     transforming the integrand.
 365
 366     Nintegrate currently supports the following integration methods:
 367
 368         method = 0  :  Gaussian quadrature (default)
 369         method = 1  :  tanh-sinh quadrature
 370
 371     Gaussian quadrature is generally very efficient if the integration
 372     interval is finite and the integrand is smooth on the entire range
 373     (including the endpoints). It may fail if the integrand has many
 374     discontinuities, is highly oscillatory, or possesses integrable
 375     singularities.
 376
 377     The tanh-sinh algorithm is often better if the integration interval
 378     is infinite or if singularities are present at the endpoints;
 379     especially at very high precision levels. It does not perform well
 380     if there are singularities between the endpoints or the integrand
 381     is bumpy or oscillatory.
 382
 383     It may help to manually transform the integrand, e.g. changing
 384     variables to remove singularities or breaking up the integration
 385     interval so that singularities appear only at the endpoints.
 386
 387     The 'verbose' flag can be set to track the computation's progress.
 388     """
 389     dps = Float.getdps()
 390     Float.store()
 391     Float.setdps(dps + 3)
 392     prec = Float.getprec()
 393
 394     # Transform infinite or semi-infinite interval to a finite interval
 395     if a == -oo or b == oo:
 396         g = f
 397         if a == -oo and b == oo:
 398             def f(x):
 399                 # make adaptive quadrature work from the left
 400                 x = 1 - x
 401                 return g(x) + g(Float(1)/x)/x**2 + g(-x) + g(Float(-1)/x)/(-x)**2
 402         elif b == oo:
 403             aa = Float(a)
 404             def f(x):
 405                 x = 1 - x
 406                 return g(x + aa) + g(Float(1)/x + aa)/x**2
 407         elif a == -oo:
 408             bb = Float(b)
 409             def f(x):
 410                 x = 1 - x
 411                 return g(-x + bb) + g(Float(-1)/x + bb)/(-x)**2
 412         a, b = Float(0), Float(1)
 413     else:
 414         a, b = Float(a), Float(b)
 415
 416     eps = Float((1, -prec+4))
 417
 418     if method == 0:
 419         degree = int(5 + dps**0.8)
 420         rule = CompositeGaussLegendre(degree, degree//2, verbose)
 421     elif method == 1:
 422         rule = AdaptiveTanhSinh(initial = int(3 + max(0, math.log(prec/30.0, 2))))
 423
 424     else:
 425         Float.revert()
 426         raise ValueError("unknown method")
 427
 428     s, err, steps = rule.adaptive(f, Float(a), Float(b), eps, steps=0, maxsteps=maxsteps, verbose=verbose)
 429
 430     Float.revert()
 431
 432     if not err.ae(0):
 433         Float.store()
 434         Float.setdps(1)
 435         print "Warning: failed to reach full accuracy.", \
 436             "Estimated magnitude of error: ", str(err)
 437         Float.revert()
 438
 439     return +s