Repository: gonum/gonum
Pull request #1129 (2019-10-13 14:19:03)

Integrate: Add Romberg integration
This PR implements the Romberg integration scheme. There are two methods: `Romberg` and `RombergHL`. The latter is simply a convenience function, which can be useful for cases where users don't want to allocate a `work` array themselves.


davidkleiven: Good points, I think all your suggestions have been addressed.
davidkleiven: @vladimir-ch 
> Is h supposed to be equal to b - a?

Yes, it is the case. And it is fixed, together with all your other suggestions.

davidkleiven: Thanks for your thorogh checking. Many of the things I should have spotted myself of course, but I missed them. @vladimir-ch Your latest suggestion should now be fixed.
kortschak: There are two approaches. Either squash all these commits into one and add another commit, "A+C: add David Kleiven" and then we merge is as a single PR. Otherwise, we can squash this in the merge and then another PR to add you to A+C.

Patch:
diff --git a/AUTHORS b/AUTHORS
index b407b2ae..132a341d 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -21,6 +21,7 @@ Clayton Northey <clayton.northey@gmail.com>
 Dan Kortschak <dan.kortschak@adelaide.edu.au> <dan@kortschak.io>
 Daniel Fireman <danielfireman@gmail.com>
 Dario Heinisch <dario.heinisch@gmail.com>
+David Kleiven <davidkleiven446@gmail.com>
 David Samborski <bloggingarrow@gmail.com>
 Davor Kapsa <davor.kapsa@gmail.com>
 DeepMind Technologies
diff --git a/CONTRIBUTORS b/CONTRIBUTORS
index 1b6e3414..e90db65c 100644
--- a/CONTRIBUTORS
+++ b/CONTRIBUTORS
@@ -29,6 +29,7 @@ Clayton Northey <clayton.northey@gmail.com>
 Dan Kortschak <dan.kortschak@adelaide.edu.au> <dan@kortschak.io>
 Daniel Fireman <danielfireman@gmail.com>
 Dario Heinisch <dario.heinisch@gmail.com>
+David Kleiven <davidkleiven446@gmail.com>
 David Samborski <bloggingarrow@gmail.com>
 Davor Kapsa <davor.kapsa@gmail.com>
 Delaney Gillilan <delaneygillilan@gmail.com>
diff --git a/integrate/romberg.go b/integrate/romberg.go
new file mode 100644
index 00000000..585b5868
--- /dev/null
+++ b/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]
+}


Test patch:
diff --git a/integrate/romberg_test.go b/integrate/romberg_test.go
new file mode 100644
index 00000000..58ac4c32
--- /dev/null
+++ b/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)
+		}
+	}
+}