202 lines
7.0 KiB
Go
202 lines
7.0 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
|
// SPDX-License-Identifier: MIT
|
|
|
|
package integrate
|
|
|
|
import (
|
|
"sourcedock.dev/petrbalvin/tensor/internal/base"
|
|
)
|
|
|
|
import (
|
|
"math"
|
|
"testing"
|
|
)
|
|
|
|
func quadValue(t *testing.T, f func(float64) (float64, error), a, b float64, opts QuadratureOptions) (float64, float64) {
|
|
t.Helper()
|
|
v, errEst, err := IntegrateFunction(f, a, b, opts)
|
|
if err != nil {
|
|
t.Fatalf("Integrate: %v", err)
|
|
}
|
|
return v, errEst
|
|
}
|
|
|
|
func TestIntegratePolynomials(t *testing.T) {
|
|
// Degree 14 on 21 points: exact by construction.
|
|
v, _ := quadValue(t, func(x float64) (float64, error) { return math.Pow(x, 7), nil }, 0, 1, QuadratureOptions{})
|
|
if math.Abs(v-1.0/8) > 1e-14 {
|
|
t.Fatalf("integral of x^7 = %.16g, want 0.125", v)
|
|
}
|
|
v, _ = quadValue(t, func(x float64) (float64, error) { return math.Pow(x, 14), nil }, -1, 1, QuadratureOptions{})
|
|
if math.Abs(v-2.0/15) > 1e-14 {
|
|
t.Fatalf("integral of x^14 = %.16g, want %g", v, 2.0/15)
|
|
}
|
|
}
|
|
|
|
func TestIntegrateSmooth(t *testing.T) {
|
|
cases := []struct {
|
|
name string
|
|
f func(float64) (float64, error)
|
|
a, b float64
|
|
want float64
|
|
}{
|
|
{"sine over a period", func(x float64) (float64, error) { return math.Sin(x), nil }, 0, math.Pi, 2},
|
|
{"exponential", func(x float64) (float64, error) { return math.Exp(x), nil }, -1, 1, 2 * math.Sinh(1)},
|
|
{"arctangent derivative", func(x float64) (float64, error) { return 1 / (1 + x*x), nil }, 0, 1, math.Pi / 4},
|
|
{"gaussian", func(x float64) (float64, error) { return math.Exp(-x * x), nil }, 0, 5, 0.5 * math.Sqrt(math.Pi)},
|
|
}
|
|
for _, c := range cases {
|
|
v, est := quadValue(t, c.f, c.a, c.b, QuadratureOptions{})
|
|
if math.Abs(v-c.want) > 1e-11 {
|
|
t.Fatalf("%s: %.14g, want %.14g", c.name, v, c.want)
|
|
}
|
|
if math.Abs(v-c.want) > 10*est+1e-14 {
|
|
t.Fatalf("%s: error estimate %g understates the true error %g", c.name, est, math.Abs(v-c.want))
|
|
}
|
|
}
|
|
}
|
|
|
|
func TestIntegrateOscillatory(t *testing.T) {
|
|
v, _ := quadValue(t, func(x float64) (float64, error) { return math.Sin(x), nil }, 0, 10*math.Pi, QuadratureOptions{})
|
|
if math.Abs(v) > 1e-9 {
|
|
t.Fatalf("ten sine periods integrate to %g, want 0", v)
|
|
}
|
|
v, _ = quadValue(t, func(x float64) (float64, error) { return math.Sin(30*x) * math.Exp(-x), nil },
|
|
0, 40, QuadratureOptions{RelTol: 1e-9})
|
|
want := 30.0 / 901.0 // exact over [0, inf): 30/(30^2 + 1)
|
|
if math.Abs(v-want) > 1e-8 {
|
|
t.Fatalf("damped oscillation = %.14g, want %.14g", v, want)
|
|
}
|
|
}
|
|
|
|
func TestIntegrateSharpPeak(t *testing.T) {
|
|
// A Lorentzian a hundred times narrower than the interval forces
|
|
// deep subdivision; the answer must land on the exact closed form
|
|
// 2·arctan(100). (A peak narrower than the first rule's node
|
|
// spacing would be invisible to any sampled scheme, which the doc
|
|
// contract states.)
|
|
v, _ := quadValue(t, func(x float64) (float64, error) {
|
|
return 10 / (1 + 100*x*x), nil
|
|
}, -10, 10, QuadratureOptions{RelTol: 1e-10})
|
|
if want := 2 * math.Atan(100); math.Abs(v-want) > 1e-9 {
|
|
t.Fatalf("narrow Lorentzian = %.12g, want %.12g", v, want)
|
|
}
|
|
}
|
|
|
|
func TestIntegrateInfiniteBounds(t *testing.T) {
|
|
cases := []struct {
|
|
name string
|
|
f func(float64) (float64, error)
|
|
a, b float64
|
|
want float64
|
|
}{
|
|
{"exponential tail", func(x float64) (float64, error) { return math.Exp(-x), nil }, 0, math.Inf(1), 1},
|
|
{"cauchy tail", func(x float64) (float64, error) { return 1 / (1 + x*x), nil }, 0, math.Inf(1), math.Pi / 2},
|
|
{"full gaussian", func(x float64) (float64, error) { return math.Exp(-x * x), nil }, math.Inf(-1), math.Inf(1), math.Sqrt(math.Pi)},
|
|
{"negative tail", func(x float64) (float64, error) { return math.Exp(x), nil }, math.Inf(-1), 0, 1},
|
|
}
|
|
for _, c := range cases {
|
|
v, _ := quadValue(t, c.f, c.a, c.b, QuadratureOptions{})
|
|
if math.Abs(v-c.want) > 1e-9 {
|
|
t.Fatalf("%s: %.14g, want %.14g", c.name, v, c.want)
|
|
}
|
|
}
|
|
}
|
|
|
|
func TestIntegrateOrientation(t *testing.T) {
|
|
sin := func(x float64) (float64, error) { return math.Sin(x), nil }
|
|
forward, _ := quadValue(t, sin, 0, math.Pi, QuadratureOptions{})
|
|
backward, _ := quadValue(t, sin, math.Pi, 0, QuadratureOptions{})
|
|
if math.Abs(backward+forward) > 1e-14 {
|
|
t.Fatalf("reversed integral = %g, want %g", backward, -forward)
|
|
}
|
|
if v, _ := quadValue(t, sin, 2, 2, QuadratureOptions{}); v != 0 {
|
|
t.Fatalf("empty interval = %g, want 0", v)
|
|
}
|
|
}
|
|
|
|
func TestIntegrateErrors(t *testing.T) {
|
|
boom := func(x float64) (float64, error) {
|
|
if x > 0.5 {
|
|
return 0, base.Errf("integrand failed at %g", x)
|
|
}
|
|
return 1, nil
|
|
}
|
|
if _, _, err := IntegrateFunction(boom, 0, 1, QuadratureOptions{}); err == nil {
|
|
t.Fatal("integrand error: want an error")
|
|
}
|
|
if _, _, err := IntegrateFunction(func(x float64) (float64, error) { return x, nil },
|
|
math.NaN(), 1, QuadratureOptions{}); err == nil {
|
|
t.Fatal("NaN bound: want an error")
|
|
}
|
|
// A flat integrand over an infinite interval diverges: the
|
|
// adaptation must report it, not return a number.
|
|
if _, _, err := IntegrateFunction(func(x float64) (float64, error) { return 1, nil },
|
|
0, math.Inf(1), QuadratureOptions{MaxIntervals: 8}); err == nil {
|
|
t.Fatal("divergent integral: want an error")
|
|
}
|
|
// A tolerance no budget can meet must be reported.
|
|
if _, _, err := IntegrateFunction(func(x float64) (float64, error) { return math.Exp(-1000 * (x - 0.5) * (x - 0.5)), nil },
|
|
0, 1, QuadratureOptions{RelTol: 1e-16, AbsTol: 0, MaxIntervals: 4}); err == nil {
|
|
t.Fatal("exhausted budget: want an error")
|
|
}
|
|
}
|
|
|
|
func TestGaussLegendreNodes(t *testing.T) {
|
|
nodes, weights, err := GaussLegendreNodes(8)
|
|
if err != nil {
|
|
t.Fatalf("GaussLegendreNodes: %v", err)
|
|
}
|
|
total := 0.0
|
|
for i := range 8 {
|
|
total += weights[i]
|
|
if math.Abs(nodes[i]+nodes[7-i]) > 1e-14 {
|
|
t.Fatalf("nodes %d and %d are not symmetric", i, 7-i)
|
|
}
|
|
if i > 0 && nodes[i] <= nodes[i-1] {
|
|
t.Fatalf("nodes are not ascending at %d", i)
|
|
}
|
|
}
|
|
if math.Abs(total-2) > 1e-14 {
|
|
t.Fatalf("weight sum = %.16g, want 2", total)
|
|
}
|
|
// A 16th degree polynomial integrates exactly on 8 points.
|
|
moment := 0.0
|
|
for i := range 8 {
|
|
moment += weights[i] * math.Pow(nodes[i], 14)
|
|
}
|
|
if math.Abs(moment-2.0/15) > 1e-14 {
|
|
t.Fatalf("moment of x^14 = %.16g, want %g", moment, 2.0/15)
|
|
}
|
|
// The cache returns the same arrays.
|
|
again, _, err := GaussLegendreNodes(8)
|
|
if err != nil {
|
|
t.Fatalf("cached call: %v", err)
|
|
}
|
|
if &again[0] != &nodes[0] {
|
|
t.Fatal("cached nodes are not the cached arrays")
|
|
}
|
|
if _, _, err := GaussLegendreNodes(0); err == nil {
|
|
t.Fatal("n = 0: want an error")
|
|
}
|
|
if _, _, err := GaussLegendreNodes(129); err == nil {
|
|
t.Fatal("n = 129: want an error")
|
|
}
|
|
}
|
|
|
|
// TestIntegrateFunctionNaNIsError pins the NaN contract: an integrand
|
|
// that yields NaN makes IntegrateFunction return an error instead of
|
|
// a silent NaN integral (the first estimate and every bisected leaf
|
|
// are checked).
|
|
func TestIntegrateFunctionNaNIsError(t *testing.T) {
|
|
f := func(x float64) (float64, error) {
|
|
if math.Abs(x) > 0.25 {
|
|
return math.NaN(), nil
|
|
}
|
|
return 1, nil
|
|
}
|
|
if _, _, err := IntegrateFunction(f, -1, 1, QuadratureOptions{}); err == nil {
|
|
t.Fatal("expected an error for an integrand that returns NaN")
|
|
}
|
|
}
|