Numerical integration: Difference between revisions

Content added Content deleted
(→‎{{header|Go}}: improved summation algorithm)
Line 1,336: Line 1,336:
"fmt"
"fmt"
"math"
"math"
"sort"
)
)


Line 1,351: Line 1,350:
var data = []spec{
var data = []spec{
spec{0, 1, 100, .25, "x^3", func(x float64) float64 { return x * x * x }},
spec{0, 1, 100, .25, "x^3", func(x float64) float64 { return x * x * x }},
spec{1, 100, 1000, float64(math.Log(100)), "1/x", func(x float64) float64 { 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 float64) float64 { return x }},
spec{0, 5000, 5e5, 12.5e6, "x", func(x float64) float64 { return x }},
spec{0, 6000, 6e6, 18e6, "x", func(x float64) float64 { return x }},
spec{0, 6000, 6e6, 18e6, "x", func(x float64) float64 { return x }},
}
}

// object for associating a printable function name with an integration method
// object for associating a printable function name with an integration method
type method struct {
type method struct {
name string
name string
integrate func(spec) float64
integrate func(spec) float64
}
}

// integration methods implemented per task description
// integration methods implemented per task description
var methods = []method{
var methods = []method{
Line 1,370: Line 1,370:
method{"Simpson's ", simpson},
method{"Simpson's ", simpson},
}
}

func rectLeft(t spec) float64 {
func rectLeft(t spec) float64 {
parts := make([]float64, t.n)
parts := make([]float64, t.n)
Line 1,456: Line 1,456:


// sum a list of numbers avoiding loss of precision
// sum a list of numbers avoiding loss of precision
⚫
func sum(v []float64) float64 {
// (there might be better ways...)
if len(v) == 0 {
⚫
func sum(parts []float64) float64 {
return 0
sort.Float64s(parts)
}
ip := sort.SearchFloat64s(parts, 0)
in := ip - 1
var parts []float64
sum := parts[ip]
for _, x := range v {
ip++
var i int
for {
for _, p := range parts {
if sum < 0 {
sum := p + x
if ip < len(parts) {
var err float64
// add the next biggest positive part
if math.Abs(x) < math.Abs(p) {
sum += parts[ip]
err = x - (sum - p)
⚫
ip++
} else {
} else {
// no more positive parts.
err = p - (sum - x)
// add remaining negatives and exit loop
for ; in >= 0; in-- {
⚫
sum += parts[in]
⚫
}
break
}
}
} else {
if err != 0 {
if in >= 0 {
parts[i] = err
// add the next biggest negative part
i++
sum += parts[in]
in--
} else {
// no more negative parts.
// add remaining positives and exit loop
for ; ip < len(parts); ip++ {
sum += parts[ip]
}
break
}
}
⚫
x = sum
}
}
parts = append(parts[:i], x)
}
var sum float64
for _, x := range parts {
⚫
sum += x
}
}
return sum
return sum
}
}

func main() {
func main() {
for _, t := range data {
for _, t := range data {
fmt.Println("Test case: f(x) =", t.fs)
fmt.Println("Test case: f(x) =", t.fs)
fmt.Println("Integration from", t.lower, "to", t.upper, "in", t.n, "parts")
fmt.Println("Integration from", t.lower, "to", t.upper,
⚫
"in", t.n, "parts")
fmt.Printf("Exact result %.7e Error\n", t.exact)
fmt.Printf("Exact result %.7e Error\n", t.exact)
for _, m := range methods {
for _, m := range methods {
Line 1,512: Line 1,504:
}</lang>
}</lang>
Output:
Output:
<pre>
<pre>Test case: f(x) = x^3
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
Line 1,519: Line 1,512:
Rectangular (midpoint) 2.4998750e-01 1.2500000e-05
Rectangular (midpoint) 2.4998750e-01 1.2500000e-05
Trapezium 2.5002500e-01 2.5000000e-05
Trapezium 2.5002500e-01 2.5000000e-05
Simpson's 2.5000000e-01 5.5511151e-17
Simpson's 2.5000000e-01 0.0000000e+00


Test case: f(x) = 1/x
Test case: f(x) = 1/x
Line 1,535: Line 1,528:
Rectangular (left) 1.2499975e+07 2.5000000e+01
Rectangular (left) 1.2499975e+07 2.5000000e+01
Rectangular (right) 1.2500025e+07 2.5000000e+01
Rectangular (right) 1.2500025e+07 2.5000000e+01
Rectangular (midpoint) 1.2500000e+07 1.6763806e-08
Rectangular (midpoint) 1.2500000e+07 0.0000000e+00
Trapezium 1.2500000e+07 3.7252903e-08
Trapezium 1.2500000e+07 0.0000000e+00
Simpson's 1.2500000e+07 2.9802322e-08
Simpson's 1.2500000e+07 0.0000000e+00


Test case: f(x) = x
Test case: f(x) = x
Line 1,544: Line 1,537:
Rectangular (left) 1.7999997e+07 3.0000000e+00
Rectangular (left) 1.7999997e+07 3.0000000e+00
Rectangular (right) 1.8000003e+07 3.0000000e+00
Rectangular (right) 1.8000003e+07 3.0000000e+00
Rectangular (midpoint) 1.8000000e+07 1.4901161e-08
Rectangular (midpoint) 1.8000000e+07 0.0000000e+00
Trapezium 1.8000000e+07 1.8626451e-08
Trapezium 1.8000000e+07 0.0000000e+00
Simpson's 1.8000000e+07 7.4505806e-09</pre>
Simpson's 1.8000000e+07 0.0000000e+00
</pre>


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