Ticket #7763: trac_7763.8.patch

File trac_7763.8.patch, 19.2 KB (added by gagansekhon, 9 years ago)
  • sage/calculus/all.py

    # HG changeset patch
    # User Gagan Sekhon <gagan.d.sekhon@gmail.com>
    # Date 1294883391 28800
    # Node ID ba8c6b036de215b12ba8c4d5a3c8959246caf25c
    # Parent  4a9f7b41ae22e247c2a6675c7854e419f7340d09
    trac 7763, unify numerical integration as a top level function
    
    diff -r 4a9f7b41ae22 -r ba8c6b036de2 sage/calculus/all.py
    a b  
    11from calculus import maxima as maxima_calculus
    22from calculus import (laplace, inverse_laplace,
    3                       limit, lim)
     3                      limit, lim, numerical_integral)
     4nintegral=numerical_integral
     5nintegrate=numerical_integral                     
     6                     
    47
    58from functional import (diff, derivative,
    69                        expand,
  • sage/calculus/calculus.py

    diff -r 4a9f7b41ae22 -r ba8c6b036de2 sage/calculus/calculus.py
    a b  
    1515
    1616- Tom Coates (2010-06-11): fixed Trac #9217
    1717
     18- Gagan Sekhon (2010-01): Added numerical_integral to fix trac 7763
     19
    1820The Sage calculus module is loosely based on the Sage Enhancement
    19 Proposal found at: http://www.sagemath.org:9001/CalculusSEP.
     21Proposal found at:http://www.sagemath.org:9001/CalculusSEP.
    2022
    2123EXAMPLES:
    2224
     
    358360    sage: integrate(sqrt(cos(x)^2 + sin(x)^2), x, 0, 2*pi)
    359361    2*pi
    360362"""
     363from  sage.symbolic.integration.external import maxima_integrator, sympy_integrator, mma_free_integrator
     364from sage.gsl.integration import gsl_numerical_integral
    361365
    362366import re
    363367from sage.rings.all import RR, Integer, CC, QQ, RealDoubleElement, algdep
     
    366370from sage.misc.latex import latex, latex_variable_name
    367371from sage.interfaces.maxima import Maxima
    368372from sage.misc.parser import Parser
     373from sage.misc.functional import numerical_approx
    369374
    370375from sage.symbolic.ring import var, SR, is_SymbolicVariable
    371 from sage.symbolic.expression import Expression
     376from sage.symbolic.expression import Expression, is_Expression
    372377from sage.symbolic.function import Function
    373378from sage.symbolic.function_factory import function_factory, function
    374379from sage.symbolic.integration.integral import indefinite_integral, \
     
    559564    else:
    560565        raise ValueError, "unknown algorithm: %s" % algorithm
    561566
    562 def nintegral(ex, x, a, b,
     567#This is so the function numerical_integral takes 2 different types of arguments, (f,a,b) which was the old syntax
     568#and (f,x,a,b), which is the new syntax to make is more like integral function.
     569#If len(arg) is 3, it uses gsl_numerical_integral, which was used before, and if len(arg)=4, then it calls _numerical_integral
     570#which has more options. 
     571def numerical_integral(*arg, **kwds):
     572    r"""
     573   
     574    Input:
     575   
     576    - 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$.
     577         
     578    - 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.
     579   
     580            - If `f` is symbolic expression then possible algorithms are "maxima" (default), "sympy","gsl".
     581            - If `f` is not symbolic expression then possible algorithms is "mpmath".
     582           
     583   
     584      - Depending on the algorithm selected, appropriate parameters
     585                    -  "maxima" - This algorithm has additional parameters
     586                          -  ``epsrel`` - desired_relative_error  (default: '1e-8') the desired relative error
     587
     588                          -  ``limit`` - maximum_num_subintervals (default: 200) maxima number of subintervals
     589                               
     590                    -  "sympy" - This algorithm has no additional parameters
     591                    -  "gsl" - This algorithm has additional parameters
     592                          -- algorithm: valid choices are
     593                            'qag' -- for an adaptive integration
     594                            'qng' -- for a non-adaptive gauss kronrod (samples at a maximum of 87pts)
     595                            -- max_points: sets the maximum number of sample points
     596                            -- params: used to pass parameters to your function
     597                            -- eps_abs, eps_rel: absolute and relative error tolerances
     598                            -- rule: This controls the Gauss-Kronrod rule used in the adaptive integration
     599                                    rule=1: 15 pt rule
     600                                    rule=2: 21 pt rule
     601                                    rule=3: 31 pt rule
     602                                    rule=4: 41 pt rule
     603                                    rule=5: 51 pt rule
     604                                    rule=6: 61 pt rule
     605                            Higher key values are more accurate for smooth functions but lower
     606                            key values deal better with discontinuities.
     607                            Note: For more detailed documentation for this option see gsl_numerical_integral
     608                    -  "mpmath" - This algorithm has additional parameters
     609                        -d: Decimal places
     610                       
     611     ALIAS: nintegrate and nintegral are the same as nintegral .
     612       
     613     Please also note that numerical_integral(f,x,a,b) is different than f.nintegral(x,a,b) (for now atleast).
     614     
     615           
     616           
     617   
     618   
     619    EXAMPLES::
     620
     621        sage: f(x) = x
     622        sage: nintegral(f,x,0,1,alg="maxima", epsrel='1e-14')
     623        0.5
     624        sage: numerical_integral(lambda x: x^2,x,0,1,alg="mpmath")
     625        mpf('0.33333333333333331')
     626        sage: numerical_integral(x^2,x,0,1,alg="gsl")
     627        0.33333333333333331
     628        sage: numerical_integral(x^2,x,0,1,alg="gsl",algorithm='qag',max_points=40)
     629        0.33333333333333331
     630        sage: numerical_integral(x^2,x,0,1)
     631        0.33333333333333343
     632        sage: numerical_integral(x^2,x,0,1,alg="sympy")
     633        0.333333333333333
     634        sage: f = x^2
     635        sage: numerical_integral(f, x,0, 1,alg="gsl", max_points=200, eps_abs=1e-7, eps_rel=1e-7, rule=4)
     636        0.33333333333333331
     637
     638        For a Python function with parameters:
     639        sage: f(x,a) = 1/(a+x^2)
     640        sage: [numerical_integral(f,x, 1, 2,alg="gsl", max_points=100, params=[n]) for n in range(10)]   # slightly random output (architecture and os dependent)
     641        [0.49999999999998657,
     642        0.32175055439664557, 
     643        0.24030098317249229,
     644        0.19253082576711697,
     645        0.16087527719832367,
     646        0.13827545676349412,
     647        0.12129975935702741,
     648        0.10806674191683065,
     649        0.09745444625548845,
     650        0.088750683050217577]       
     651
     652        sage:  numerical_integral(x^2,x,0, 1,alg="maxima",epsrel=1e-8,limit=300)
     653        0.33333333333333343
     654       
     655       
     656    In this example we try to use gsl algorithm even though the function is non-symbolic.::
     657       
     658        sage: numerical_integral(lambda x: x^2,x,0,1,alg="gsl")
     659        'Algorithm either not implemented or misspelled'
     660
     661
     662    """
     663    if len(arg)==3:
     664        #old numerical_integral is now renamed to gsl_numerical_integral
     665        return gsl_numerical_integral(*arg, **kwds);
     666    elif len(arg)==4:
     667        return _numerical_integral(*arg, **kwds)
     668    else:
     669        return "Check you input, you have either not enough arguements or too many arguements. The number of arguements you entered is ", len(arg), "arguement, but only 3 or 4 are allowed"
     670       
     671def _numerical_integral(f,x,a,b,alg=None, **kwds):
     672    r"""
     673    This is the new numerical integration function.
     674
     675    Input:
     676 
     677    - arg - General arguments for an integral function such as (f,x,a,b).
     678 
     679    - alg - Algorithm to be used. For example: alg="maxima". Default algorithm "maxima" for symbolic expression and "mpmath" for non-symbolic expression.
     680            - If `f` is symbolic expression then possible algorithms are "maxima" (default), "sympy", "mathematica","gsl".
     681            - If `f` is not symbolic expression then possible algorithms is "mpmath".
     682            Note: Reason for calling it alg instead of algorithm is that "gsl" algorithm already has a parameter called algorithm.
     683    -kwds -
     684
     685   - Depending on the algorithm selected, appropriate parameters
     686                    -  "maxima" - This algorithm has additional parameters
     687                          -  ``epsrel`` - desired_relative_error  (default: '1e-8') the desired relative error
     688
     689                          -  ``limit`` - maximum_num_subintervals (default: 200) maxima number of subintervals
     690                               
     691                    -  "sympy" - This algorithm has no additional parameters
     692                    -  "gsl" - This algorithm has additional parameters
     693                          -- algorithm: valid choices are
     694                            'qag' -- for an adaptive integration
     695                            'qng' -- for a non-adaptive gauss kronrod (samples at a maximum of 87pts)
     696                            -- max_points: sets the maximum number of sample points
     697                            -- params: used to pass parameters to your function
     698                            -- eps_abs, eps_rel: absolute and relative error tolerances
     699                            -- rule: This controls the Gauss-Kronrod rule used in the adaptive integration
     700                                    rule=1: 15 pt rule
     701                                    rule=2: 21 pt rule
     702                                    rule=3: 31 pt rule
     703                                    rule=4: 41 pt rule
     704                                    rule=5: 51 pt rule
     705                                    rule=6: 61 pt rule
     706                            Higher key values are more accurate for smooth functions but lower
     707                            key values deal better with discontinuities.
     708                    -  "mpmath" - This algorithm has additional parameters
     709                        -d: Decimal places
     710                       
     711     ALIAS: nintegrate and nintegral are the same as nintegral .
     712       
     713     Please also note that numerical_integral(f,x,a,b) is different than f.nintegral(x,a,b) (for now atleast).
     714
     715    """
     716
     717    if is_Expression(f):
     718        if alg==None or alg=="maxima":
     719            try:
     720                v = f._maxima_().quad_qags(x, a, b,
     721                                    epsrel=1e-8,
     722                                    limit=200)
     723            except TypeError, err:
     724                if "ERROR" in str(err):
     725                    raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
     726                else:
     727                    raise TypeError, err
     728            if 'quad_qags' in str(v):
     729                raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
     730
     731            return float(v[0])
     732
     733
     734        elif alg=="sympy":
     735            return numerical_approx(sympy_integrator(f,x,a,b))
     736        elif alg=="gsl":
     737            return gsl_numerical_integral(f,a,b, **kwds)[0]
     738        else:
     739            return "input function is symbolic, and currently maxima and gsl is the option for algorithm"
     740    elif alg=="mpmath" or alg==None:
     741        from  mpmath import quad,mp, mpf
     742        def mpintegral(f,a,b,d):
     743            if d==None:
     744                d=15
     745            mp.dps=d; mp.pretty=True
     746        return quad(f, [a, b])
     747        return mpintegral(f,a,b,d)
     748       
     749    else:
     750        return "Algorithm either not implemented or misspelled"
     751
     752
     753def _nintegral_sym(ex, x, a, b,
    563754              desired_relative_error='1e-8',
    564755              maximum_num_subintervals=200):
    565756    r"""
     
    621812    Also, there are limits to the precision to which Maxima can compute
    622813    the integral to due to limitations in quadpack.
    623814
    624     ::
     815     EXAMPLES::
     816
     817        sage: f(x) = exp(-sqrt(x))
     818        sage: f.nintegral(x, 0, 1)
     819        (0.52848223531423055, 4.1633141378838452e-11, 231, 0)
     820
     821
     822
     823    Note that in exotic cases where floating point evaluation of the
     824    expression leads to the wrong value, then the output can be
     825    completely wrong::
     826
     827        sage: f = exp(pi*sqrt(163)) - 262537412640768744
     828
     829    Despite appearance, `f` is really very close to 0, but one
     830    gets a nonzero value since the definition of
     831    ``float(f)`` is that it makes all constants inside the
     832    expression floats, then evaluates each function and each arithmetic
     833    operation using float arithmetic::
     834
     835        sage: float(f)
     836        -480.0
     837
     838    Computing to higher precision we see the truth::
     839
     840        sage: f.n(200)
     841        -7.4992740280181431112064614366622348652078895136533593355718e-13
     842        sage: f.n(300)
     843        -7.49927402801814311120646143662663009137292462589621789352095066181709095575681963967103004e-13
     844
     845    Now numerically integrating, we see why the answer is wrong::
     846
     847        sage: f.nintegrate(x,0,1)
     848        (-480.00000000000011, 5.3290705182007538e-12, 21, 0)
     849         
     850    It is just because every floating point evaluation of return -480.0
     851    in floating point.
     852
     853    Important note: using PARI/GP one can compute numerical integrals
     854    to high precision::
     855
     856        sage: gp.eval('intnum(x=17,42,exp(-x^2)*log(x))')
     857        '2.565728500561051482917356396 E-127'        # 32-bit
     858        '2.5657285005610514829173563961304785900 E-127'    # 64-bit
     859        sage: old_prec = gp.set_real_precision(50)
     860        sage: gp.eval('intnum(x=17,42,exp(-x^2)*log(x))')
     861        '2.5657285005610514829173563961304785900147709554020 E-127'
     862        sage: gp.set_real_precision(old_prec)
     863        50
     864
     865    Note that the input function above is a string in PARI syntax.::
    625866
    626867        sage: f = x
    627868        sage: f = f.nintegral(x,0,1,1e-14)
     
    635876        sage: f.nintegral(x, 0, 1)
    636877        (0.52848223531423055, 4.163...e-11, 231, 0)
    637878
    638     We can also use the ``numerical_integral`` function,
    639     which calls the GSL C library.
    640 
    641     ::
    642 
    643         sage: numerical_integral(f, 0, 1)
    644         (0.52848223225314706, 6.83928460...e-07)
    645 
    646879    Note that in exotic cases where floating point evaluation of the
    647880    expression leads to the wrong value, then the output can be
    648881    completely wrong::
     
    687920
    688921    Note that the input function above is a string in PARI syntax.
    689922    """
     923   
    690924    try:
    691925        v = ex._maxima_().quad_qags(x, a, b,
    692926                                    epsrel=desired_relative_error,
     
    703937        raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
    704938
    705939    return float(v[0]), float(v[1]), Integer(v[2]), Integer(v[3])
    706 
    707 nintegrate = nintegral
     940   
    708941
    709942def minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0):
    710943    r"""
  • sage/functions/transcendental.py

    diff -r 4a9f7b41ae22 -r ba8c6b036de2 sage/functions/transcendental.py
    a b  
    2424import sage.rings.complex_field as complex_field
    2525import sage.rings.real_double as real_double
    2626import sage.rings.complex_number
    27 from sage.gsl.integration import numerical_integral
     27from sage.gsl.integration import gsl_numerical_integral
    2828from sage.structure.parent import Parent
    2929from sage.structure.coerce import parent
    3030from sage.symbolic.expression import Expression
     
    344344            sage: t.prec()
    345345            100
    346346
    347             sage: Ei(1.1, prec=300)
     347            sage: Ei(1.1, prec=30)
    348348            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.
    349             2.16737827956340306615064476647912607220394065907142504328679588538509331805598360907980986
     349            2.1673783
     350
    350351        """
    351352        if prec is not None:
    352353            from sage.misc.misc import deprecation
     
    450451    if x < 2:
    451452        raise ValueError, "Li only defined for x at least 2."
    452453    if eps_rel:
    453         ans = numerical_integral(_one_over_log, 2, float(x),
     454        ans = gsl_numerical_integral(_one_over_log, 2, float(x),
    454455                             eps_rel=eps_rel)
    455456    else:
    456         ans = numerical_integral(_one_over_log, 2, float(x))
     457        ans = gsl_numerical_integral(_one_over_log, 2, float(x))
    457458    if err_bound:
    458459        return real_double.RDF(ans[0]), ans[1]
    459460    else:
  • sage/gsl/all.py

    diff -r 4a9f7b41ae22 -r ba8c6b036de2 sage/gsl/all.py
    a b  
    1111from ode import ode_solver
    1212from ode import ode_system
    1313from probability_distribution import RealDistribution
    14 from integration import numerical_integral
    15 integral_numerical = numerical_integral
     14#from integration import gsl_numerical_integral
     15#integral_numerical = numerical_integral
    1616from probability_distribution import SphericalDistribution
    1717from probability_distribution import GeneralDiscreteDistribution
  • sage/gsl/integration.pyx

    diff -r 4a9f7b41ae22 -r ba8c6b036de2 sage/gsl/integration.pyx
    a b  
    6262    return (<FastDoubleFunc>params)._call_c(&t)
    6363
    6464
    65 def numerical_integral(func, a, b=None,
     65def gsl_numerical_integral(func, a, b=None,
    6666                       algorithm='qag',
    6767                       max_points=87, params=[], eps_abs=1e-6,
    6868                       eps_rel=1e-6, rule=6):
     
    126126
    127127   For a Python function with parameters:
    128128      sage: f(x,a) = 1/(a+x^2)
    129       sage: [numerical_integral(f, 1, 2, max_points=100, params=[n]) for n in range(10)]   # slightly random output (architecture and os dependent)
     129      sage: [gsl_numerical_integral(f, 1, 2, max_points=100, params=[n]) for n in range(10)]   # slightly random output (architecture and os dependent)
    130130      [(0.49999999999998657, 5.5511151231256336e-15),
    131131       (0.32175055439664557, 3.5721487367706477e-15),
    132132       (0.24030098317249229, 2.6678768435816325e-15),
  • sage/symbolic/expression.pyx

    diff -r 4a9f7b41ae22 -r ba8c6b036de2 sage/symbolic/expression.pyx
    a b  
    81598159            sage: sin(x).nintegral(x,0,3)
    81608160            (1.989992496600..., 2.209335488557...e-14, 21, 0)
    81618161        """
    8162         from sage.calculus.calculus import nintegral
     8162        from sage.calculus.calculus import _nintegral_sym as nintegral
    81638163        return nintegral(self, *args, **kwds)
    81648164
    81658165    nintegrate = nintegral
  • sage/symbolic/integration/integral.py

    diff -r 4a9f7b41ae22 -r ba8c6b036de2 sage/symbolic/integration/integral.py
    a b  
    195195            sage: integrate(x^2.7 * e^(-2.4*x), x, 0, 3).n()
    196196            0.154572952320790
    197197        """
    198         from sage.gsl.integration import numerical_integral
     198        from sage.calculus.calculus import numerical_integral
    199199        # The gsl routine returns a tuple, which also contains the error.
    200200        # We only return the result.
    201201        return numerical_integral(f, a, b)[0]