Ticket #7763: integration.patch

sage/calculus/all.py

@@ -1 +1 @@
 from calculus import maxima as maxima_calculus
 from calculus import (laplace, inverse_laplace,
-                      limit, lim)
+                      limit, lim, numerical_integral, nintegral, nintegrate)
 
 from functional import (diff, derivative,
                         expand,

sage/calculus/calculus.py

@@ -15 +15 @@
 
 - Tom Coates (2010-06-11): fixed Trac #9217
 
+- Gagan Sekhon (2010-01): Added numerical_integral to fix trac 7763
+
 The Sage calculus module is loosely based on the Sage Enhancement
 Proposal found at: http://www.sagemath.org:9001/CalculusSEP.
 
@@ -358 +360 @@
     sage: integrate(sqrt(cos(x)^2 + sin(x)^2), x, 0, 2*pi)
     2*pi
 """
+from  sage.symbolic.integration.external import maxima_integrator, sympy_integrator, mma_free_integrator
+from sage.gsl.integration import gsl_numerical_integral
+from mpmath import quad, mp
 
 import re
 from sage.rings.all import RR, Integer, CC, QQ, RealDoubleElement, algdep
@@ -366 +371 @@
 from sage.misc.latex import latex, latex_variable_name
 from sage.interfaces.maxima import Maxima
 from sage.misc.parser import Parser
+from sage.misc.functional import numerical_approx
 
 from sage.symbolic.ring import var, SR, is_SymbolicVariable
-from sage.symbolic.expression import Expression
+from sage.symbolic.expression import Expression, is_Expression
 from sage.symbolic.function import Function
 from sage.symbolic.function_factory import function_factory, function
 from sage.symbolic.integration.integral import indefinite_integral, \
@@ -559 +565 @@
     else:
         raise ValueError, "unknown algorithm: %s" % algorithm
 
-def nintegral(ex, x, a, b,
+#This is so the old numerical_integral will still work
+def numerical_integral(*arg, **kwds):
+    r"""
+    
+    Input:
+    
+    - arg - General arguments for an integral function such as (f,x,a,b). However, (f,a,b) is also valid, 
+    it just calls the old numerical_integral which is now called old gsl_numerical_integral. `f` can be 
+    symbolic expression such as `sin(x)` or not symbolic `sin`. 
+   
+   Note: There are other methods, but I am still trying to figure out how to include them
+        Possible algorithms are gp.intnum, Gi/Pynac
+   
+    - kwds - First specify algorithm specifies which algorithm is to be used. For example: alg="maxima".
+             Default algorithm "maxima" for symbolic expression and "mpmath" for non-symbolic expression.
+    
+            - If `f` is symbolic expression then possible algorithms are "maxima" (default), "sympy",
+             "mathematica","gsl".
+            - If `f` is not symbolic expression then possible algorithms is "mpmath".
+            
+    
+      - Depending on the algorithm selected, appropriate parameters
+                    -  "maxima" - This algorithm has additional parameters
+                          -  ``epsrel`` - desired_relative_error  (default: '1e-8') the
+                                desired relative error
+
+                          -  ``limit`` - maximum_num_subintervals (default: 200)
+                                maxima number of subintervals
+                                
+                    -  "sympy" - This algorithm has no additional parameters
+                    -  "mathematica" - This algorithm has no additional parameters
+                    -  "gsl" - This algorithm has additional parameters
+                          -- algorithm: valid choices are
+                            'qag' -- for an adaptive integration
+                            '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 pt rule
+                                    rule=2: 21 pt rule
+                                    rule=3: 31 pt rule
+                                    rule=4: 41 pt rule
+                                    rule=5: 51 pt rule
+                                    rule=6: 61 pt rule
+                            Higher key values are more accurate for smooth functions but lower
+                            key values deal better with discontinuities.
+                    -  "mpmath" - This algorithm has additional parameters
+                        -d: Decimal places
+                        
+     ALIAS: nintegrate and nintegral are the same as nintegral    
+     
+            
+    ::
+
+        sage: f(x) = x
+        sage: nintegral(f,x,0,1, epsrel='1e-14')
+        0.5
+        
+    EXAMPLES::
+
+        sage: f(x) = exp(-sqrt(x))
+        sage: f.nintegral(x, 0, 1)
+        (0.52848223531423055, 4.1633141378838452e-11, 231, 0)
+
+
+    We can also use the ``numerical_integral`` function,
+    which calls the GSL C library.
+
+    ::
+
+        sage: numerical_integral(f, 0, 1)
+        (0.52848223225314706, 6.83928460...e-07)
+
+    Note that in exotic cases where floating point evaluation of the
+    expression leads to the wrong value, then the output can be
+    completely wrong::
+
+        sage: f = exp(pi*sqrt(163)) - 262537412640768744
+
+    Despite appearance, `f` is really very close to 0, but one
+    gets a nonzero value since the definition of
+    ``float(f)`` is that it makes all constants inside the
+    expression floats, then evaluates each function and each arithmetic
+    operation using float arithmetic::
+
+        sage: float(f)
+        -480.0
+
+    Computing to higher precision we see the truth::
+
+        sage: f.n(200)
+        -7.4992740280181431112064614366622348652078895136533593355718e-13
+        sage: f.n(300)
+        -7.49927402801814311120646143662663009137292462589621789352095066181709095575681963967103004e-13
+
+    Now numerically integrating, we see why the answer is wrong::
+
+        sage: f.nintegrate(x,0,1)
+        (-480.00000000000011, 5.3290705182007538e-12, 21, 0)
+         
+    It is just because every floating point evaluation of return -480.0
+    in floating point.
+
+    Important note: using PARI/GP one can compute numerical integrals
+    to high precision::
+
+        sage: gp.eval('intnum(x=17,42,exp(-x^2)*log(x))')
+        '2.565728500561051482917356396 E-127'        # 32-bit
+        '2.5657285005610514829173563961304785900 E-127'    # 64-bit
+        sage: old_prec = gp.set_real_precision(50)
+        sage: gp.eval('intnum(x=17,42,exp(-x^2)*log(x))')
+        '2.5657285005610514829173563961304785900147709554020 E-127'
+        sage: gp.set_real_precision(old_prec)
+        50
+
+    Note that the input function above is a string in PARI syntax.
+
+    """
+    if len(arg)==3:
+        #old numerical_integral is now renamed to gsl_numerical_integral
+        return gsl_numerical_integral(*arg, **kwds);
+    elif len(arg)==4:
+        return _numerical_integral(*arg, **kwds)
+    else:
+        return "Check you input, you have eithr not enough arguements or you entered a keyword as an arguement. You have entered ", len(arg), "arguement, but only 3 or 4 are allowed"
+        
+def _numerical_integral(f,x,a,b,alg=None, **kwds):
+    if is_Expression(f):
+        if alg==None or alg=="maxima":
+            try:
+                v = f._maxima_().quad_qags(x, a, b,
+                                    epsrel=1e-8,
+                                    limit=200)
+            except TypeError, err:
+                if "ERROR" in str(err):
+                    raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
+                else:
+                    raise TypeError, err
+            if 'quad_qags' in str(v):
+                raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
+
+            return float(v[0])
+
+
+        elif alg=="sympy":
+            return numerical_approx(sympy_integrator(f,x,a,b))
+        elif alg=="mathematica":
+            return numerical_approx(mma_free_integrator(f,x,a,b))
+           
+        elif alg=="gsl":
+            return gsl_numerical_integral(f,a,b, algorithm='qag', max_points=87, params=[], eps_abs=1e-6, eps_rel=1e-6, rule=6)[0]
+        else:
+            return "input function is symbolic, and currently maxima and gsl is the option for algorithm"
+    elif alg=="mpmath" or alg==None:
+    
+        def mpintegral(f,a,b,d=15):
+            mp.dps = d; mp.pretty = True
+        return quad(f, [a, b])
+        return mpintegral(f,a,b,d)
+        
+    else:
+        return "Algorithm either not implemented or misspelled"
+
+
+def _nintegral_sym(ex, x, a, b,
               desired_relative_error='1e-8',
               maximum_num_subintervals=200):
     r"""
@@ -687 +858 @@
 
     Note that the input function above is a string in PARI syntax.
     """
+    
     try:
         v = ex._maxima_().quad_qags(x, a, b,
                                     epsrel=desired_relative_error,
@@ -703 +875 @@
         raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
 
     return float(v[0]), float(v[1]), Integer(v[2]), Integer(v[3])
-
-nintegrate = nintegral
+    
+nintegrate = numerical_integral
+nintegral = numerical_integral
 
 def minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0):
     r"""

sage/functions/transcendental.py

@@ -24 +24 @@
 import sage.rings.complex_field as complex_field
 import sage.rings.real_double as real_double
 import sage.rings.complex_number
-from sage.gsl.integration import numerical_integral
+from sage.gsl.integration import gsl_numerical_integral
 from sage.structure.parent import Parent
 from sage.structure.coerce import parent
 from sage.symbolic.expression import Expression
@@ -344 +344 @@
             sage: t.prec()
             100
 
-            sage: Ei(1.1, prec=300)
+            sage: Ei(1.1, prec=30)
             doctest:...: DeprecationWarning: The prec keyword argument is deprecated. Explicitly set the precision of the input, for example Ei(RealField(300)(1)), or use the prec argument to .n() for exact inputs, e.g., Ei(1).n(300), instead.
-            2.16737827956340306615064476647912607220394065907142504328679588538509331805598360907980986
+            2.1673783
+
         """
         if prec is not None:
             from sage.misc.misc import deprecation
@@ -450 +451 @@
     if x < 2:
         raise ValueError, "Li only defined for x at least 2."
     if eps_rel:
-        ans = numerical_integral(_one_over_log, 2, float(x),
+        ans = gsl_numerical_integral(_one_over_log, 2, float(x),
                              eps_rel=eps_rel)
     else:
-        ans = numerical_integral(_one_over_log, 2, float(x))
+        ans = gsl_numerical_integral(_one_over_log, 2, float(x))
     if err_bound:
         return real_double.RDF(ans[0]), ans[1]
     else:

sage/gsl/all.py

@@ -11 +11 @@
 from ode import ode_solver
 from ode import ode_system
 from probability_distribution import RealDistribution
-from integration import numerical_integral
-integral_numerical = numerical_integral
+from integration import gsl_numerical_integral
+#integral_numerical = numerical_integral
 from probability_distribution import SphericalDistribution
 from probability_distribution import GeneralDiscreteDistribution

sage/gsl/integration.pyx

@@ -62 +62 @@
     return (<FastDoubleFunc>params)._call_c(&t)
 
 
-def numerical_integral(func, a, b=None,
+def gsl_numerical_integral(func, a, b=None,
                        algorithm='qag',
                        max_points=87, params=[], eps_abs=1e-6,
                        eps_rel=1e-6, rule=6):
@@ -126 +126 @@
 
    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)]   # slightly random output (architecture and os dependent)
+      sage: [gsl_numerical_integral(f, 1, 2, max_points=100, params=[n]) for n in range(10)]   # slightly random output (architecture and os dependent)
       [(0.49999999999998657, 5.5511151231256336e-15),
        (0.32175055439664557, 3.5721487367706477e-15),
        (0.24030098317249229, 2.6678768435816325e-15),

sage/symbolic/expression.pyx

@@ -8158 +8158 @@
             sage: sin(x).nintegral(x,0,3)
             (1.989992496600..., 2.209335488557...e-14, 21, 0)
         """
-        from sage.calculus.calculus import nintegral
+        from sage.calculus.calculus import _nintegral_sym as nintegral
         return nintegral(self, *args, **kwds)
 
     nintegrate = nintegral

sage/symbolic/integration/integral.py

@@ -195 +195 @@
             sage: integrate(x^2.7 * e^(-2.4*x), x, 0, 3).n()
             0.154572952320790
         """
-        from sage.gsl.integration import numerical_integral
+        from sage.calculus.calculus import numerical_integral
         # The gsl routine returns a tuple, which also contains the error.
         # We only return the result.
         return numerical_integral(f, a, b)[0]