mirror of
https://github.com/gonum/gonum.git
synced 2025-09-27 03:26:04 +08:00
5f0141ca4c
Changes made in dsp/fourier/internal/fftpack break the formatting used there, so these are reverted. There will be complaints in CI. [git-generate] gofmt -w . go generate gonum.org/v1/gonum/blas go generate gonum.org/v1/gonum/blas/gonum go generate gonum.org/v1/gonum/unit go generate gonum.org/v1/gonum/unit/constant go generate gonum.org/v1/gonum/graph/formats/dot go generate gonum.org/v1/gonum/graph/formats/rdf go generate gonum.org/v1/gonum/stat/card git checkout -- dsp/fourier/internal/fftpack
70 lines
1.5 KiB
Go
70 lines
1.5 KiB
Go
// 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]
|
|
}
|