1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
|
// 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
}
|
