Numerical integration: Difference between revisions

Content added Content deleted
m ('round' defaults to 3)
(→‎{{header|Go}}: fix float type per language change)
Line 1,007: Line 1,007:
end module FunctionHolder</lang>
end module FunctionHolder</lang>
=={{header|Go}}==
=={{header|Go}}==
⚫
<lang go>package main
The float type used here is 32 bits.
⚫
<lang go>
⚫
package main


import (
import (
Line 1,019: Line 1,017:
// specification for an integration
// specification for an integration
type spec struct {
type spec struct {
lower, upper float // bounds for integration
lower, upper float64 // bounds for integration
n int // number of parts
n int // number of parts
exact float // expected answer
exact float64 // expected answer
fs string // mathematical description of function
fs string // mathematical description of function
f func(float) float // function to integrate
f func(float64) float64 // function to integrate
}
}


// test cases per task description
// test cases per task description
var data = []spec{
var data = []spec{
spec{0, 1, 100, .25, "x^3", func(x float) float { return x * x * x }},
spec{0, 1, 100, .25, "x^3", func(x float64) float64 { return x * x * x }},
spec{1, 100, 1000, float(math.Log(100)), "1/x", func(x float) float { return 1 / x }},
spec{1, 100, 1000, float64(math.Log(100)), "1/x", func(x float64) float64 { return 1 / x }},
spec{0, 5000, 5e5, 12.5e6, "x", func(x float) float { return x }},
spec{0, 5000, 5e5, 12.5e6, "x", func(x float64) float64 { return x }},
spec{0, 6000, 6e6, 18e6, "x", func(x float) float { return x }},
spec{0, 6000, 6e6, 18e6, "x", func(x float64) float64 { return x }},
}
}


Line 1,037: Line 1,035:
type method struct {
type method struct {
name string
name string
integrate func(spec) float
integrate func(spec) float64
}
}


Line 1,049: Line 1,047:
}
}


func rectLeft(t spec) float {
func rectLeft(t spec) float64 {
parts := make([]float, t.n)
parts := make([]float64, t.n)
r := t.upper - t.lower
r := t.upper - t.lower
nf := float(t.n)
nf := float64(t.n)
x0 := t.lower
x0 := t.lower
for i := range parts {
for i := range parts {
x1 := t.lower + float(i+1)*r/nf
x1 := t.lower + float64(i+1)*r/nf
// x1-x0 better than r/nf.
// x1-x0 better than r/nf.
// (with r/nf, the represenation error accumulates)
// (with r/nf, the represenation error accumulates)
Line 1,064: Line 1,062:
}
}


func rectRight(t spec) float {
func rectRight(t spec) float64 {
parts := make([]float, t.n)
parts := make([]float64, t.n)
r := t.upper - t.lower
r := t.upper - t.lower
nf := float(t.n)
nf := float64(t.n)
x0 := t.lower
x0 := t.lower
for i := range parts {
for i := range parts {
x1 := t.lower + float(i+1)*r/nf
x1 := t.lower + float64(i+1)*r/nf
parts[i] = t.f(x1) * (x1 - x0)
parts[i] = t.f(x1) * (x1 - x0)
x0 = x1
x0 = x1
Line 1,077: Line 1,075:
}
}


func rectMid(t spec) float {
func rectMid(t spec) float64 {
parts := make([]float, t.n)
parts := make([]float64, t.n)
r := t.upper - t.lower
r := t.upper - t.lower
nf := float(t.n)
nf := float64(t.n)
// there's a tiny gloss in the x1-x0 trick here. the correct way
// 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
// would be to compute x's at division boundaries, but we don't need
Line 1,089: Line 1,087:
x0 := t.lower - .5*r/nf
x0 := t.lower - .5*r/nf
for i := range parts {
for i := range parts {
x1 := t.lower + (float(i)+.5)*r/nf
x1 := t.lower + (float64(i)+.5)*r/nf
parts[i] = t.f(x1) * (x1 - x0)
parts[i] = t.f(x1) * (x1 - x0)
x0 = x1
x0 = x1
Line 1,096: Line 1,094:
}
}


func trap(t spec) float {
func trap(t spec) float64 {
parts := make([]float, t.n)
parts := make([]float64, t.n)
r := t.upper - t.lower
r := t.upper - t.lower
nf := float(t.n)
nf := float64(t.n)
x0 := t.lower
x0 := t.lower
f0 := t.f(x0)
f0 := t.f(x0)
for i := range parts {
for i := range parts {
x1 := t.lower + float(i+1)*r/nf
x1 := t.lower + float64(i+1)*r/nf
f1 := t.f(x1)
f1 := t.f(x1)
parts[i] = (f0 + f1) * .5 * (x1 - x0)
parts[i] = (f0 + f1) * .5 * (x1 - x0)
Line 1,111: Line 1,109:
}
}


func simpson(t spec) float {
func simpson(t spec) float64 {
parts := make([]float, 2*t.n+1)
parts := make([]float64, 2*t.n+1)
r := t.upper - t.lower
r := t.upper - t.lower
nf := float(t.n)
nf := float64(t.n)
// similar to the rectangle midpoint logic explained above,
// similar to the rectangle midpoint logic explained above,
// we play a little loose with the values used for dx and dx0.
// we play a little loose with the values used for dx and dx0.
Line 1,122: Line 1,120:
x0 := t.lower + dx0
x0 := t.lower + dx0
for i := 1; i < t.n; i++ {
for i := 1; i < t.n; i++ {
x1 := t.lower + float(i+1)*r/nf
x1 := t.lower + float64(i+1)*r/nf
xmid := (x0 + x1) * .5
xmid := (x0 + x1) * .5
dx := x1 - x0
dx := x1 - x0
Line 1,131: Line 1,129:
parts[2*t.n] = t.f(t.upper) * dx0
parts[2*t.n] = t.f(t.upper) * dx0
return sum(parts) / 6
return sum(parts) / 6
}
}

// sum a list of numbers avoiding loss of precision
// sum a list of numbers avoiding loss of precision
// (there might be better ways...)
func sum(parts []float) float {
func sum(parts []float64) float64 {
sort.SortFloats(parts)
ip := sort.SearchFloats(parts, 0)
sort.SortFloat64s(parts)
in := ip - 1
ip := sort.SearchFloat64s(parts, 0)
in := ip - 1
sum := parts[ip]
sum := parts[ip]
ip++
ip++
Line 1,187: Line 1,186:
fmt.Println("")
fmt.Println("")
}
}
⚫
}</lang>
}
</lang>
Output:
Output:
⚫
<pre>Test case: f(x) = x^3
<pre>
⚫
Test case: f(x) = x^3
Integration from 0 to 1 in 100 parts
Integration from 0 to 1 in 100 parts
Exact result 2.5000000e-01 Error
Exact result 2.5000000e-01 Error
Rectangular (left) 2.4502498e-01 4.9750209e-03
Rectangular (left) 2.4502500e-01 4.9750000e-03
Rectangular (right) 2.5502500e-01 5.0249994e-03
Rectangular (right) 2.5502500e-01 5.0250000e-03
Rectangular (midpoint) 2.4998754e-01 1.2457371e-05
Rectangular (midpoint) 2.4998750e-01 1.2500000e-05
Trapezium 2.5002500e-01 2.5004148e-05
Trapezium 2.5002500e-01 2.5000000e-05
Simpson's 2.4999999e-01 1.4901161e-08
Simpson's 2.5000000e-01 5.5511151e-17


Test case: f(x) = 1/x
Test case: f(x) = 1/x
Integration from 1 to 100 in 1000 parts
Integration from 1 to 100 in 1000 parts
Exact result 4.6051702e+00 Error
Exact result 4.6051702e+00 Error
skeys byron 88> t) 4.6549916e+00 4.9821377e-02
Rectangular (left) 4.6549911e+00 4.9820872e-02
Rectangular (right) 4.5569830e+00 4.8187256e-02
Rectangular (right) 4.5569811e+00 4.8189128e-02
Rectangular (midpoint) 4.6047630e+00 4.0721893e-04
Rectangular (midpoint) 4.6047625e+00 4.0763731e-04
Trapezium 4.6059861e+00 8.1586838e-04
Trapezium 4.6059861e+00 8.1587153e-04
Simpson's 4.6051679e+00 2.3841858e-06
Simpson's 4.6051704e+00 1.9896905e-07


Test case: f(x) = x
Test case: f(x) = x
Integration from 0 to 5000 in 500000 parts
Integration from 0 to 5000 in 500000 parts
Exact result 1.2500000e+07 Error
Exact result 1.2500000e+07 Error
Rectangular (left) 1.2500435e+07 4.3500000e+02
Rectangular (left) 1.2499975e+07 2.5000000e+01
Rectangular (right) 1.2500486e+07 4.8600000e+02
Rectangular (right) 1.2500025e+07 2.5000000e+01
Rectangular (midpoint) 1.2500465e+07 4.6500000e+02
Rectangular (midpoint) 1.2500000e+07 1.6763806e-08
Trapezium 1.2500461e+07 4.6100000e+02
Trapezium 1.2500000e+07 3.7252903e-08
Simpson's 1.2498455e+07 1.5450000e+03
Simpson's 1.2500000e+07 2.9802322e-08


Test case: f(x) = x
Test case: f(x) = x
Integration from 0 to 6000 in 6000000 parts
Integration from 0 to 6000 in 6000000 parts
Exact result 1.8000000e+07 Error
Exact result 1.8000000e+07 Error
Rectangular (left) 1.7911916e+07 8.8084000e+04
Rectangular (left) 1.7999997e+07 3.0000000e+00
Rectangular (right) 1.7911926e+07 8.8074000e+04
Rectangular (right) 1.8000003e+07 3.0000000e+00
Rectangular (midpoint) 1.7911904e+07 8.8096000e+04
Rectangular (midpoint) 1.8000000e+07 1.4901161e-08
Trapezium 1.7911922e+07 8.8078000e+04
Trapezium 1.8000000e+07 1.8626451e-08
Simpson's 1.7913002e+07 8.6998000e+04
Simpson's 1.8000000e+07 7.4505806e-09</pre>
</pre>


=={{header|Haskell}}==
=={{header|Haskell}}==
Retrieved from "https://rosettacode.org/wiki/Numerical_integration"