blob: ab221c7e47bd5a5cefaa018cb9779cfa4c511e97 [file] [log] [blame]
 // 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 "sort" // Simpsons returns an approximate value of the integral // \int_a^b f(x)dx // computed using the Simpsons's method. The function f is given as a slice of // samples evaluated at locations in x, that is, // f[i] = f(x[i]), x[0] = a, x[len(x)-1] = b // The slice x must be sorted in strictly increasing order. x and f must be of // equal length and the length must be at least 3. // // See https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule_for_irregularly_spaced_data // for more information. func Simpsons(x, f []float64) float64 { n := len(x) switch { case len(f) != n: panic("integrate: slice length mismatch") case n < 3: panic("integrate: input data too small") case !sort.Float64sAreSorted(x): panic("integrate: must be sorted") } var integral float64 for i := 1; i < n-1; i += 2 { h0 := x[i] - x[i-1] h1 := x[i+1] - x[i] if h0 == 0 || h1 == 0 { panic("integrate: repeated abscissa") } h0p2 := h0 * h0 h0p3 := h0 * h0 * h0 h1p2 := h1 * h1 h1p3 := h1 * h1 * h1 hph := h0 + h1 a0 := (2*h0p3 - h1p3 + 3*h1*h0p2) / (6 * h0 * hph) a1 := (h0p3 + h1p3 + 3*h0*h1*hph) / (6 * h0 * h1) a2 := (-h0p3 + 2*h1p3 + 3*h0*h1p2) / (6 * h1 * hph) integral += a0 * f[i-1] integral += a1 * f[i] integral += a2 * f[i+1] } if n%2 == 0 { h0 := x[n-2] - x[n-3] h1 := x[n-1] - x[n-2] if h0 == 0 || h1 == 0 { panic("integrate: repeated abscissa") } h1p2 := h1 * h1 h1p3 := h1 * h1 * h1 hph := h0 + h1 a0 := -1 * h1p3 / (6 * h0 * hph) a1 := (h1p2 + 3*h0*h1) / (6 * h0) a2 := (2*h1p2 + 3*h0*h1) / (6 * hph) integral += a0 * f[n-3] integral += a1 * f[n-2] integral += a2 * f[n-1] } return integral }