diff options
Diffstat (limited to 'vquad/test/test_core.py')
| -rw-r--r-- | vquad/test/test_core.py | 207 |
1 files changed, 160 insertions, 47 deletions
diff --git a/vquad/test/test_core.py b/vquad/test/test_core.py index 46fd823..f5fce6d 100644 --- a/vquad/test/test_core.py +++ b/vquad/test/test_core.py @@ -3,41 +3,170 @@ # This file is part of Vquad. It is subject to the license terms in the file # LICENSE.rst found in the top-level directory of this distribution. +from math import exp, log, log10 import numpy as np from numpy.testing import assert_allclose from pytest import raises from .. import core -def f0(x): - return x * np.sin(1/x) * np.sqrt(abs(1 - x)) +def f1(x): + return np.exp(x) + +def f2(x): + return x >= 0.3 + +def f3(x): + return np.sqrt(x) + +def f4(x): + return 23/25 * np.cosh(x) - np.cos(x) + +def f5(x): + xx = x*x + return 1 / ( xx * (xx + 1) + 0.9) + +def f6(x): + return x * np.sqrt(x) def f7(x): - return x**-0.5 + return 1 / np.sqrt(x) +def f8(x): + xx = x*x + return 1 / (1 + xx*xx) -def f24(x): - return np.floor(np.exp(x)) +def f9(x): + return 2 / (2 + np.sin(10 * np.pi * x)) + +def f10(x): + return 1 / (1 + x) + +def f11(x): + return 1 / (1 + np.exp(x)) + +def f12(x): + return x / (np.exp(x) - 1) + +def f13(x): + return np.sin(100 * np.pi * x) / (np.pi * x) + +def f14(x): + return np.sqrt(50) * np.exp(-50 * np.pi * x * x) +def f15(x): + return 25 * np.exp(-25 * x) + +def f16(x): + return 50 / np.pi * (2500 * x*x + 1) + +def f17(x): + t = 50 * np.pi * x + t = np.sin(t) / t + return 50 * t * t + +def f18(x): + return np.cos( np.cos(x) + 3 * np.sin(x) + 2 * np.cos(2*x) + + 3 * np.sin(2*x) + 3 * np.cos(3*x) ) + +def f19(x): + return np.log(x) + +def f20(x): + return 1 / (1.005 + x*x) def f21(x): - y = 0 - for i in range(1, 4): - y += 1 / np.cosh(20**i * (x - 2 * i / 10)) - return y + return sum(1 / np.cosh(20**i * (x - 2*i/10)) for i in range(1, 4)) +def f22(x): + return 4 * np.pi**2 * x * np.sin(20 * np.pi * x) * np.cos(2 * np.pi * x) -def f63(x, alpha, beta): - return abs(x - beta) ** alpha +def f23(x): + t = 230 * x - 30 + return 1 / (1 + t*t) +def f24(x): + return np.floor(np.exp(x)) -def F63(x, alpha, beta): - return (x - beta) * abs(x - beta) ** alpha / (alpha + 1) +def f25(x): + return ((x + 1) * (x < 1) + + (3 - x) * ((1 <= x) & (x <= 3)) + + 2 * (x > 3)) + +battery = [ + (f1, (0, 1), 1.7182818284590452354), + (f2, (0, 1), 0.7), + (f3, (0, 1), 2/3), + (f4, (-1, 1), 0.4794282266888016674), + (f5, (-1, 1), 1.5822329637296729331), + (f6, (0, 1), 0.4), + (f7, (0, 1), 2), + (f8, (0, 1), 0.86697298733991103757), + (f9, (0, 1), 1.1547005383792515290), + (f10, (0, 1), 0.69314718055994530942), + (f11, (0, 1), 0.3798854930417224753), + (f12, (0, 1), 0.77750463411224827640), + (f13, (0, 1), 0.49898680869304550249), + (f14, (0, 10), 0.5), + (f15, (0, 10), 1), + (f16, (0, 10), 0.13263071079267703209e+08), + (f17, (0, 1), 0.49898680869304550249), + (f18, (0, np.pi), 0.83867634269442961454), + (f19, (0, 1), -1), + (f20, (-1, 1), 1.5643964440690497731), + (f21, (0, 1), 0.16349494301863722618), + (f22, (0, 1), -0.63466518254339257343), + (f23, (0, 1), 0.013492485649467772692), + (f24, (0, 3), 17.664383539246514971), + (f25, (0, 5), 7.5), +] + +def test_battery(): + """Test and gather statistics on a battery of difficult integrands + + The battery is used as listed on page 168 of Pedro Gonnet's thesis. + + Reference: + W. Gander and W. Gautschi. Adaptive quadrature — revisited. + Technical Report 306, Department of Computer Science, ETH + Zurich, Switzerland, 1998. + """ + def ff(x): + nonlocal neval + neval += len(x) + return f(x) + old_settings = np.seterr(all='ignore') -def fdiv(x): - return abs(x - 0.987654321) ** -1.1 + rtols = [1e-3, 1e-6, 1e-9, 1e-12] + sum_log_neval = 0 + sum_extra_digits = 0 + sum_nonrep_digits = 0 + for f, (a, b), exact in battery: + for rtol in rtols: + neval = 0 + igral, err = core.vquad(ff, a, b, rtol) + sum_log_neval += log(neval) + extra_digits = min(16, log10(rtol / abs((igral - exact) / exact))) + sum_extra_digits += extra_digits + nonrep_digits = min(16, log10(err / abs(igral - exact))) + sum_nonrep_digits += nonrep_digits + if not (f is f21 and rtol > 1e-11): + assert extra_digits > 0 + assert nonrep_digits > 0 + + n = len(rtols) * len(battery) + print("geom. mean of number of evaluations:", + round(exp(sum_log_neval / n), 2), "(cquad: 237.92)") + print("mean non-requested significant digits:", + round(sum_extra_digits / n, 3), + "(cquad: 5.989)") + print("mean non-reported significant digits:", + round(sum_nonrep_digits / n, 3), + "(cquad: 4.056)") + + np.seterr(**old_settings) def f_one_with_nan(x): @@ -47,6 +176,14 @@ def f_one_with_nan(x): return result +def test_one_with_nan(): + rtol = 1e-12 + exact = 2 + igral, err = core.vquad(f_one_with_nan, -1, 1, rtol) + assert_allclose(igral, exact, rtol) + assert err / exact < rtol + + def test_coeffs(level=3): vals = np.abs(core.tbls.nodes[level]) vals[1::2] = np.nan @@ -61,38 +198,6 @@ def test_coeffs(level=3): assert_allclose(vals_from_c, vals, rtol=0, atol=1e-15) -def test_integration(): - old_settings = np.seterr(all='ignore') - - igral, err = core.vquad(f0, 0, 3, 1e-5) - assert_allclose(igral, 1.98194117954329, 1e-15) - assert_allclose(err, 1.9563545589988155e-05, 1e-10) - - igral, err = core.vquad(f7, 0, 1, 1e-6) - assert_allclose(igral, 1.9999998579359648, 1e-15) - assert_allclose(err, 1.8561437334964041e-06, 1e-10) - - igral, err = core.vquad(f24, 0, 3, 1e-3) - assert_allclose(igral, 17.664696186312934, 1e-15) - assert_allclose(err, 0.017602618074957457, 1e-10) - - igral, err = core.vquad(f21, 0, 1, 1e-3) - assert_allclose(igral, 0.16310022131213361, 1e-15) - assert_allclose(err, 0.00011848806384952786, 1e-10) - - igral, err = core.vquad(f_one_with_nan, -1, 1, 1e-12) - assert_allclose(igral, 2, 1e-15) - assert_allclose(err, 2.4237853822937613e-15, 1e-7) - - try: - igral, err = core.vquad(fdiv, 0, 1, 1e-6) - except core.DivergentIntegralError as e: - assert e.igral == np.inf - assert e.err is None - - np.seterr(**old_settings) - - def test_interpolation(): vquad = core.Vquad(f24, 0, 3) vquad.improve_until(1e-6) @@ -108,6 +213,14 @@ def test_interpolation(): raises(ValueError, vquad, x) +def f63(x, alpha, beta): + return abs(x - beta) ** alpha + + +def F63(x, alpha, beta): + return (x - beta) * abs(x - beta) ** alpha / (alpha + 1) + + def test_analytic(n=200): def f(x): return f63(x, alpha, beta) |