gonum.org/v1/gonum@v0.14.0/integrate/romberg_test.go (about)

     1  // Copyright ©2019 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  	"gonum.org/v1/gonum/floats"
    12  	"gonum.org/v1/gonum/integrate/testquad"
    13  )
    14  
    15  func TestRomberg(t *testing.T) {
    16  	t.Parallel()
    17  	for i, test := range []struct {
    18  		integral testquad.Integral
    19  		n        int
    20  		tol      float64
    21  	}{
    22  		{integral: testquad.Constant(0), n: 3, tol: 0},
    23  		{integral: testquad.Constant(0), n: 1<<5 + 1, tol: 0},
    24  		{integral: testquad.Poly(0), n: 3, tol: 1e-14},
    25  		{integral: testquad.Poly(0), n: 1<<5 + 1, tol: 1e-14},
    26  		{integral: testquad.Poly(1), n: 3, tol: 1e-14},
    27  		{integral: testquad.Poly(1), n: 1<<5 + 1, tol: 1e-14},
    28  		{integral: testquad.Poly(2), n: 3, tol: 1e-14},
    29  		{integral: testquad.Poly(2), n: 1<<5 + 1, tol: 1e-14},
    30  		{integral: testquad.Poly(3), n: 3, tol: 1e-14},
    31  		{integral: testquad.Poly(3), n: 1<<5 + 1, tol: 1e-14},
    32  		{integral: testquad.Poly(4), n: 5, tol: 1e-14},
    33  		{integral: testquad.Poly(4), n: 1<<5 + 1, tol: 1e-14},
    34  		{integral: testquad.Poly(5), n: 5, tol: 1e-14},
    35  		{integral: testquad.Poly(5), n: 1<<5 + 1, tol: 1e-14},
    36  		{integral: testquad.Sin(), n: 1<<3 + 1, tol: 1e-10},
    37  		{integral: testquad.Sin(), n: 1<<5 + 1, tol: 1e-14},
    38  		{integral: testquad.XExpMinusX(), n: 1<<3 + 1, tol: 1e-9},
    39  		{integral: testquad.XExpMinusX(), n: 1<<5 + 1, tol: 1e-14},
    40  		{integral: testquad.Sqrt(), n: 1<<10 + 1, tol: 1e-5},
    41  		{integral: testquad.ExpOverX2Plus1(), n: 1<<4 + 1, tol: 1e-7},
    42  		{integral: testquad.ExpOverX2Plus1(), n: 1<<6 + 1, tol: 1e-14},
    43  	} {
    44  		n := test.n
    45  		a := test.integral.A
    46  		b := test.integral.B
    47  
    48  		x := make([]float64, n)
    49  		floats.Span(x, a, b)
    50  
    51  		y := make([]float64, n)
    52  		for i, xi := range x {
    53  			y[i] = test.integral.F(xi)
    54  		}
    55  
    56  		dx := (b - a) / float64(n-1)
    57  		got := Romberg(y, dx)
    58  
    59  		want := test.integral.Value
    60  		diff := math.Abs(got - want)
    61  		if diff > test.tol {
    62  			t.Errorf("Test #%d: %v, n=%v: unexpected result; got=%v want=%v diff=%v",
    63  				i, test.integral.Name, n, got, want, diff)
    64  		}
    65  	}
    66  }