Ticket #7763: integration.patch
File integration.patch, 14.3 KB (added by , 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 1 1 from calculus import maxima as maxima_calculus 2 2 from calculus import (laplace, inverse_laplace, 3 limit, lim )3 limit, lim, numerical_integral, nintegral, nintegrate) 4 4 5 5 from functional import (diff, derivative, 6 6 expand, -
sage/calculus/calculus.py
diff -r c58f82863c75 -r 57fc8c33a0d6 sage/calculus/calculus.py
a b 15 15 16 16 - Tom Coates (2010-06-11): fixed Trac #9217 17 17 18 - Gagan Sekhon (2010-01): Added numerical_integral to fix trac 7763 19 18 20 The Sage calculus module is loosely based on the Sage Enhancement 19 21 Proposal found at: http://www.sagemath.org:9001/CalculusSEP. 20 22 … … 358 360 sage: integrate(sqrt(cos(x)^2 + sin(x)^2), x, 0, 2*pi) 359 361 2*pi 360 362 """ 363 from sage.symbolic.integration.external import maxima_integrator, sympy_integrator, mma_free_integrator 364 from sage.gsl.integration import gsl_numerical_integral 365 from mpmath import quad, mp 361 366 362 367 import re 363 368 from sage.rings.all import RR, Integer, CC, QQ, RealDoubleElement, algdep … … 366 371 from sage.misc.latex import latex, latex_variable_name 367 372 from sage.interfaces.maxima import Maxima 368 373 from sage.misc.parser import Parser 374 from sage.misc.functional import numerical_approx 369 375 370 376 from sage.symbolic.ring import var, SR, is_SymbolicVariable 371 from sage.symbolic.expression import Expression 377 from sage.symbolic.expression import Expression, is_Expression 372 378 from sage.symbolic.function import Function 373 379 from sage.symbolic.function_factory import function_factory, function 374 380 from sage.symbolic.integration.integral import indefinite_integral, \ … … 559 565 else: 560 566 raise ValueError, "unknown algorithm: %s" % algorithm 561 567 562 def nintegral(ex, x, a, b, 568 #This is so the old numerical_integral will still work 569 def 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 695 def _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 733 def _nintegral_sym(ex, x, a, b, 563 734 desired_relative_error='1e-8', 564 735 maximum_num_subintervals=200): 565 736 r""" … … 687 858 688 859 Note that the input function above is a string in PARI syntax. 689 860 """ 861 690 862 try: 691 863 v = ex._maxima_().quad_qags(x, a, b, 692 864 epsrel=desired_relative_error, … … 703 875 raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision" 704 876 705 877 return float(v[0]), float(v[1]), Integer(v[2]), Integer(v[3]) 706 707 nintegrate = nintegral 878 879 nintegrate = numerical_integral 880 nintegral = numerical_integral 708 881 709 882 def minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0): 710 883 r""" -
sage/functions/transcendental.py
diff -r c58f82863c75 -r 57fc8c33a0d6 sage/functions/transcendental.py
a b 24 24 import sage.rings.complex_field as complex_field 25 25 import sage.rings.real_double as real_double 26 26 import sage.rings.complex_number 27 from sage.gsl.integration import numerical_integral27 from sage.gsl.integration import gsl_numerical_integral 28 28 from sage.structure.parent import Parent 29 29 from sage.structure.coerce import parent 30 30 from sage.symbolic.expression import Expression … … 344 344 sage: t.prec() 345 345 100 346 346 347 sage: Ei(1.1, prec=30 0)347 sage: Ei(1.1, prec=30) 348 348 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 350 351 """ 351 352 if prec is not None: 352 353 from sage.misc.misc import deprecation … … 450 451 if x < 2: 451 452 raise ValueError, "Li only defined for x at least 2." 452 453 if eps_rel: 453 ans = numerical_integral(_one_over_log, 2, float(x),454 ans = gsl_numerical_integral(_one_over_log, 2, float(x), 454 455 eps_rel=eps_rel) 455 456 else: 456 ans = numerical_integral(_one_over_log, 2, float(x))457 ans = gsl_numerical_integral(_one_over_log, 2, float(x)) 457 458 if err_bound: 458 459 return real_double.RDF(ans[0]), ans[1] 459 460 else: -
sage/gsl/all.py
diff -r c58f82863c75 -r 57fc8c33a0d6 sage/gsl/all.py
a b 11 11 from ode import ode_solver 12 12 from ode import ode_system 13 13 from probability_distribution import RealDistribution 14 from integration import numerical_integral15 integral_numerical = numerical_integral14 from integration import gsl_numerical_integral 15 #integral_numerical = numerical_integral 16 16 from probability_distribution import SphericalDistribution 17 17 from probability_distribution import GeneralDiscreteDistribution -
sage/gsl/integration.pyx
diff -r c58f82863c75 -r 57fc8c33a0d6 sage/gsl/integration.pyx
a b 62 62 return (<FastDoubleFunc>params)._call_c(&t) 63 63 64 64 65 def numerical_integral(func, a, b=None,65 def gsl_numerical_integral(func, a, b=None, 66 66 algorithm='qag', 67 67 max_points=87, params=[], eps_abs=1e-6, 68 68 eps_rel=1e-6, rule=6): … … 126 126 127 127 For a Python function with parameters: 128 128 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) 130 130 [(0.49999999999998657, 5.5511151231256336e-15), 131 131 (0.32175055439664557, 3.5721487367706477e-15), 132 132 (0.24030098317249229, 2.6678768435816325e-15), -
sage/symbolic/expression.pyx
diff -r c58f82863c75 -r 57fc8c33a0d6 sage/symbolic/expression.pyx
a b 8158 8158 sage: sin(x).nintegral(x,0,3) 8159 8159 (1.989992496600..., 2.209335488557...e-14, 21, 0) 8160 8160 """ 8161 from sage.calculus.calculus import nintegral8161 from sage.calculus.calculus import _nintegral_sym as nintegral 8162 8162 return nintegral(self, *args, **kwds) 8163 8163 8164 8164 nintegrate = nintegral -
sage/symbolic/integration/integral.py
diff -r c58f82863c75 -r 57fc8c33a0d6 sage/symbolic/integration/integral.py
a b 195 195 sage: integrate(x^2.7 * e^(-2.4*x), x, 0, 3).n() 196 196 0.154572952320790 197 197 """ 198 from sage. gsl.integrationimport numerical_integral198 from sage.calculus.calculus import numerical_integral 199 199 # The gsl routine returns a tuple, which also contains the error. 200 200 # We only return the result. 201 201 return numerical_integral(f, a, b)[0]