integrate: Add Romberg integration

integrate/romberg.go

@@ -0,0 +1,65 @@
+// Copyright ©2019 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package integrate
+
+import (
+	"math"
+	"math/bits"
+)
+
+// Romberg returns an approximate value of the integral
+//  \int_a^b f(x)dx
+// computed using the Romberg's method. The function f is given
+// as a slice of equally-spaced samples, that is,
+//  f[i] = f(a + i*dx)
+// and dx is the spacing between the samples.
+//
+// The length of f must be 2^k + 1, where k is a positive integer,
+// and dx must be positive.
+//
+// See https://en.wikipedia.org/wiki/Romberg%27s_method for a description of
+// the algorithm
+func Romberg(f []float64, dx float64) float64 {
+	if len(f) < 3 {
+		panic("integral: invalid slice length: must be at least 3")
+	}
+
+	if dx <= 0 {
+		panic("integral: invalid spacing: must be larger than 0")
+	}
+
+	n := len(f) - 1
+	k := bits.Len(uint(n - 1))
+
+	if len(f) != 1<<uint(k)+1 {
+		panic("integral: invalid slice length: must be 2^k + 1")
+	}
+
+	work := make([]float64, 2*(k+1))
+	prev := work[:k+1]
+	curr := work[k+1:]
+
+	h := dx * float64(n)
+	prev[0] = (f[0] + f[len(f)-1]) * 0.5 * h
+
+	step := n
+	for i := 1; i <= k; i++ {
+		h /= 2
+		step /= 2
+		var estimate float64
+		for j := 0; j < n/2; j += step {
+			estimate += f[2*j+step]
+		}
+
+		curr[0] = estimate*h + 0.5*prev[0]
+		for j := 1; j <= i; j++ {
+			factor := math.Pow(4, float64(j))
+			curr[j] = (factor*curr[j-1] - prev[j-1]) / (factor - 1)
+		}
+
+		prev, curr = curr, prev
+	}
+	return prev[k]
+}

integrate/romberg_test.go

@@ -0,0 +1,86 @@
+// Copyright ©2019 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package integrate
+
+import (
+	"math"
+	"testing"
+
+	"gonum.org/v1/gonum/floats"
+)
+
+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
+	}{
+		{
+			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,
+		},
+	} {
+		n := len(test.x)
+		y := make([]float64, n)
+		for i, v := range test.x {
+			y[i] = test.f(v)
+		}
+
+		dx := (test.x[n-1] - test.x[0]) / float64(n-1)
+		v := Romberg(y, dx)
+		diff := math.Abs(v - test.want)
+		if diff > test.tol {
+			t.Errorf("test #%d: got=%v want=%v diff=%v\n", i, v, test.want, diff)
+		}
+	}
+}