if the method='scipy' is specified, the scipy.integrate.quad function is called on the lambdified argument of the integral.
Tests were added to sympy/test_external/test_scipy.py Signed-off-by: Ondrej Certik <ond...@certik.cz> --- sympy/core/evalf.py | 43 +++++++++++++++++++++------- sympy/core/numbers.py | 2 +- sympy/integrals/tests/test_integrals.py | 6 ++++ sympy/test_external/test_scipy.py | 46 +++++++++++++++++++++++++++++++ 4 files changed, 85 insertions(+), 12 deletions(-) create mode 100644 sympy/test_external/test_scipy.py diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py index a32d9b8..0314a03 100644 --- a/sympy/core/evalf.py +++ b/sympy/core/evalf.py @@ -280,7 +280,6 @@ def add_terms(terms, prec, target_prec): sum_accuracy = sum_exp + sum_bc - absolute_error r = normalize(sum_sign, sum_man, sum_exp, sum_bc, target_prec, round_nearest), sum_accuracy - #print "returning", to_str(r[0],50), r[1] return r def evalf_add(v, prec, options): @@ -730,16 +729,34 @@ def f(t): return result def evalf_integral(expr, prec, options): - workprec = prec - i = 0 - maxprec = options.get('maxprec', INF) - while 1: - result = do_integral(expr, workprec, options) - accuracy = complex_accuracy(result) - if accuracy >= prec or workprec >= maxprec: - return result - workprec += prec - max(-2**i, accuracy) - i += 1 + method = options.get('method', 'mpmath') + if method == 'scipy': + from scipy.integrate import quad + func = expr.args[0] + x, (xlow, xhigh) = expr.args[1][0] + from sympy import lambdify + f = lambdify(x, func, ["math"]) + # scipy only operates on the standard floating point precision: + eps = 1e-9 + re, im = quad(f, xlow, xhigh, epsabs=eps, epsrel=eps) + # chop imaginary parts that are below the precision + if abs(im) < eps: + im = 0 + re = from_float(re) + im = from_float(im) + return re, im, re[3], im[3] + elif method == 'mpmath': + workprec = prec + i = 0 + maxprec = options.get('maxprec', INF) + while 1: + result = do_integral(expr, workprec, options) + accuracy = complex_accuracy(result) + if accuracy >= prec or workprec >= maxprec: + return result + workprec += prec - max(-2**i, accuracy) + i += 1 + raise ValueError("Unknown method: %s" % method) def check_convergence(numer, denom, n): """ @@ -1003,6 +1020,10 @@ def Basic_evalf(x, n=15, **options): verbose=<bool> Print debug information (default=False) + method=<str> + mpmath ... uses mpmath for evaluation + scipy .... uses scipy where possible (e.g. evaluation of integrals) + """ if not evalf_table: _create_evalf_table() diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py index f9cc5f2..b6a74c6 100644 --- a/sympy/core/numbers.py +++ b/sympy/core/numbers.py @@ -642,7 +642,7 @@ def __ne__(self, other): if other.is_comparable and not isinstance(other, Rational): other = other.evalf() if isinstance(other, Number): if isinstance(other, Real): - return bool(not mlib.feq(self._as_mpf_val(other._prec), other._mpf_)) + return bool(not mlib.mpf_eq(self._as_mpf_val(other._prec), other._mpf_)) return bool(self.p!=other.p or self.q!=other.q) return True # Rational != non-Number diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py index 301090c..9a7bc63 100644 --- a/sympy/integrals/tests/test_integrals.py +++ b/sympy/integrals/tests/test_integrals.py @@ -339,3 +339,9 @@ def test_subs5(): def test_integration_variable(): raises(ValueError, "Integral(exp(-x**2), 3)") raises(ValueError, "Integral(exp(-x**2), (3, -oo, oo))") + +def test_integrals_evalf(): + i = Integral(exp(-x**2), (x, 0, oo)) + assert NS(i.evalf()) == "0.886226925452758" + assert NS(i.evalf(method="mpmath")) == "0.886226925452758" + raises(ValueError, "i.evalf(method='nonsense5')") diff --git a/sympy/test_external/test_scipy.py b/sympy/test_external/test_scipy.py new file mode 100644 index 0000000..3525f9b --- /dev/null +++ b/sympy/test_external/test_scipy.py @@ -0,0 +1,46 @@ +# This testfile tests SymPy <-> SciPy compatibility + +# Don't test any SymPy features here. Just pure interaction with SymPy. +# Always write regular SymPy tests for anything, that can be tested in pure +# Python (without numpy). Here we test everything, that a user may need when +# using SymPy with NumPy + +try: + import scipy +except ImportError: + #py.test will not execute any tests now + disabled = True + + +from sympy import (Integral, exp, oo, Symbol, sstr, sympify, raises) + +def NS(e, n=15, **options): + return sstr(sympify(e).evalf(n, **options), full_prec=True) + +def test_integrals_evalf_scipy1(): + x = Symbol("x") + i = Integral(exp(-x**2), (x, 0, oo)) + assert NS(i.evalf()) == "0.886226925452758" + assert NS(i.evalf(method="mpmath")) == "0.886226925452758" + assert NS(i.evalf(method="scipy")) == "0.886226925452758" + +def test_integrals_evalf_scipy2(): + x = Symbol("x") + i = Integral(exp(-x**2), (x, -oo, 0)) + assert NS(i.evalf()) == "0.886226925452758" + assert NS(i.evalf(method="mpmath")) == "0.886226925452758" + assert NS(i.evalf(method="scipy")) == "0.886226925452758" + +def test_integrals_evalf_scipy3(): + x = Symbol("x") + i = Integral(exp(-x**2), (x, -oo, oo)) + assert NS(i.evalf()) == "1.77245385090552" + assert NS(i.evalf(method="mpmath")) == "1.77245385090552" + assert NS(i.evalf(method="scipy")) == "1.77245385090552" + +def test_integrals_evalf_scipy4(): + x = Symbol("x") + i = Integral(exp(-x**2), (x, -0.1, 0.1)) + assert NS(i.evalf()) == "0.199335328580673" + assert NS(i.evalf(method="mpmath")) == "0.199335328580673" + assert NS(i.evalf(method="scipy")) == "0.199335328580673" -- 1.6.2.1 --~--~---------~--~----~------------~-------~--~----~ You received this message because you are subscribed to the Google Groups "sympy-patches" group. To post to this group, send email to sympy-patches@googlegroups.com To unsubscribe from this group, send email to sympy-patches+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sympy-patches?hl=en -~----------~----~----~----~------~----~------~--~---