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
-~----------~----~----~----~------~----~------~--~---
