github.com/gopherd/gonum@v0.0.4/integrate/trapezoidal_test.go (about)

     1  // Copyright ©2016 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package integrate
     6  
     7  import (
     8  	"math"
     9  	"testing"
    10  
    11  	"math/rand"
    12  
    13  	"github.com/gopherd/gonum/integrate/testquad"
    14  )
    15  
    16  func TestTrapezoidal(t *testing.T) {
    17  	t.Parallel()
    18  	rnd := rand.New(rand.NewSource(1))
    19  	for i, test := range []struct {
    20  		integral testquad.Integral
    21  		n        int
    22  		tol      float64
    23  	}{
    24  		{integral: testquad.Constant(0), n: 2, tol: 0},
    25  		{integral: testquad.Constant(0), n: 10, tol: 0},
    26  		{integral: testquad.Poly(0), n: 2, tol: 1e-14},
    27  		{integral: testquad.Poly(0), n: 10, tol: 1e-14},
    28  		{integral: testquad.Poly(1), n: 2, tol: 1e-14},
    29  		{integral: testquad.Poly(1), n: 10, tol: 1e-14},
    30  		{integral: testquad.Poly(2), n: 1e5, tol: 1e-8},
    31  		{integral: testquad.Poly(3), n: 1e5, tol: 1e-8},
    32  		{integral: testquad.Poly(4), n: 1e5, tol: 1e-7},
    33  		{integral: testquad.Poly(5), n: 1e5, tol: 1e-7},
    34  		{integral: testquad.Sin(), n: 1e5, tol: 1e-11},
    35  		{integral: testquad.XExpMinusX(), n: 1e5, tol: 1e-10},
    36  		{integral: testquad.Sqrt(), n: 1e5, tol: 1e-8},
    37  		{integral: testquad.ExpOverX2Plus1(), n: 1e5, tol: 1e-10},
    38  	} {
    39  		n := test.n
    40  		a := test.integral.A
    41  		b := test.integral.B
    42  
    43  		x := jitterSpan(n, a, b, rnd)
    44  		y := make([]float64, n)
    45  		for i, xi := range x {
    46  			y[i] = test.integral.F(xi)
    47  		}
    48  
    49  		got := Trapezoidal(x, y)
    50  
    51  		want := test.integral.Value
    52  		diff := math.Abs(got - want)
    53  		if diff > test.tol {
    54  			t.Errorf("Test #%d: %v, n=%v: unexpected result; got=%v want=%v diff=%v",
    55  				i, test.integral.Name, n, got, want, diff)
    56  		}
    57  	}
    58  }
    59  
    60  func jitterSpan(n int, a, b float64, rnd *rand.Rand) []float64 {
    61  	dx := (b - a) / float64(n-1)
    62  	x := make([]float64, n)
    63  	x[0] = a
    64  	for i := 1; i < n-1; i++ {
    65  		// Set x[i] to its regular location.
    66  		x[i] = a + float64(i)*dx
    67  		// Generate a random number in [-1,1).
    68  		jitter := 2*rnd.Float64() - 1
    69  		// Jitter x[i] without crossing over its neighbors.
    70  		x[i] += 0.4 * jitter * dx
    71  	}
    72  	x[n-1] = b
    73  	return x
    74  }