%maxima contour_plot(w^2+y^2,[w,-4,4],[y,-4,4])$
%maxima plot2d(r^2,[r,-2,2])
["/projects/c866c77f-ee7b-4d60-a777-22d3492d10b1/maxout.gnuplot_pipes"]
%
reset()
#import the needed libraries import scipy
def ff(x): return x^2
a = vector([1.0,2.0,3.0,4.0,5.0]) a
(1.00000000000000, 2.00000000000000, 3.00000000000000, 4.00000000000000, 5.00000000000000)
c = srange(1,5,0.01)
scipy.integrate.quadpack(map(ff,c),c)
Error in lines 1-1
Traceback (most recent call last):
File "/usr/local/lib/python2.7/dist-packages/smc_sagews/sage_server.py", line 957, in execute
exec compile(block+'\n', '', 'single') in namespace, locals
File "", line 1, in <module>
TypeError: 'module' object is not callable
scipy.integrate.quad(ff,1,5)
(41.33333333333333, 4.588921835117313e-13)
scipy.integrate.simps(map(ff,c),c)
41.083833166665109
scipy.integrate.trapz(map(ff,c),c)
41.083899499998438
integral_numerical(ff,1,5)
(41.333333333333336, 4.588921835117314e-13)
db
File: /projects/sage/sage-dev/src/sage/gsl/integration.pyx Signature : integral_numerical(func, a, b=None, algorithm='qag', max_points=87, params=[], eps_abs=1e-06, eps_rel=1e-06, rule=6) Docstring : Returns the numerical integral of the function on the interval from a to b and an error bound. INPUT: * "a", "b" -- The interval of integration, specified as two numbers or as a tuple/list with the first element the lower bound and the second element the upper bound. Use "+Infinity" and "-Infinity" for plus or minus infinity. * "algorithm" -- valid choices are: * 'qag' -- for an adaptive integration * 'qags' -- for an adaptive integration with (integrable) singularities * 'qng' -- for a non-adaptive Gauss-Kronrod (samples at a maximum of 87pts) * "max_points" -- sets the maximum number of sample points * "params" -- used to pass parameters to your function * "eps_abs", "eps_rel" -- absolute and relative error tolerances * "rule" -- This controls the Gauss-Kronrod rule used in the adaptive integration: * rule=1 -- 15 point rule * rule=2 -- 21 point rule * rule=3 -- 31 point rule * rule=4 -- 41 point rule * rule=5 -- 51 point rule * rule=6 -- 61 point rule Higher key values are more accurate for smooth functions but lower key values deal better with discontinuities. OUTPUT: A tuple whose first component is the answer and whose second component is an error estimate. REMARK: There is also a method "nintegral" on symbolic expressions that implements numerical integration using Maxima. It is potentially very useful for symbolic expressions. EXAMPLES: To integrate the function x^2 from 0 to 1, we do sage: numerical_integral(x^2, 0, 1, max_points=100) (0.3333333333333333, 3.700743415417188e-15) To integrate the function sin(x)^3 + sin(x) we do sage: numerical_integral(sin(x)^3 + sin(x), 0, pi) (3.333333333333333, 3.700743415417188e-14) The input can be any callable: sage: numerical_integral(lambda x: sin(x)^3 + sin(x), 0, pi) (3.333333333333333, 3.700743415417188e-14) We check this with a symbolic integration: sage: (sin(x)^3+sin(x)).integral(x,0,pi) 10/3 If we want to change the error tolerances and gauss rule used: sage: f = x^2 sage: numerical_integral(f, 0, 1, max_points=200, eps_abs=1e-7, eps_rel=1e-7, rule=4) (0.3333333333333333, 3.700743415417188e-15) For a Python function with parameters: sage: f(x,a) = 1/(a+x^2) sage: [numerical_integral(f, 1, 2, max_points=100, params=[n]) for n in range(10)] # random output (architecture and os dependent) [(0.49999999999998657, 5.5511151231256336e-15), (0.32175055439664557, 3.5721487367706477e-15), (0.24030098317249229, 2.6678768435816325e-15), (0.19253082576711697, 2.1375215571674764e-15), (0.16087527719832367, 1.7860743683853337e-15), (0.13827545676349412, 1.5351659583939151e-15), (0.12129975935702741, 1.3466978571966261e-15), (0.10806674191683065, 1.1997818507228991e-15), (0.09745444625548845, 1.0819617008493815e-15), (0.088750683050217577, 9.8533051773561173e-16)] sage: y = var('y') sage: numerical_integral(x*y, 0, 1) Traceback (most recent call last): ... ValueError: The function to be integrated depends on 2 variables (x, y), and so cannot be integrated in one dimension. Please fix additional variables with the 'params' argument Note the parameters are always a tuple even if they have one component. It is possible to integrate on infinite intervals as well by using +Infinity or -Infinity in the interval argument. For example: sage: f = exp(-x) sage: numerical_integral(f, 0, +Infinity) # random output (0.99999999999957279, 1.8429811298996553e-07) Note the coercion to the real field RR, which prevents underflow: sage: f = exp(-x**2) sage: numerical_integral(f, -Infinity, +Infinity) # random output (1.7724538509060035, 3.4295192165889879e-08) One can integrate any real-valued callable function: sage: numerical_integral(lambda x: abs(zeta(x)), [1.1,1.5]) # random output (1.8488570602160455, 2.052643677492633e-14) We can also numerically integrate symbolic expressions using either this function (which uses GSL) or the native integration (which uses Maxima): sage: exp(-1/x).nintegral(x, 1, 2) # via maxima (0.50479221787318..., 5.60431942934407...e-15, 21, 0) sage: numerical_integral(exp(-1/x), 1, 2) (0.50479221787318..., 5.60431942934407...e-15) We can also integrate constant expressions: sage: numerical_integral(2, 1, 7) (12.0, 0.0) If the interval of integration is a point, then the result is always zero (this makes sense within the Lebesgue theory of integration), see https://trac.sagemath.org/12047: sage: numerical_integral(log, 0, 0) (0.0, 0.0) sage: numerical_integral(lambda x: sqrt(x), (-2.0, -2.0) ) (0.0, 0.0) In the presence of integrable singularity, the default adaptive method might fail and it is advised to use "'qags'": sage: b = 1.81759643554688 sage: F(x) = sqrt((-x + b)/((x - 1.0)*x)) sage: numerical_integral(F, 1, b) (inf, nan) sage: numerical_integral(F, 1, b, algorithm='qags') # abs tol 1e-10 (1.1817104238446596, 3.387268288079781e-07) AUTHORS: * Josh Kantor * William Stein * Robert Bradshaw * Jeroen Demeyer ALGORITHM: Uses calls to the GSL (GNU Scientific Library) C library.
def trapz1(y, x=None, dx=1.0, axis=-1): y = asanyarray(y) if x is None: d = dx else: x = asanyarray(x) if x.ndim == 1: d = diff(x) # reshape to correct shape shape = [1]*y.ndim shape[axis] = d.shape[0] d = d.reshape(shape) else: d = diff(x, axis=axis) nd = len(y.shape) slice1 = [slice(None)]*nd slice2 = [slice(None)]*nd slice1[axis] = slice(1, None) slice2[axis] = slice(None, -1) try: ret = (d * (y[slice1] + y[slice2]) / 2.0).sum(axis) except ValueError: # Operations didn't work, cast to ndarray d = np.asarray(d) y = np.asarray(y) ret = add.reduce(d * (y[slice1]+y[slice2])/2.0, axis) return ret