1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
|
// 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[n]) * 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]
}
|
