// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

// Package integrate solves differential equations, integrates functions
// and evolves partial differential equations. It carries five families:
// ordinary differential equations as initial value problems, with event
// detection on the way; two-point boundary value problems by shooting
// and by collocation, and Hamiltonian systems by symplectic schemes;
// adaptive quadrature and cubature; turnkey heat, wave and advection
// solvers in one and two space dimensions; and a piecewise-linear
// finite element Poisson solver on triangular and tetrahedral meshes.
//
// # The state contract
//
// An ordinary differential equation is y' = f(t, y), where f returns the
// derivative of the state at a time and a state. The state is a rank-1
// array: a system of higher rank flattens to its leading-axis vector
// first. The ODE family reads its elements one by one and widens them to
// float64, so an int or float32 state integrates there; the symplectic
// family accepts float64 and float32 positions and momenta only, and a
// complex state is refused everywhere. The returned trajectories are
// float64 arrays, freshly allocated, and the inputs are never written
// to.
//
// Every solver refuses rather than guesses. An exhausted step budget, a
// step size that has collapsed below the resolution of t, an f that
// returns a wrongly shaped state, a non-finite value or a violated CFL
// budget under an explicit stencil is an error naming itself, never a
// silently truncated or silently wrong trajectory.
//
// # Initial value problems
//
// IntegrateODE is the adaptive workhorse (an embedded Dormand-Prince
// 4(5) pair); IntegrateRK4 is the classical fixed-step scheme;
// IntegrateBackwardEuler, IntegrateBDF2, IntegrateBDFVar and
// IntegrateROS4 cover the stiff regime, from the entry-level implicit
// Euler to variable-order BDF and an L-stable Rosenbrock-Wanner method.
// IntegrateODEPath samples the trajectory on an even time grid,
// IntegrateODESteps records every accepted step, and IntegrateODEEvents
// additionally reports where a list of watches crosses zero, filtered by
// direction. IntegrateDAE takes the semi-explicit index-1 mass-matrix
// form M·y' = f(t, y). Backward integration works throughout: a t1 < t0
// integrates in the negative direction.
//
// # Beyond the initial value problem
//
// IntegrateBoundary shoots a two-point boundary value problem, choosing
// the free initial components so the trajectory lands on the prescribed
// end values; SolveBoundaryCollocation solves the same problem by
// three-point Lobatto IIIA collocation on an adaptively refined mesh and
// returns the mesh with the nodal states and slopes. IntegrateVerlet and
// IntegrateYoshida4 integrate a separable Hamiltonian system at a fixed
// step, and IntegrateMidpoint does the same for a general, non-separable
// one through an implicit stage.
//
// # Quadrature and cubature
//
// IntegrateFunction integrates a scalar function over a finite or
// infinite interval and reports an estimate of its own absolute error;
// GaussLegendreNodes hands out the nodes and weights of a fixed rule.
// IntegrateFilon integrates a smooth amplitude against a high-frequency
// sine or cosine carrier, whose cost tracks the amplitude alone rather
// than the carrier a sampled rule must resolve. IntegrateND integrates
// over a hyperrectangle by globally adaptive bisection with product
// Gauss-Legendre rules.
//
// # PDE evolution
//
// IntegrateHeat1D, IntegrateWave1D, IntegrateUpwindAdvection1D,
// IntegrateAdvection1D and IntegrateAdvectionDiffusion1D run on the
// interior grid of a rank-1 initial state, while IntegrateHeat2D and
// IntegrateWave2D run on the rank-2 grid of a rectangle. All of them
// return the trajectory sampled on a time grid, endpoints included.
//
// # Finite elements
//
// GridTriangleMesh2D and BoxTetraMesh3D build structured meshes;
// NewTriangleMesh2D and NewTetraMesh3D accept general conforming ones,
// refusing degenerate elements. SolvePoissonFEM2D and SolvePoissonFEM3D
// assemble and solve -∇·(κ∇u) = f with P1 elements, a conductivity that
// may vary in space, Dirichlet values eliminated by lifting and Neumann
// fluxes integrated on prescribed boundary edges or faces.
//
// # What it deliberately does not do
//
// There is no dense-output object: IntegrateODEPath and
// IntegrateODESteps return the states a caller asked for, and event
// times are narrowed by re-integrating the accepted step rather than
// through a continuous extension. The symplectic family takes a fixed
// step by design, because adaptivity would destroy the property the
// methods exist for. IntegrateDAE is first order and does not project an
// inconsistent start onto the constraint manifold; consistent initial
// values are the caller's contract. The finite element surface is P1 on
// conforming meshes only, and the collocation solver factors a dense
// Newton matrix, so its mesh size is bounded by CollocationOptions. The
// gradient of a trajectory with respect to its parameters is the grad
// package's to compute.
//
// A tour:
//
//	end, _ := integrate.IntegrateODE(f, 0, 1, y0, integrate.ODEOptions{})
//	hits, end, _ := integrate.IntegrateODEEvents(f, 0, 5, y0, watches, integrate.ODEOptions{})
//	area, _ := integrate.IntegrateFunction(g, 0, 1, integrate.QuadratureOptions{})
//	history, _ := integrate.IntegrateHeat1D(u0, 1, dx, 0.1, 1e-4, 5, 0, 0)
//	u, _ := integrate.SolvePoissonFEM2D(mesh, f, opts) // opts carries κ and the Dirichlet set
//
// The examples in this documentation are executable and checked by the
// test suite.
package integrate