Numerical integration: Difference between revisions

→‎{{header|Go}}: fix float type per language change
m ('round' defaults to 3)
(→‎{{header|Go}}: fix float type per language change)
Line 1,007:
end module FunctionHolder</lang>
=={{header|Go}}==
<lang go>package main
The float type used here is 32 bits.
<lang go>
package main
 
import (
Line 1,019 ⟶ 1,017:
// specification for an integration
type spec struct {
lower, upper floatfloat64 // bounds for integration
n int // number of parts
exact floatfloat64 // expected answer
fs string // mathematical description of function
f func(floatfloat64) floatfloat64 // function to integrate
}
 
// test cases per task description
var data = []spec{
spec{0, 1, 100, .25, "x^3", func(x floatfloat64) floatfloat64 { return x * x * x }},
spec{1, 100, 1000, floatfloat64(math.Log(100)), "1/x", func(x floatfloat64) floatfloat64 { return 1 / x }},
spec{0, 5000, 5e5, 12.5e6, "x", func(x floatfloat64) floatfloat64 { return x }},
spec{0, 6000, 6e6, 18e6, "x", func(x floatfloat64) floatfloat64 { return x }},
}
 
Line 1,037 ⟶ 1,035:
type method struct {
name string
integrate func(spec) floatfloat64
}
 
Line 1,049 ⟶ 1,047:
}
 
func rectLeft(t spec) floatfloat64 {
parts := make([]floatfloat64, t.n)
r := t.upper - t.lower
nf := floatfloat64(t.n)
x0 := t.lower
for i := range parts {
x1 := t.lower + floatfloat64(i+1)*r/nf
// x1-x0 better than r/nf.
// (with r/nf, the represenation error accumulates)
Line 1,064 ⟶ 1,062:
}
 
func rectRight(t spec) floatfloat64 {
parts := make([]floatfloat64, t.n)
r := t.upper - t.lower
nf := floatfloat64(t.n)
x0 := t.lower
for i := range parts {
x1 := t.lower + floatfloat64(i+1)*r/nf
parts[i] = t.f(x1) * (x1 - x0)
x0 = x1
Line 1,077 ⟶ 1,075:
}
 
func rectMid(t spec) floatfloat64 {
parts := make([]floatfloat64, t.n)
r := t.upper - t.lower
nf := floatfloat64(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
Line 1,089 ⟶ 1,087:
x0 := t.lower - .5*r/nf
for i := range parts {
x1 := t.lower + (floatfloat64(i)+.5)*r/nf
parts[i] = t.f(x1) * (x1 - x0)
x0 = x1
Line 1,096 ⟶ 1,094:
}
 
func trap(t spec) floatfloat64 {
parts := make([]floatfloat64, t.n)
r := t.upper - t.lower
nf := floatfloat64(t.n)
x0 := t.lower
f0 := t.f(x0)
for i := range parts {
x1 := t.lower + floatfloat64(i+1)*r/nf
f1 := t.f(x1)
parts[i] = (f0 + f1) * .5 * (x1 - x0)
Line 1,111 ⟶ 1,109:
}
 
func simpson(t spec) floatfloat64 {
parts := make([]floatfloat64, 2*t.n+1)
r := t.upper - t.lower
nf := floatfloat64(t.n)
// similar to the rectangle midpoint logic explained above,
// we play a little loose with the values used for dx and dx0.
Line 1,122 ⟶ 1,120:
x0 := t.lower + dx0
for i := 1; i < t.n; i++ {
x1 := t.lower + floatfloat64(i+1)*r/nf
xmid := (x0 + x1) * .5
dx := x1 - x0
Line 1,131 ⟶ 1,129:
parts[2*t.n] = t.f(t.upper) * dx0
return sum(parts) / 6
}
 
// sum a list of numbers avoiding loss of precision
// (there might be better ways...)
func sum(parts []floatfloat64) floatfloat64 {
sort.SortFloats(parts)
ip := sort.SearchFloatsSortFloat64s(parts, 0)
inip := ip - 1sort.SearchFloat64s(parts, 0)
in := ip - 1
sum := parts[ip]
ip++
Line 1,187 ⟶ 1,186:
fmt.Println("")
}
}</lang go>
}
</lang>
Output:
<pre>Test case: f(x) = x^3
<pre>
Test case: f(x) = x^3
Integration from 0 to 1 in 100 parts
Exact result 2.5000000e-01 Error
Rectangular (left) 2.4502498e4502500e-01 4.9750209e9750000e-03
Rectangular (right) 2.5502500e-01 5.0249994e0250000e-03
Rectangular (midpoint) 2.4998754e4998750e-01 1.2457371e2500000e-05
Trapezium 2.5002500e-01 2.5004148e5000000e-05
Simpson's 2.4999999e5000000e-01 15.4901161e5511151e-0817
 
Test case: f(x) = 1/x
Integration from 1 to 100 in 1000 parts
Exact result 4.6051702e+00 Error
skeysRectangular byron 88> t(left) 4.6549916e6549911e+00 4.9821377e9820872e-02
Rectangular (right) 4.5569830e5569811e+00 4.8187256e8189128e-02
Rectangular (midpoint) 4.6047630e6047625e+00 4.0721893e0763731e-04
Trapezium 4.6059861e+00 8.1586838e1587153e-04
Simpson's 4.6051679e6051704e+00 21.3841858e9896905e-0607
 
Test case: f(x) = x
Integration from 0 to 5000 in 500000 parts
Exact result 1.2500000e+07 Error
Rectangular (left) 1.2500435e2499975e+07 42.3500000e5000000e+0201
Rectangular (right) 1.2500486e2500025e+07 42.8600000e5000000e+0201
Rectangular (midpoint) 1.2500465e2500000e+07 41.6500000e+026763806e-08
Trapezium 1.2500461e2500000e+07 43.6100000e+027252903e-08
Simpson's 1.2498455e2500000e+07 12.5450000e+039802322e-08
 
Test case: f(x) = x
Integration from 0 to 6000 in 6000000 parts
Exact result 1.8000000e+07 Error
Rectangular (left) 1.7911916e7999997e+07 83.8084000e0000000e+0400
Rectangular (right) 1.7911926e8000003e+07 83.8074000e0000000e+0400
Rectangular (midpoint) 1.7911904e8000000e+07 81.8096000e+044901161e-08
Trapezium 1.7911922e8000000e+07 81.8078000e+048626451e-08
Simpson's 1.7913002e8000000e+07 87.6998000e+044505806e-09</pre>
</pre>
 
=={{header|Haskell}}==
1,707

edits