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