Thierry Dumont wrote: > I have looked at the quadrature routines: > > -in gsl: it seems that qags routine is called: this is a sophisticated > procedure with step adaptation, and convergence acceleration with the > epsilon-algorithm. This should integrate some singular functions and > discontinuous functions. > -in scipy: things look not so sophisticated, but it's difficult to > understand what is really used! > > *All* this is from quadpack (you can find quadpack on netlib), > transcripted to C for the gsl. The gsl routine is certainly expensive, > but robust. This is more or less the most "universal" method (but a > universal method of integration cannot exist...), certainly too > powerfull in many cases... I would recommend to keep it. >
Scipy uses quadpack too. Below is an excerpt of the code that is called from the quad routine (http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html#scipy.integrate.quad), which is what I was planning on using for the scipy algorithm option. You can see this file at http://svn.scipy.org/svn/scipy/trunk/scipy/integrate/quadpack.py, and you can see they include quadpack at http://svn.scipy.org/svn/scipy/trunk/scipy/integrate/quadpack/. From the code below, depending on the options, they call qagse, qagie, qagpe, qawoe, qawfe, qawse, or qawce, all from quadpack. Does that influence your decision any? Note that with scipy, we call the fortran version of quadpack. Thanks, Jason def _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points): infbounds = 0 if (b != Inf and a != -Inf): pass # standard integration elif (b == Inf and a != -Inf): infbounds = 1 bound = a elif (b == Inf and a == -Inf): infbounds = 2 bound = 0 # ignored elif (b != Inf and a == -Inf): infbounds = -1 bound = b else: raise RuntimeError, "Infinity comparisons don't work for you." if points is None: if infbounds == 0: return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit) else: return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit) else: if infbounds !=0: raise ValueError, "Infinity inputs cannot be used with break points." else: nl = len(points) the_points = numpy.zeros((nl+2,), float) the_points[:nl] = points return _quadpack._qagpe(func,a,b,the_points,args,full_output,epsabs,epsrel,limit) def _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts): if weight not in ['cos','sin','alg','alg-loga','alg-logb','alg-log','cauchy']: raise ValueError, "%s not a recognized weighting function." % weight strdict = {'cos':1,'sin':2,'alg':1,'alg-loga':2,'alg-logb':3,'alg-log':4} if weight in ['cos','sin']: integr = strdict[weight] if (b != Inf and a != -Inf): # finite limits if wopts is None: # no precomputed chebyshev moments return _quadpack._qawoe(func,a,b,wvar,integr,args,full_output,epsabs,epsrel,limit,maxp1,1) else: # precomputed chebyshev moments momcom = wopts[0] chebcom = wopts[1] return _quadpack._qawoe(func,a,b,wvar,integr,args,full_output,epsabs,epsrel,limit,maxp1,2,momcom,chebcom) elif (b == Inf and a != -Inf): return _quadpack._qawfe(func,a,wvar,integr,args,full_output,epsabs,limlst,limit,maxp1) elif (b != Inf and a == -Inf): # remap function and interval if weight == 'cos': def thefunc(x,*myargs): y = -x func = myargs[0] myargs = (y,) + myargs[1:] return apply(func,myargs) else: def thefunc(x,*myargs): y = -x func = myargs[0] myargs = (y,) + myargs[1:] return -apply(func,myargs) args = (func,) + args return _quadpack._qawfe(thefunc,-b,wvar,integr,args,full_output,epsabs,limlst,limit,maxp1) else: raise ValueError, "Cannot integrate with this weight from -Inf to +Inf." else: if a in [-Inf,Inf] or b in [-Inf,Inf]: raise ValueError, "Cannot integrate with this weight over an infinite interval." if weight[:3] == 'alg': integr = strdict[weight] return _quadpack._qawse(func,a,b,wvar,integr,args,full_output,epsabs,epsrel,limit) else: # weight == 'cauchy' return _quadpack._qawce(func,a,b,wvar,args,full_output,epsabs,epsrel,limit) > t.d. --~--~---------~--~----~------------~-------~--~----~ 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 -~----------~----~----~----~------~----~------~--~---
