integrate: rework tests to use testquad integrals
diff --git a/integrate/romberg_test.go b/integrate/romberg_test.go
index 58ac4c3..00c34ad 100644
--- a/integrate/romberg_test.go
+++ b/integrate/romberg_test.go
@@ -9,78 +9,57 @@
"testing"
"gonum.org/v1/gonum/floats"
+ "gonum.org/v1/gonum/integrate/testquad"
)
func TestRomberg(t *testing.T) {
- const n = 1<<8 + 1
- x := floats.Span(make([]float64, n), 0, 1)
-
for i, test := range []struct {
- x []float64
- f func(x float64) float64
- want float64
- tol float64
+ integral testquad.Integral
+ n int
+ tol float64
}{
- {
- x: x,
- f: func(x float64) float64 { return x },
- want: 0.5,
- tol: 1e-12,
- },
- {
- x: x,
- f: func(x float64) float64 { return x * x },
- want: 1.0 / 3.0,
- tol: 1e-12,
- },
- {
- x: x,
- f: func(x float64) float64 { return x * x * x },
- want: 1.0 / 4.0,
- tol: 1e-12,
- },
- {
- x: x,
- f: func(x float64) float64 { return math.Sqrt(x) },
- want: 2.0 / 3.0,
- tol: 1e-4,
- },
- {
- x: x,
- f: func(x float64) float64 { return math.Sin(math.Pi * x) },
- want: 2.0 / math.Pi,
- tol: 1e-12,
- },
- {
- x: floats.Span(make([]float64, 3), 0, 1),
- f: func(x float64) float64 { return x * x },
- want: 1.0 / 3.0,
- tol: 1e-12,
- },
- {
- x: floats.Span(make([]float64, 3), 0, 1),
- f: func(x float64) float64 { return x * x * x },
- want: 1.0 / 4.0,
- tol: 1e-12,
- },
- {
- x: x,
- f: func(x float64) float64 { return x * math.Exp(-x) },
- want: (math.Exp(1) - 2) / math.Exp(1),
- tol: 1e-12,
- },
+ {integral: testquad.Constant(0), n: 3, tol: 0},
+ {integral: testquad.Constant(0), n: 1<<5 + 1, tol: 0},
+ {integral: testquad.Poly(0), n: 3, tol: 1e-14},
+ {integral: testquad.Poly(0), n: 1<<5 + 1, tol: 1e-14},
+ {integral: testquad.Poly(1), n: 3, tol: 1e-14},
+ {integral: testquad.Poly(1), n: 1<<5 + 1, tol: 1e-14},
+ {integral: testquad.Poly(2), n: 3, tol: 1e-14},
+ {integral: testquad.Poly(2), n: 1<<5 + 1, tol: 1e-14},
+ {integral: testquad.Poly(3), n: 3, tol: 1e-14},
+ {integral: testquad.Poly(3), n: 1<<5 + 1, tol: 1e-14},
+ {integral: testquad.Poly(4), n: 5, tol: 1e-14},
+ {integral: testquad.Poly(4), n: 1<<5 + 1, tol: 1e-14},
+ {integral: testquad.Poly(5), n: 5, tol: 1e-14},
+ {integral: testquad.Poly(5), n: 1<<5 + 1, tol: 1e-14},
+ {integral: testquad.Sin(), n: 1<<3 + 1, tol: 1e-10},
+ {integral: testquad.Sin(), n: 1<<5 + 1, tol: 1e-14},
+ {integral: testquad.XExpMinusX(), n: 1<<3 + 1, tol: 1e-9},
+ {integral: testquad.XExpMinusX(), n: 1<<5 + 1, tol: 1e-14},
+ {integral: testquad.Sqrt(), n: 1<<10 + 1, tol: 1e-5},
+ {integral: testquad.ExpOverX2Plus1(), n: 1<<4 + 1, tol: 1e-7},
+ {integral: testquad.ExpOverX2Plus1(), n: 1<<6 + 1, tol: 1e-14},
} {
- n := len(test.x)
+ n := test.n
+ a := test.integral.A
+ b := test.integral.B
+
+ x := make([]float64, n)
+ floats.Span(x, a, b)
+
y := make([]float64, n)
- for i, v := range test.x {
- y[i] = test.f(v)
+ for i, xi := range x {
+ y[i] = test.integral.F(xi)
}
- dx := (test.x[n-1] - test.x[0]) / float64(n-1)
- v := Romberg(y, dx)
- diff := math.Abs(v - test.want)
+ dx := (b - a) / float64(n-1)
+ got := Romberg(y, dx)
+
+ want := test.integral.Value
+ diff := math.Abs(got - want)
if diff > test.tol {
- t.Errorf("test #%d: got=%v want=%v diff=%v\n", i, v, test.want, diff)
+ t.Errorf("Test #%d: %v, n=%v: unexpected result; got=%v want=%v diff=%v",
+ i, test.integral.Name, n, got, want, diff)
}
}
}
diff --git a/integrate/simpsons_test.go b/integrate/simpsons_test.go
index cdaff0c..cb6dc54 100644
--- a/integrate/simpsons_test.go
+++ b/integrate/simpsons_test.go
@@ -8,165 +8,51 @@
"math"
"testing"
- "gonum.org/v1/gonum/floats"
+ "golang.org/x/exp/rand"
+
+ "gonum.org/v1/gonum/integrate/testquad"
)
func TestSimpsons(t *testing.T) {
- x := floats.Span(make([]float64, 1e6), 0, 1)
+ rnd := rand.New(rand.NewSource(1))
for i, test := range []struct {
- x []float64
- f func(x float64) float64
- want float64
- tol float64
+ integral testquad.Integral
+ n int
+ tol float64
}{
- {
- x: floats.Span(make([]float64, 1e6), 0, 1),
- f: func(x float64) float64 { return math.Pi },
- want: math.Pi,
- tol: 1e-10,
- },
- {
- x: floats.Span(make([]float64, 1e6), 0, 1),
- f: func(x float64) float64 { return 1.0 },
- want: 1,
- tol: 1e-10,
- },
- {
- x: floats.Span(make([]float64, 3), 0, 2),
- f: func(x float64) float64 { return 2*x + 0.5 },
- want: 5,
- tol: 1e-12,
- },
- {
- x: floats.Span(make([]float64, 10), 0, 2),
- f: func(x float64) float64 { return 2*x + 0.5 },
- want: 5,
- tol: 1e-12,
- },
- {
- x: floats.Span(make([]float64, 1e6), 0, 2),
- f: func(x float64) float64 { return 2*x + 0.5 },
- want: 5,
- tol: 1e-12,
- },
- {
- x: x,
- f: func(x float64) float64 { return x },
- want: 0.5,
- tol: 1e-12,
- },
- {
- x: floats.Span(make([]float64, 1e6), -1, 1),
- f: func(x float64) float64 { return x },
- want: 0,
- tol: 1e-12,
- },
- {
- x: x,
- f: func(x float64) float64 { return x + 10 },
- want: 10.5,
- tol: 1e-12,
- },
- {
- x: x,
- f: func(x float64) float64 { return 3*x*x + 10 },
- want: 11,
- tol: 1e-12,
- },
- {
- x: x,
- f: func(x float64) float64 { return math.Exp(x) },
- want: 1.7182818284591876,
- tol: 1e-12,
- },
- {
- x: floats.Span(make([]float64, 1e6), 0, math.Pi),
- f: func(x float64) float64 { return math.Cos(x) },
- want: 0,
- tol: 1e-12,
- },
- {
- x: floats.Span(make([]float64, 1e6), 0, 2*math.Pi),
- f: func(x float64) float64 { return math.Cos(x) },
- want: 0,
- tol: 1e-12,
- },
- {
- x: floats.Span(make([]float64, 1e6*10), 0, math.Pi),
- f: func(x float64) float64 { return math.Sin(x) },
- want: 2,
- tol: 1e-12,
- },
- {
- x: floats.Span(make([]float64, 1e6*10), 0, 0.5*math.Pi),
- f: func(x float64) float64 { return math.Sin(x) },
- want: 1,
- tol: 1e-12,
- },
- {
- x: floats.Span(make([]float64, 1e6), 0, 2*math.Pi),
- f: func(x float64) float64 { return math.Sin(x) },
- want: 0,
- tol: 1e-12,
- },
- {
- x: join(floats.Span(make([]float64, 3), 0, math.Pi/3),
- []float64{4 * math.Pi / 10},
- floats.Span(make([]float64, 4), 3*math.Pi/4, math.Pi)),
- f: func(x float64) float64 { return math.Sin(x) },
- want: 2,
- tol: 1e-2,
- },
- {
- x: join(floats.Span(make([]float64, 30), 0, 5*math.Pi/16),
- floats.Span(make([]float64, 100), 3*math.Pi/8, math.Pi)),
- f: func(x float64) float64 { return math.Sin(x) },
- want: 2,
- tol: 1e-4,
- },
- {
- x: join(floats.Span(make([]float64, 1e5), 0, 15),
- floats.Span(make([]float64, 1e5), 23, 40),
- floats.Span(make([]float64, 1e5), 50, 80),
- floats.Span(make([]float64, 1e5), 90, 100)),
- f: func(x float64) float64 { return 2 * x },
- want: 10000,
- tol: 1e-9,
- },
- {
- x: []float64{0, 1, 2},
- f: func(x float64) float64 { return 2 * x },
- want: 4,
- tol: 1e-12,
- },
- {
- x: []float64{0, 0.2, 2},
- f: func(x float64) float64 { return 2 * x },
- want: 4,
- tol: 1e-12,
- },
- {
- x: []float64{0, 1.89, 2},
- f: func(x float64) float64 { return 2 * x },
- want: 4,
- tol: 1e-12,
- },
+ {integral: testquad.Constant(0), n: 3, tol: 0},
+ {integral: testquad.Constant(0), n: 10, tol: 0},
+ {integral: testquad.Poly(0), n: 3, tol: 1e-14},
+ {integral: testquad.Poly(0), n: 10, tol: 1e-14},
+ {integral: testquad.Poly(1), n: 3, tol: 1e-14},
+ {integral: testquad.Poly(1), n: 10, tol: 1e-14},
+ {integral: testquad.Poly(2), n: 3, tol: 1e-14},
+ {integral: testquad.Poly(2), n: 10, tol: 1e-14},
+ {integral: testquad.Poly(3), n: 1e3, tol: 1e-8},
+ {integral: testquad.Poly(4), n: 1e3, tol: 1e-8},
+ {integral: testquad.Poly(5), n: 1e3, tol: 1e-7},
+ {integral: testquad.Sin(), n: 1e2, tol: 1e-8},
+ {integral: testquad.XExpMinusX(), n: 1e2, tol: 1e-8},
+ {integral: testquad.Sqrt(), n: 1e4, tol: 1e-6},
+ {integral: testquad.ExpOverX2Plus1(), n: 1e2, tol: 1e-7},
} {
- y := make([]float64, len(test.x))
- for i, v := range test.x {
- y[i] = test.f(v)
- }
- v := Simpsons(test.x, y)
- if !floats.EqualWithinAbs(v, test.want, test.tol) {
- t.Errorf("test #%d: got=%v want=%f\n", i, v, float64(test.want))
- }
- }
-}
+ n := test.n
+ a := test.integral.A
+ b := test.integral.B
-func join(slices ...[]float64) []float64 {
- var c []float64
- for _, s := range slices {
- c = append(c, s...)
+ x := jitterSpan(n, a, b, rnd)
+ y := make([]float64, n)
+ for i, xi := range x {
+ y[i] = test.integral.F(xi)
+ }
+
+ got := Simpsons(x, y)
+
+ want := test.integral.Value
+ diff := math.Abs(got - want)
+ if diff > test.tol {
+ t.Errorf("Test #%d: %v, n=%v: unexpected result; got=%v want=%v diff=%v",
+ i, test.integral.Name, n, got, want, diff)
+ }
}
- return c
}
diff --git a/integrate/trapezoidal_test.go b/integrate/trapezoidal_test.go
index 1a4451e..68b7379 100644
--- a/integrate/trapezoidal_test.go
+++ b/integrate/trapezoidal_test.go
@@ -8,75 +8,66 @@
"math"
"testing"
- "gonum.org/v1/gonum/floats"
+ "golang.org/x/exp/rand"
+
+ "gonum.org/v1/gonum/integrate/testquad"
)
func TestTrapezoidal(t *testing.T) {
- const N = 1e6
- x := floats.Span(make([]float64, N), 0, 1)
+ rnd := rand.New(rand.NewSource(1))
for i, test := range []struct {
- x []float64
- f func(x float64) float64
- want float64
+ integral testquad.Integral
+ n int
+ tol float64
}{
- {
- x: x,
- f: func(x float64) float64 { return x },
- want: 0.5,
- },
- {
- x: floats.Span(make([]float64, N), -1, 1),
- f: func(x float64) float64 { return x },
- want: 0,
- },
- {
- x: x,
- f: func(x float64) float64 { return x + 10 },
- want: 10.5,
- },
- {
- x: x,
- f: func(x float64) float64 { return 3*x*x + 10 },
- want: 11,
- },
- {
- x: x,
- f: func(x float64) float64 { return math.Exp(x) },
- want: 1.7182818284591876,
- },
- {
- x: floats.Span(make([]float64, N), 0, math.Pi),
- f: func(x float64) float64 { return math.Cos(x) },
- want: 0,
- },
- {
- x: floats.Span(make([]float64, N), 0, 2*math.Pi),
- f: func(x float64) float64 { return math.Cos(x) },
- want: 0,
- },
- {
- x: floats.Span(make([]float64, N*10), 0, math.Pi),
- f: func(x float64) float64 { return math.Sin(x) },
- want: 2,
- },
- {
- x: floats.Span(make([]float64, N*10), 0, 0.5*math.Pi),
- f: func(x float64) float64 { return math.Sin(x) },
- want: 1,
- },
- {
- x: floats.Span(make([]float64, N), 0, 2*math.Pi),
- f: func(x float64) float64 { return math.Sin(x) },
- want: 0,
- },
+ {integral: testquad.Constant(0), n: 2, tol: 0},
+ {integral: testquad.Constant(0), n: 10, tol: 0},
+ {integral: testquad.Poly(0), n: 2, tol: 1e-14},
+ {integral: testquad.Poly(0), n: 10, tol: 1e-14},
+ {integral: testquad.Poly(1), n: 2, tol: 1e-14},
+ {integral: testquad.Poly(1), n: 10, tol: 1e-14},
+ {integral: testquad.Poly(2), n: 1e5, tol: 1e-8},
+ {integral: testquad.Poly(3), n: 1e5, tol: 1e-8},
+ {integral: testquad.Poly(4), n: 1e5, tol: 1e-7},
+ {integral: testquad.Poly(5), n: 1e5, tol: 1e-7},
+ {integral: testquad.Sin(), n: 1e5, tol: 1e-11},
+ {integral: testquad.XExpMinusX(), n: 1e5, tol: 1e-10},
+ {integral: testquad.Sqrt(), n: 1e5, tol: 1e-8},
+ {integral: testquad.ExpOverX2Plus1(), n: 1e5, tol: 1e-10},
} {
- y := make([]float64, len(test.x))
- for i, v := range test.x {
- y[i] = test.f(v)
+ n := test.n
+ a := test.integral.A
+ b := test.integral.B
+
+ x := jitterSpan(n, a, b, rnd)
+ y := make([]float64, n)
+ for i, xi := range x {
+ y[i] = test.integral.F(xi)
}
- v := Trapezoidal(test.x, y)
- if !floats.EqualWithinAbs(v, test.want, 1e-12) {
- t.Errorf("test #%d: got=%v want=%v\n", i, v, test.want)
+
+ got := Trapezoidal(x, y)
+
+ want := test.integral.Value
+ diff := math.Abs(got - want)
+ if diff > test.tol {
+ t.Errorf("Test #%d: %v, n=%v: unexpected result; got=%v want=%v diff=%v",
+ i, test.integral.Name, n, got, want, diff)
}
}
}
+
+func jitterSpan(n int, a, b float64, rnd *rand.Rand) []float64 {
+ dx := (b - a) / float64(n-1)
+ x := make([]float64, n)
+ x[0] = a
+ for i := 1; i < n-1; i++ {
+ // Set x[i] to its regular location.
+ x[i] = a + float64(i)*dx
+ // Generate a random number in [-1,1).
+ jitter := 2*rnd.Float64() - 1
+ // Jitter x[i] without crossing over its neighbors.
+ x[i] += 0.4 * jitter * dx
+ }
+ x[n-1] = b
+ return x
+}
