gonum.org/v1/gonum@v0.14.0/integrate/simpsons.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 "sort"
     8  
     9  // Simpsons returns an approximate value of the integral
    10  //
    11  //	\int_a^b f(x)dx
    12  //
    13  // computed using the Simpsons's method. The function f is given as a slice of
    14  // samples evaluated at locations in x, that is,
    15  //
    16  //	f[i] = f(x[i]), x[0] = a, x[len(x)-1] = b
    17  //
    18  // The slice x must be sorted in strictly increasing order. x and f must be of
    19  // equal length and the length must be at least 3.
    20  //
    21  // See https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule_for_irregularly_spaced_data
    22  // for more information.
    23  func Simpsons(x, f []float64) float64 {
    24  	n := len(x)
    25  	switch {
    26  	case len(f) != n:
    27  		panic("integrate: slice length mismatch")
    28  	case n < 3:
    29  		panic("integrate: input data too small")
    30  	case !sort.Float64sAreSorted(x):
    31  		panic("integrate: must be sorted")
    32  	}
    33  
    34  	var integral float64
    35  	for i := 1; i < n-1; i += 2 {
    36  		h0 := x[i] - x[i-1]
    37  		h1 := x[i+1] - x[i]
    38  		if h0 == 0 || h1 == 0 {
    39  			panic("integrate: repeated abscissa")
    40  		}
    41  		h0p2 := h0 * h0
    42  		h0p3 := h0 * h0 * h0
    43  		h1p2 := h1 * h1
    44  		h1p3 := h1 * h1 * h1
    45  		hph := h0 + h1
    46  		a0 := (2*h0p3 - h1p3 + 3*h1*h0p2) / (6 * h0 * hph)
    47  		a1 := (h0p3 + h1p3 + 3*h0*h1*hph) / (6 * h0 * h1)
    48  		a2 := (-h0p3 + 2*h1p3 + 3*h0*h1p2) / (6 * h1 * hph)
    49  		integral += a0 * f[i-1]
    50  		integral += a1 * f[i]
    51  		integral += a2 * f[i+1]
    52  	}
    53  
    54  	if n%2 == 0 {
    55  		h0 := x[n-2] - x[n-3]
    56  		h1 := x[n-1] - x[n-2]
    57  		if h0 == 0 || h1 == 0 {
    58  			panic("integrate: repeated abscissa")
    59  		}
    60  		h1p2 := h1 * h1
    61  		h1p3 := h1 * h1 * h1
    62  		hph := h0 + h1
    63  		a0 := -1 * h1p3 / (6 * h0 * hph)
    64  		a1 := (h1p2 + 3*h0*h1) / (6 * h0)
    65  		a2 := (2*h1p2 + 3*h0*h1) / (6 * hph)
    66  		integral += a0 * f[n-3]
    67  		integral += a1 * f[n-2]
    68  		integral += a2 * f[n-1]
    69  	}
    70  
    71  	return integral
    72  }