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