Ticket #7763: trac_7763.7.patch

File trac_7763.7.patch, 15.3 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 57113673f32b479ae5c815a491530d6a0abbd1b1
    # Parent  4a9f7b41ae22e247c2a6675c7854e419f7340d09
    [mq]: numerical_integral
    
    diff -r 4a9f7b41ae22 -r 57113673f32b 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 57113673f32b 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
    1921Proposal found at: http://www.sagemath.org:9001/CalculusSEP.
    2022
     
    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 old numerical_integral will still work
     568def numerical_integral(*arg, **kwds):
     569    r"""
     570   
     571    Input:
     572   
     573    - 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$.
     574         
     575    - 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.
     576   
     577            - If `f` is symbolic expression then possible algorithms are "maxima" (default), "sympy", "mathematica","gsl".
     578            - If `f` is not symbolic expression then possible algorithms is "mpmath".
     579           
     580   
     581      - Depending on the algorithm selected, appropriate parameters
     582                    -  "maxima" - This algorithm has additional parameters
     583                          -  ``epsrel`` - desired_relative_error  (default: '1e-8') the desired relative error
     584
     585                          -  ``limit`` - maximum_num_subintervals (default: 200) maxima number of subintervals
     586                               
     587                    -  "sympy" - This algorithm has no additional parameters
     588                    -  "mathematica" - This algorithm has no additional parameters
     589                    -  "gsl" - This algorithm has additional parameters
     590                          -- algorithm: valid choices are
     591                            'qag' -- for an adaptive integration
     592                            'qng' -- for a non-adaptive gauss kronrod (samples at a maximum of 87pts)
     593                            -- max_points: sets the maximum number of sample points
     594                            -- params: used to pass parameters to your function
     595                            -- eps_abs, eps_rel: absolute and relative error tolerances
     596                            -- rule: This controls the Gauss-Kronrod rule used in the adaptive integration
     597                                    rule=1: 15 pt rule
     598                                    rule=2: 21 pt rule
     599                                    rule=3: 31 pt rule
     600                                    rule=4: 41 pt rule
     601                                    rule=5: 51 pt rule
     602                                    rule=6: 61 pt rule
     603                            Higher key values are more accurate for smooth functions but lower
     604                            key values deal better with discontinuities.
     605                    -  "mpmath" - This algorithm has additional parameters
     606                        -d: Decimal places
     607                       
     608     ALIAS: nintegrate and nintegral are the same as nintegral .
     609       
     610     Please also note that numerical_integral(f,x,a,b) is different than f.nintegral(x,a,b) (for now atleast).
     611     
     612           
     613           
     614   
     615   
     616    Examples::
     617
     618        sage: f(x) = x
     619        sage: nintegral(f,x,0,1, epsrel='1e-14')
     620        0.5
     621        sage: numerical_integral(lambda x: x^2,x,0,1,alg="mpmath")
     622        mpf('0.33333333333333331')
     623        sage: numerical_integral(x^2,x,0,1,alg="gsl")
     624        0.33333333333333331
     625        sage: numerical_integral(x^2,x,0,1)
     626        0.33333333333333343
     627
     628       
     629       
     630    In this example we try to use gsl algorithm even though the function is non-symbolic.::
     631       
     632        sage: numerical_integral(lambda x: x^2,x,0,1,alg="gsl")
     633        'Algorithm either not implemented or misspelled'
     634
     635
     636    """
     637    if len(arg)==3:
     638        #old numerical_integral is now renamed to gsl_numerical_integral
     639        return gsl_numerical_integral(*arg, **kwds);
     640    elif len(arg)==4:
     641        return _numerical_integral(*arg, **kwds)
     642    else:
     643        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"
     644       
     645def _numerical_integral(f,x,a,b,alg=None, **kwds):
     646    if is_Expression(f):
     647        if alg==None or alg=="maxima":
     648            try:
     649                v = f._maxima_().quad_qags(x, a, b,
     650                                    epsrel=1e-8,
     651                                    limit=200)
     652            except TypeError, err:
     653                if "ERROR" in str(err):
     654                    raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
     655                else:
     656                    raise TypeError, err
     657            if 'quad_qags' in str(v):
     658                raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
     659
     660            return float(v[0])
     661
     662
     663        elif alg=="sympy":
     664            return numerical_approx(sympy_integrator(f,x,a,b))
     665        elif alg=="mathematica":
     666            return numerical_approx(mma_free_integrator(f,x,a,b))
     667           
     668        elif alg=="gsl":
     669            return gsl_numerical_integral(f,a,b, algorithm='qag', max_points=87, params=[], eps_abs=1e-6, eps_rel=1e-6, rule=6)[0]
     670        else:
     671            return "input function is symbolic, and currently maxima and gsl is the option for algorithm"
     672    elif alg=="mpmath" or alg==None:
     673        from  mpmath import quad,mp, mpf
     674        def mpintegral(f,a,b,d):
     675            if d==None:
     676                d=15
     677            mp.dps=d; mp.pretty=True
     678        return quad(f, [a, b])
     679        return mpintegral(f,a,b,d)
     680       
     681    else:
     682        return "Algorithm either not implemented or misspelled"
     683
     684
     685def _nintegral_sym(ex, x, a, b,
    563686              desired_relative_error='1e-8',
    564687              maximum_num_subintervals=200):
    565688    r"""
     
    621744    Also, there are limits to the precision to which Maxima can compute
    622745    the integral to due to limitations in quadpack.
    623746
    624     ::
     747     EXAMPLES::
     748
     749        sage: f(x) = exp(-sqrt(x))
     750        sage: f.nintegral(x, 0, 1)
     751        (0.52848223531423055, 4.1633141378838452e-11, 231, 0)
     752
     753
     754
     755    Note that in exotic cases where floating point evaluation of the
     756    expression leads to the wrong value, then the output can be
     757    completely wrong::
     758
     759        sage: f = exp(pi*sqrt(163)) - 262537412640768744
     760
     761    Despite appearance, `f` is really very close to 0, but one
     762    gets a nonzero value since the definition of
     763    ``float(f)`` is that it makes all constants inside the
     764    expression floats, then evaluates each function and each arithmetic
     765    operation using float arithmetic::
     766
     767        sage: float(f)
     768        -480.0
     769
     770    Computing to higher precision we see the truth::
     771
     772        sage: f.n(200)
     773        -7.4992740280181431112064614366622348652078895136533593355718e-13
     774        sage: f.n(300)
     775        -7.49927402801814311120646143662663009137292462589621789352095066181709095575681963967103004e-13
     776
     777    Now numerically integrating, we see why the answer is wrong::
     778
     779        sage: f.nintegrate(x,0,1)
     780        (-480.00000000000011, 5.3290705182007538e-12, 21, 0)
     781         
     782    It is just because every floating point evaluation of return -480.0
     783    in floating point.
     784
     785    Important note: using PARI/GP one can compute numerical integrals
     786    to high precision::
     787
     788        sage: gp.eval('intnum(x=17,42,exp(-x^2)*log(x))')
     789        '2.565728500561051482917356396 E-127'        # 32-bit
     790        '2.5657285005610514829173563961304785900 E-127'    # 64-bit
     791        sage: old_prec = gp.set_real_precision(50)
     792        sage: gp.eval('intnum(x=17,42,exp(-x^2)*log(x))')
     793        '2.5657285005610514829173563961304785900147709554020 E-127'
     794        sage: gp.set_real_precision(old_prec)
     795        50
     796
     797    Note that the input function above is a string in PARI syntax.::
    625798
    626799        sage: f = x
    627800        sage: f = f.nintegral(x,0,1,1e-14)
     
    635808        sage: f.nintegral(x, 0, 1)
    636809        (0.52848223531423055, 4.163...e-11, 231, 0)
    637810
    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 
    646811    Note that in exotic cases where floating point evaluation of the
    647812    expression leads to the wrong value, then the output can be
    648813    completely wrong::
     
    687852
    688853    Note that the input function above is a string in PARI syntax.
    689854    """
     855   
    690856    try:
    691857        v = ex._maxima_().quad_qags(x, a, b,
    692858                                    epsrel=desired_relative_error,
     
    703869        raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
    704870
    705871    return float(v[0]), float(v[1]), Integer(v[2]), Integer(v[3])
    706 
    707 nintegrate = nintegral
     872   
    708873
    709874def minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0):
    710875    r"""
  • sage/functions/transcendental.py

    diff -r 4a9f7b41ae22 -r 57113673f32b 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 57113673f32b 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
     14from 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 57113673f32b 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 57113673f32b 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 57113673f32b 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]