Ticket #7763: integration.patch

File integration.patch, 14.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 57fc8c33a0d676189977babd349a7bd7ed8e50da
    # Parent  c58f82863c75b2e53662e08cd550bb398653174c
    [mq]: numerical_integral
    
    diff -r c58f82863c75 -r 57fc8c33a0d6 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, nintegral, nintegrate)
    44
    55from functional import (diff, derivative,
    66                        expand,
  • sage/calculus/calculus.py

    diff -r c58f82863c75 -r 57fc8c33a0d6 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
     365from mpmath import quad, mp
    361366
    362367import re
    363368from sage.rings.all import RR, Integer, CC, QQ, RealDoubleElement, algdep
     
    366371from sage.misc.latex import latex, latex_variable_name
    367372from sage.interfaces.maxima import Maxima
    368373from sage.misc.parser import Parser
     374from sage.misc.functional import numerical_approx
    369375
    370376from sage.symbolic.ring import var, SR, is_SymbolicVariable
    371 from sage.symbolic.expression import Expression
     377from sage.symbolic.expression import Expression, is_Expression
    372378from sage.symbolic.function import Function
    373379from sage.symbolic.function_factory import function_factory, function
    374380from sage.symbolic.integration.integral import indefinite_integral, \
     
    559565    else:
    560566        raise ValueError, "unknown algorithm: %s" % algorithm
    561567
    562 def nintegral(ex, x, a, b,
     568#This is so the old numerical_integral will still work
     569def numerical_integral(*arg, **kwds):
     570    r"""
     571   
     572    Input:
     573   
     574    - arg - General arguments for an integral function such as (f,x,a,b). However, (f,a,b) is also valid,
     575    it just calls the old numerical_integral which is now called old gsl_numerical_integral. `f` can be
     576    symbolic expression such as `sin(x)` or not symbolic `sin`.
     577   
     578   Note: There are other methods, but I am still trying to figure out how to include them
     579        Possible algorithms are gp.intnum, Gi/Pynac
     580   
     581    - kwds - First specify algorithm specifies which algorithm is to be used. For example: alg="maxima".
     582             Default algorithm "maxima" for symbolic expression and "mpmath" for non-symbolic expression.
     583   
     584            - If `f` is symbolic expression then possible algorithms are "maxima" (default), "sympy",
     585             "mathematica","gsl".
     586            - If `f` is not symbolic expression then possible algorithms is "mpmath".
     587           
     588   
     589      - Depending on the algorithm selected, appropriate parameters
     590                    -  "maxima" - This algorithm has additional parameters
     591                          -  ``epsrel`` - desired_relative_error  (default: '1e-8') the
     592                                desired relative error
     593
     594                          -  ``limit`` - maximum_num_subintervals (default: 200)
     595                                maxima number of subintervals
     596                               
     597                    -  "sympy" - This algorithm has no additional parameters
     598                    -  "mathematica" - This algorithm has no additional parameters
     599                    -  "gsl" - This algorithm has additional parameters
     600                          -- algorithm: valid choices are
     601                            'qag' -- for an adaptive integration
     602                            'qng' -- for a non-adaptive gauss kronrod (samples at a maximum of 87pts)
     603                            -- max_points: sets the maximum number of sample points
     604                            -- params: used to pass parameters to your function
     605                            -- eps_abs, eps_rel: absolute and relative error tolerances
     606                            -- rule: This controls the Gauss-Kronrod rule used in the adaptive integration
     607                                    rule=1: 15 pt rule
     608                                    rule=2: 21 pt rule
     609                                    rule=3: 31 pt rule
     610                                    rule=4: 41 pt rule
     611                                    rule=5: 51 pt rule
     612                                    rule=6: 61 pt rule
     613                            Higher key values are more accurate for smooth functions but lower
     614                            key values deal better with discontinuities.
     615                    -  "mpmath" - This algorithm has additional parameters
     616                        -d: Decimal places
     617                       
     618     ALIAS: nintegrate and nintegral are the same as nintegral   
     619     
     620           
     621    ::
     622
     623        sage: f(x) = x
     624        sage: nintegral(f,x,0,1, epsrel='1e-14')
     625        0.5
     626       
     627    EXAMPLES::
     628
     629        sage: f(x) = exp(-sqrt(x))
     630        sage: f.nintegral(x, 0, 1)
     631        (0.52848223531423055, 4.1633141378838452e-11, 231, 0)
     632
     633
     634    We can also use the ``numerical_integral`` function,
     635    which calls the GSL C library.
     636
     637    ::
     638
     639        sage: numerical_integral(f, 0, 1)
     640        (0.52848223225314706, 6.83928460...e-07)
     641
     642    Note that in exotic cases where floating point evaluation of the
     643    expression leads to the wrong value, then the output can be
     644    completely wrong::
     645
     646        sage: f = exp(pi*sqrt(163)) - 262537412640768744
     647
     648    Despite appearance, `f` is really very close to 0, but one
     649    gets a nonzero value since the definition of
     650    ``float(f)`` is that it makes all constants inside the
     651    expression floats, then evaluates each function and each arithmetic
     652    operation using float arithmetic::
     653
     654        sage: float(f)
     655        -480.0
     656
     657    Computing to higher precision we see the truth::
     658
     659        sage: f.n(200)
     660        -7.4992740280181431112064614366622348652078895136533593355718e-13
     661        sage: f.n(300)
     662        -7.49927402801814311120646143662663009137292462589621789352095066181709095575681963967103004e-13
     663
     664    Now numerically integrating, we see why the answer is wrong::
     665
     666        sage: f.nintegrate(x,0,1)
     667        (-480.00000000000011, 5.3290705182007538e-12, 21, 0)
     668         
     669    It is just because every floating point evaluation of return -480.0
     670    in floating point.
     671
     672    Important note: using PARI/GP one can compute numerical integrals
     673    to high precision::
     674
     675        sage: gp.eval('intnum(x=17,42,exp(-x^2)*log(x))')
     676        '2.565728500561051482917356396 E-127'        # 32-bit
     677        '2.5657285005610514829173563961304785900 E-127'    # 64-bit
     678        sage: old_prec = gp.set_real_precision(50)
     679        sage: gp.eval('intnum(x=17,42,exp(-x^2)*log(x))')
     680        '2.5657285005610514829173563961304785900147709554020 E-127'
     681        sage: gp.set_real_precision(old_prec)
     682        50
     683
     684    Note that the input function above is a string in PARI syntax.
     685
     686    """
     687    if len(arg)==3:
     688        #old numerical_integral is now renamed to gsl_numerical_integral
     689        return gsl_numerical_integral(*arg, **kwds);
     690    elif len(arg)==4:
     691        return _numerical_integral(*arg, **kwds)
     692    else:
     693        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"
     694       
     695def _numerical_integral(f,x,a,b,alg=None, **kwds):
     696    if is_Expression(f):
     697        if alg==None or alg=="maxima":
     698            try:
     699                v = f._maxima_().quad_qags(x, a, b,
     700                                    epsrel=1e-8,
     701                                    limit=200)
     702            except TypeError, err:
     703                if "ERROR" in str(err):
     704                    raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
     705                else:
     706                    raise TypeError, err
     707            if 'quad_qags' in str(v):
     708                raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
     709
     710            return float(v[0])
     711
     712
     713        elif alg=="sympy":
     714            return numerical_approx(sympy_integrator(f,x,a,b))
     715        elif alg=="mathematica":
     716            return numerical_approx(mma_free_integrator(f,x,a,b))
     717           
     718        elif alg=="gsl":
     719            return gsl_numerical_integral(f,a,b, algorithm='qag', max_points=87, params=[], eps_abs=1e-6, eps_rel=1e-6, rule=6)[0]
     720        else:
     721            return "input function is symbolic, and currently maxima and gsl is the option for algorithm"
     722    elif alg=="mpmath" or alg==None:
     723   
     724        def mpintegral(f,a,b,d=15):
     725            mp.dps = d; mp.pretty = True
     726        return quad(f, [a, b])
     727        return mpintegral(f,a,b,d)
     728       
     729    else:
     730        return "Algorithm either not implemented or misspelled"
     731
     732
     733def _nintegral_sym(ex, x, a, b,
    563734              desired_relative_error='1e-8',
    564735              maximum_num_subintervals=200):
    565736    r"""
     
    687858
    688859    Note that the input function above is a string in PARI syntax.
    689860    """
     861   
    690862    try:
    691863        v = ex._maxima_().quad_qags(x, a, b,
    692864                                    epsrel=desired_relative_error,
     
    703875        raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
    704876
    705877    return float(v[0]), float(v[1]), Integer(v[2]), Integer(v[3])
    706 
    707 nintegrate = nintegral
     878   
     879nintegrate = numerical_integral
     880nintegral = numerical_integral
    708881
    709882def minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0):
    710883    r"""
  • sage/functions/transcendental.py

    diff -r c58f82863c75 -r 57fc8c33a0d6 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 c58f82863c75 -r 57fc8c33a0d6 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 c58f82863c75 -r 57fc8c33a0d6 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 c58f82863c75 -r 57fc8c33a0d6 sage/symbolic/expression.pyx
    a b  
    81588158            sage: sin(x).nintegral(x,0,3)
    81598159            (1.989992496600..., 2.209335488557...e-14, 21, 0)
    81608160        """
    8161         from sage.calculus.calculus import nintegral
     8161        from sage.calculus.calculus import _nintegral_sym as nintegral
    81628162        return nintegral(self, *args, **kwds)
    81638163
    81648164    nintegrate = nintegral
  • sage/symbolic/integration/integral.py

    diff -r c58f82863c75 -r 57fc8c33a0d6 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]