github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/integrate/romberg.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  	"math/bits"
    10  )
    11  
    12  // Romberg returns an approximate value of the integral
    13  //  \int_a^b f(x)dx
    14  // computed using the Romberg's method. The function f is given
    15  // as a slice of equally-spaced samples, that is,
    16  //  f[i] = f(a + i*dx)
    17  // and dx is the spacing between the samples.
    18  //
    19  // The length of f must be 2^k + 1, where k is a positive integer,
    20  // and dx must be positive.
    21  //
    22  // See https://en.wikipedia.org/wiki/Romberg%27s_method for a description of
    23  // the algorithm.
    24  func Romberg(f []float64, dx float64) float64 {
    25  	if len(f) < 3 {
    26  		panic("integral: invalid slice length: must be at least 3")
    27  	}
    28  
    29  	if dx <= 0 {
    30  		panic("integral: invalid spacing: must be larger than 0")
    31  	}
    32  
    33  	n := len(f) - 1
    34  	k := bits.Len(uint(n - 1))
    35  
    36  	if len(f) != 1<<uint(k)+1 {
    37  		panic("integral: invalid slice length: must be 2^k + 1")
    38  	}
    39  
    40  	work := make([]float64, 2*(k+1))
    41  	prev := work[:k+1]
    42  	curr := work[k+1:]
    43  
    44  	h := dx * float64(n)
    45  	prev[0] = (f[0] + f[n]) * 0.5 * h
    46  
    47  	step := n
    48  	for i := 1; i <= k; i++ {
    49  		h /= 2
    50  		step /= 2
    51  		var estimate float64
    52  		for j := 0; j < n/2; j += step {
    53  			estimate += f[2*j+step]
    54  		}
    55  
    56  		curr[0] = estimate*h + 0.5*prev[0]
    57  		for j := 1; j <= i; j++ {
    58  			factor := math.Pow(4, float64(j))
    59  			curr[j] = (factor*curr[j-1] - prev[j-1]) / (factor - 1)
    60  		}
    61  
    62  		prev, curr = curr, prev
    63  	}
    64  	return prev[k]
    65  }