226 lines
8.4 KiB
Go
226 lines
8.4 KiB
Go
|
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
|||
|
|
// SPDX-License-Identifier: MIT
|
|||
|
|
|
|||
|
|
package integrate_test
|
|||
|
|
|
|||
|
|
// Runnable examples for the package: the flagship workflows, each
|
|||
|
|
// with a fixed output that `go test` checks, so the printed
|
|||
|
|
// documentation cannot drift from the code.
|
|||
|
|
|
|||
|
|
import (
|
|||
|
|
"fmt"
|
|||
|
|
"log"
|
|||
|
|
"math"
|
|||
|
|
|
|||
|
|
tensor "sourcedock.dev/petrbalvin/tensor"
|
|||
|
|
"sourcedock.dev/petrbalvin/tensor/integrate"
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
// The stiff scalar problem y' = −1000·(y − cos t) − sin t with
|
|||
|
|
// y(0) = 1, whose exact solution is y = cos t. The variable-order,
|
|||
|
|
// variable-step BDF scheme takes the long steps the solution's
|
|||
|
|
// smoothness allows where a fixed small step would be forced by the
|
|||
|
|
// fast transient, and BDFVarStats reports what it did.
|
|||
|
|
func ExampleIntegrateBDFVar() {
|
|||
|
|
f := func(t float64, y *tensor.Array) (*tensor.Array, error) {
|
|||
|
|
return tensor.FromFloats([]float64{-1000*(y.FloatAt(0)-math.Cos(t)) - math.Sin(t)}, 1)
|
|||
|
|
}
|
|||
|
|
y0, _ := tensor.FromFloats([]float64{1}, 1)
|
|||
|
|
var stats integrate.BDFVarStats
|
|||
|
|
y, err := integrate.IntegrateBDFVar(f, 0, 1, y0, integrate.BDFVarOptions{Stats: &stats})
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Printf("y(1) = %.6f, exact %.6f\n", y.FloatAt(0), math.Cos(1))
|
|||
|
|
fmt.Printf("accepted %d steps, rejected %d, highest order %d\n", stats.Steps, stats.Rejected, stats.MaxOrder)
|
|||
|
|
// Output:
|
|||
|
|
// y(1) = 0.540302, exact 0.540302
|
|||
|
|
// accepted 21 steps, rejected 1, highest order 5
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Event detection along a trajectory: the oscillator y″ = −y started
|
|||
|
|
// at y = (1, 0) passes the level y = 0.5 falling at t = π/3 and rising
|
|||
|
|
// at t = 5π/3. Each watch carries its own direction filter, and
|
|||
|
|
// IntegrateODEEvents returns the crossings sorted by time alongside
|
|||
|
|
// the final state.
|
|||
|
|
func ExampleIntegrateODEEvents() {
|
|||
|
|
f := func(t float64, y *tensor.Array) (*tensor.Array, error) {
|
|||
|
|
return tensor.FromFloats([]float64{y.FloatAt(1), -y.FloatAt(0)}, 2)
|
|||
|
|
}
|
|||
|
|
y0, _ := tensor.FromFloats([]float64{1, 0}, 2)
|
|||
|
|
// The same level function twice, with opposite direction filters;
|
|||
|
|
// Direction 0 would record both crossings on one watch.
|
|||
|
|
level := func(t float64, y *tensor.Array) (float64, error) { return y.FloatAt(0) - 0.5, nil }
|
|||
|
|
watches := []integrate.ODEWatch{
|
|||
|
|
{Function: level, Direction: -1},
|
|||
|
|
{Function: level, Direction: +1},
|
|||
|
|
}
|
|||
|
|
hits, final, err := integrate.IntegrateODEEvents(f, 0, 7, y0, watches, integrate.ODEOptions{})
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
for _, h := range hits {
|
|||
|
|
direction := "falling"
|
|||
|
|
if h.Rising {
|
|||
|
|
direction = "rising"
|
|||
|
|
}
|
|||
|
|
fmt.Printf("watch %d fired at t = %.4f (%s), y = %.4f\n", h.Watch, h.Time, direction, h.State.FloatAt(0))
|
|||
|
|
}
|
|||
|
|
fmt.Printf("y(7) = %.4f\n", final.FloatAt(0))
|
|||
|
|
// Output:
|
|||
|
|
// watch 0 fired at t = 1.0472 (falling), y = 0.5000
|
|||
|
|
// watch 1 fired at t = 5.2360 (rising), y = 0.5000
|
|||
|
|
// y(7) = 0.7539
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// A symplectic integrator on a separable Hamiltonian: the harmonic
|
|||
|
|
// oscillator q″ = −q with unit mass, q(0) = 1 and p(0) = 0, whose
|
|||
|
|
// energy ½(p² + q²) stays in a bounded band instead of drifting.
|
|||
|
|
// The step stays fixed by design; only the number of steps is chosen.
|
|||
|
|
func ExampleIntegrateVerlet() {
|
|||
|
|
accel := func(q *tensor.Array) (*tensor.Array, error) {
|
|||
|
|
return tensor.FromFloats([]float64{-q.FloatAt(0)}, 1)
|
|||
|
|
}
|
|||
|
|
q0, _ := tensor.FromFloats([]float64{1}, 1)
|
|||
|
|
p0, _ := tensor.FromFloats([]float64{0}, 1)
|
|||
|
|
const steps = 4000
|
|||
|
|
positions, momenta, err := integrate.IntegrateVerlet(accel, 0, 2*math.Pi, q0, p0, steps)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
worst := 0.0
|
|||
|
|
for s := range steps + 1 {
|
|||
|
|
q, p := positions[s].FloatAt(0), momenta[s].FloatAt(0)
|
|||
|
|
drift := math.Abs(0.5*(p*p+q*q) - 0.5)
|
|||
|
|
worst = math.Max(worst, drift)
|
|||
|
|
}
|
|||
|
|
fmt.Printf("q(2π) = %.6f, p(2π) = %.2e\n", positions[steps].FloatAt(0), momenta[steps].FloatAt(0))
|
|||
|
|
fmt.Printf("worst energy deviation over the period: %.2e\n", worst)
|
|||
|
|
// Output:
|
|||
|
|
// q(2π) = 1.000000, p(2π) = -6.46e-07
|
|||
|
|
// worst energy deviation over the period: 3.08e-07
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Quadrature and cubature: a Gauss-Legendre rule read from the
|
|||
|
|
// package's cache, an adaptive integral over an infinite range, and a
|
|||
|
|
// two-dimensional integral by globally adaptive bisection.
|
|||
|
|
func Example_quadratureAndCubature() {
|
|||
|
|
nodes, weights, err := integrate.GaussLegendreNodes(3)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
for i := range nodes {
|
|||
|
|
fmt.Printf("node %.6f, weight %.6f\n", nodes[i], weights[i])
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
value, errEst, err := integrate.IntegrateFunction(func(x float64) (float64, error) {
|
|||
|
|
return math.Exp(-x * x), nil
|
|||
|
|
}, 0, math.Inf(1), integrate.QuadratureOptions{})
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Printf("the Gaussian tail integrates to %.6f (error estimate %.1e)\n", value, errEst)
|
|||
|
|
|
|||
|
|
area, err := integrate.IntegrateND(func(x []float64) float64 {
|
|||
|
|
return x[0] * x[1]
|
|||
|
|
}, []float64{0, 0}, []float64{1, 1}, integrate.CubatureOptions{})
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Printf("x·y over the unit square integrates to %.6f\n", area)
|
|||
|
|
// Output:
|
|||
|
|
// node -0.774597, weight 0.555556
|
|||
|
|
// node 0.000000, weight 0.888889
|
|||
|
|
// node 0.774597, weight 0.555556
|
|||
|
|
// the Gaussian tail integrates to 0.886227 (error estimate 9.8e-12)
|
|||
|
|
// x·y over the unit square integrates to 0.250000
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Heat evolution in one dimension: u_t = u_xx on [0, 1] from
|
|||
|
|
// u = sin(πx), Dirichlet ends held at zero. The sampled history is a
|
|||
|
|
// (samples, n) array of interior states, and the centre decays as the
|
|||
|
|
// exact e^(−π²t)·sin(π/2) predicts.
|
|||
|
|
func ExampleIntegrateHeat1D() {
|
|||
|
|
const n, samples = 399, 5
|
|||
|
|
const dx, tFinal = 1.0 / 400, 0.1
|
|||
|
|
u0 := make([]float64, n)
|
|||
|
|
for i := range u0 {
|
|||
|
|
u0[i] = math.Sin(math.Pi * float64(i+1) * dx)
|
|||
|
|
}
|
|||
|
|
state, err := tensor.FromFloats(u0, n)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
history, err := integrate.IntegrateHeat1D(state, 1, dx, tFinal, 1e-4, samples, 0, 0)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
centre := history.FloatAt((samples-1)*n + n/2)
|
|||
|
|
exact := math.Exp(-math.Pi * math.Pi * tFinal)
|
|||
|
|
fmt.Printf("history shape %v\n", history.Shape())
|
|||
|
|
fmt.Printf("u(1/2, 0.1) = %.6f, exact %.6f\n", centre, exact)
|
|||
|
|
// Output:
|
|||
|
|
// history shape [5 399]
|
|||
|
|
// u(1/2, 0.1) = 0.372710, exact 0.372708
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// The finite element Poisson solve: −∇·(κ∇u) = f on the unit square
|
|||
|
|
// with κ = 1 and Dirichlet data on the boundary ring. The manufactured
|
|||
|
|
// solution u = sin(πx)·sin(πy) makes f = 2π²·sin(πx)·sin(πy), and the
|
|||
|
|
// P1 solution reproduces it to the mesh's accuracy at the centre.
|
|||
|
|
func ExampleSolvePoissonFEM2D() {
|
|||
|
|
const cells = 16
|
|||
|
|
mesh, err := integrate.GridTriangleMesh2D(0, 0, 1, 1, cells, cells)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
var nodes []int
|
|||
|
|
var values []float64
|
|||
|
|
for v := range mesh.Vertices2() {
|
|||
|
|
x, y := mesh.Vertices[2*v], mesh.Vertices[2*v+1]
|
|||
|
|
onEdge := x == 0 || x == 1 || y == 0 || y == 1
|
|||
|
|
if onEdge {
|
|||
|
|
nodes = append(nodes, v)
|
|||
|
|
values = append(values, math.Sin(math.Pi*x)*math.Sin(math.Pi*y))
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
u, err := integrate.SolvePoissonFEM2D(mesh, func(x, y float64) float64 {
|
|||
|
|
return 2 * math.Pi * math.Pi * math.Sin(math.Pi*x) * math.Sin(math.Pi*y)
|
|||
|
|
}, integrate.FEMPoissonOptions{Kappa: 1, DirichletNodes: nodes, DirichletValues: values})
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
centre := (cells/2)*(cells+1) + cells/2
|
|||
|
|
fmt.Printf("mesh of %d vertices, %d triangles, %d boundary edges\n",
|
|||
|
|
mesh.Vertices2(), mesh.Triangles3(), len(mesh.BoundaryEdges())/2)
|
|||
|
|
fmt.Printf("u(1/2, 1/2) = %.4f on this mesh, exact 1.0000\n", u.FloatAt(centre))
|
|||
|
|
// Output:
|
|||
|
|
// mesh of 289 vertices, 512 triangles, 64 boundary edges
|
|||
|
|
// u(1/2, 1/2) = 0.9946 on this mesh, exact 1.0000
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// The two-point boundary value problem: y″ = −y with y(0) = 0 and
|
|||
|
|
// y(π/2) = 1, solved by shooting on the free initial slope. The slope
|
|||
|
|
// comes out as 1 and the sampled trajectory traces y = sin t.
|
|||
|
|
func ExampleIntegrateBoundary() {
|
|||
|
|
f := func(t float64, y *tensor.Array) (*tensor.Array, error) {
|
|||
|
|
return tensor.FromFloats([]float64{y.FloatAt(1), -y.FloatAt(0)}, 2)
|
|||
|
|
}
|
|||
|
|
y0, _ := tensor.FromFloats([]float64{0, 0.5}, 2)
|
|||
|
|
bc := integrate.BoundaryConditions{Start: []int{0}, End: []int{0}, EndValues: []float64{1}}
|
|||
|
|
times, states, err := integrate.IntegrateBoundary(f, 0, math.Pi/2, y0, bc, 3,
|
|||
|
|
integrate.ODEOptions{RelTol: 1e-10, AbsTol: 1e-13})
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Printf("shooting slope y'(0) = %.6f\n", states[0].FloatAt(1))
|
|||
|
|
for i := range times {
|
|||
|
|
fmt.Printf("y(%.4f) = %.6f, exact %.6f\n", times[i], states[i].FloatAt(0), math.Sin(times[i]))
|
|||
|
|
}
|
|||
|
|
// Output:
|
|||
|
|
// shooting slope y'(0) = 1.000000
|
|||
|
|
// y(0.0000) = 0.000000, exact 0.000000
|
|||
|
|
// y(0.7854) = 0.707107, exact 0.707107
|
|||
|
|
// y(1.5708) = 1.000000, exact 1.000000
|
|||
|
|
}
|