Numerical integration: Difference between revisions
Go solution
(Go solution) |
|||
Line 1,006:
end module FunctionHolder</lang>
=={{header|Go}}==
The float type used here is 32 bits.
<lang go>
package main
import (
"fmt"
"math"
"sort"
)
// specification for an integration
type spec struct {
lower, upper float // bounds for integration
n int // number of parts
exact float // expected answer
fs string // mathematical description of function
f func(float) float // function to integrate
}
// test cases per task description
var data = []spec{
spec{0, 1, 100, .25, "x^3", func(x float) float { return x * x * x }},
spec{1, 100, 1000, float(math.Log(100)), "1/x", func(x float) float { return 1 / x }},
spec{0, 5000, 5e5, 12.5e6, "x", func(x float) float { return x }},
spec{0, 6000, 6e6, 18e6, "x", func(x float) float { return x }},
}
// object for associating a printable function name with an integration method
type method struct {
name string
integrate func(spec) float
}
// integration methods implemented per task description
var methods = []method{
method{"Rectangular (left) ", rectLeft},
method{"Rectangular (right) ", rectRight},
method{"Rectangular (midpoint)", rectMid},
method{"Trapezium ", trap},
method{"Simpson's ", simpson},
}
func rectLeft(t spec) float {
parts := make([]float, t.n)
r := t.upper - t.lower
nf := float(t.n)
x0 := t.lower
for i := range parts {
x1 := t.lower + float(i+1)*r/nf
// x1-x0 better than r/nf.
// (with r/nf, the represenation error accumulates)
parts[i] = t.f(x0) * (x1 - x0)
x0 = x1
}
return sum(parts)
}
func rectRight(t spec) float {
parts := make([]float, t.n)
r := t.upper - t.lower
nf := float(t.n)
x0 := t.lower
for i := range parts {
x1 := t.lower + float(i+1)*r/nf
parts[i] = t.f(x1) * (x1 - x0)
x0 = x1
}
return sum(parts)
}
func rectMid(t spec) float {
parts := make([]float, t.n)
r := t.upper - t.lower
nf := float(t.n)
// there's a tiny gloss in the x1-x0 trick here. the correct way
// would be to compute x's at division boundaries, but we don't need
// those x's for anything else. (the function is evaluated on x's
// at division midpoints rather than division boundaries.) so, we
// reuse the midpoint x's, knowing that they will average out just
// as well. we just need one extra point, so we use lower-.5.
x0 := t.lower - .5*r/nf
for i := range parts {
x1 := t.lower + (float(i)+.5)*r/nf
parts[i] = t.f(x1) * (x1 - x0)
x0 = x1
}
return sum(parts)
}
func trap(t spec) float {
parts := make([]float, t.n)
r := t.upper - t.lower
nf := float(t.n)
x0 := t.lower
f0 := t.f(x0)
for i := range parts {
x1 := t.lower + float(i+1)*r/nf
f1 := t.f(x1)
parts[i] = (f0 + f1) * .5 * (x1 - x0)
x0, f0 = x1, f1
}
return sum(parts)
}
func simpson(t spec) float {
parts := make([]float, 2*t.n+1)
r := t.upper - t.lower
nf := float(t.n)
// similar to the rectangle midpoint logic explained above,
// we play a little loose with the values used for dx and dx0.
dx0 := r / nf
parts[0] = t.f(t.lower) * dx0
parts[1] = t.f(t.lower+dx0*.5) * dx0 * 4
x0 := t.lower + dx0
for i := 1; i < t.n; i++ {
x1 := t.lower + float(i+1)*r/nf
xmid := (x0 + x1) * .5
dx := x1 - x0
parts[2*i] = t.f(x0) * dx * 2
parts[2*i+1] = t.f(xmid) * dx * 4
x0 = x1
}
parts[2*t.n] = t.f(t.upper) * dx0
return sum(parts) / 6
}
// sum a list of numbers avoiding loss of precision
func sum(parts []float) float {
sort.SortFloats(parts)
ip := sort.SearchFloats(parts, 0)
in := ip - 1
sum := parts[ip]
ip++
for {
if sum < 0 {
if ip < len(parts) {
// add the next biggest positive part
sum += parts[ip]
ip++
} else {
// no more positive parts.
// add remaining negatives and exit loop
for ; in >= 0; in-- {
sum += parts[in]
}
break
}
} else {
if in >= 0 {
// add the next biggest negative part
sum += parts[in]
in--
} else {
// no more negative parts.
// add remaining positives and exit loop
for ; ip < len(parts); ip++ {
sum += parts[ip]
}
break
}
}
}
return sum
}
func main() {
for _, t := range data {
fmt.Println("Test case: f(x) =", t.fs)
fmt.Println("Integration from", t.lower, "to", t.upper, "in", t.n, "parts")
fmt.Printf("Exact result %.7e Error\n", t.exact)
for _, m := range methods {
a := m.integrate(t)
e := a - t.exact
if e < 0 {
e = -e
}
fmt.Printf("%s %.7e %.7e\n", m.name, a, e)
}
fmt.Println("")
}
}
</lang>
Output:
<pre>
Test case: f(x) = x^3
Integration from 0 to 1 in 100 parts
Exact result 2.5000000e-01 Error
Rectangular (left) 2.4502498e-01 4.9750209e-03
Rectangular (right) 2.5502500e-01 5.0249994e-03
Rectangular (midpoint) 2.4998754e-01 1.2457371e-05
Trapezium 2.5002500e-01 2.5004148e-05
Simpson's 2.4999999e-01 1.4901161e-08
Test case: f(x) = 1/x
Integration from 1 to 100 in 1000 parts
Exact result 4.6051702e+00 Error
skeys byron 88> t) 4.6549916e+00 4.9821377e-02
Rectangular (right) 4.5569830e+00 4.8187256e-02
Rectangular (midpoint) 4.6047630e+00 4.0721893e-04
Trapezium 4.6059861e+00 8.1586838e-04
Simpson's 4.6051679e+00 2.3841858e-06
Test case: f(x) = x
Integration from 0 to 5000 in 500000 parts
Exact result 1.2500000e+07 Error
Rectangular (left) 1.2500435e+07 4.3500000e+02
Rectangular (right) 1.2500486e+07 4.8600000e+02
Rectangular (midpoint) 1.2500465e+07 4.6500000e+02
Trapezium 1.2500461e+07 4.6100000e+02
Simpson's 1.2498455e+07 1.5450000e+03
Test case: f(x) = x
Integration from 0 to 6000 in 6000000 parts
Exact result 1.8000000e+07 Error
Rectangular (left) 1.7911916e+07 8.8084000e+04
Rectangular (right) 1.7911926e+07 8.8074000e+04
Rectangular (midpoint) 1.7911904e+07 8.8096000e+04
Trapezium 1.7911922e+07 8.8078000e+04
Simpson's 1.7913002e+07 8.6998000e+04
</pre>
=={{header|Haskell}}==
| |||