NumericsLib

Numerical methods domain for Free Pascal — root finding, quadrature, ODE solvers, and interpolation.

For the 1.7 adaptive, vector, fitting, derivative, and reusable-interpolant APIs, see the numerical modelling guide. The APIs on this page remain the concise compatibility and teaching path.

Depends on: MathBase

Learning routes

Beginner route

Copy and run the double-real numerical quick start. It prints checked root, integral, ODE, and spline values using plain callbacks and allocating result records. Start with Brent for a bracketed scalar root.

Common tasks and algorithm choice

Advanced route

Run example 17 for reusable interpolants, fitting, bounded options, event-aware adaptive ODEs, and explicit statuses. It retains the double-real array path; coefficient/result arrays are allocated explicitly and no private dense conversion is needed.

Units

All methods are class-static — no instantiation required:

Root := TNumericsKit.Brent(@F, 1.0, 2.0);

---

Function Types

TScalarFunc = function(X: Double): Double;        // f(x) — used by root finders and integrators

TODEFunc    = function(T, Y: Double): Double;     // dy/dt = f(t,y) — used by ODE solvers

Result Records

TODESolution = record

  T: TDoubleArray;   // independent variable values (length N+1)

  Y: TDoubleArray;   // solution values             (length N+1)

end;



TCubicSpline = record

  X: TDoubleArray;   // knot x-values (sorted ascending)

  A, B, C, D: TDoubleArray;  // polynomial coefficients per interval

end;

---

Root Finding

All methods find x such that f(x) ≈ 0. Tolerance parameters control convergence. Each scalar-returning method has a corresponding ...Result overload returning:

TRootResult = record

  Root: Double;

  Residual: Double;

  Iterations: Integer;

  Converged: Boolean;

  Status: TIterationStatus;

end;

Use the detailed result when iteration counts or explicit convergence reporting matter. Scalar methods raise ENumericsConvergenceError when their iteration limit is exhausted. Status distinguishes convergence, numerical breakdown, and iteration exhaustion; Converged remains for source compatibility.

Bisection

class function Bisection(F: TScalarFunc; A, B: Double;

  Tol: Double = 1E-10; MaxIter: Integer = 100): Double;

class function BisectionResult(F: TScalarFunc; A, B: Double;

  Tol: Double = 1E-10; MaxIter: Integer = 100): TRootResult;

Guaranteed convergence for bracketed roots. Requires opposite endpoint signs, unless either endpoint is itself a root, in which case it is returned immediately. Raises EInvalidArgument if the bracket is invalid.

Best for: when robustness matters more than speed; verifying root existence.

NewtonRaphson

class function NewtonRaphson(F, DF: TScalarFunc; X0: Double;

  Tol: Double = 1E-10; MaxIter: Integer = 100): Double;

class function NewtonRaphsonResult(F, DF: TScalarFunc; X0: Double;

  Tol: Double = 1E-10; MaxIter: Integer = 100): TRootResult;

Quadratic convergence near the root using f and its derivative df/dx. Raises EInvalidArgument if df/dx is near zero (flat region).

Best for: smooth functions where the derivative is cheap to evaluate.

Brent

class function Brent(F: TScalarFunc; A, B: Double;

  Tol: Double = 1E-10; MaxIter: Integer = 100): Double;

class function BrentResult(F: TScalarFunc; A, B: Double;

  Tol: Double = 1E-10; MaxIter: Integer = 100): TRootResult;

Hybrid method combining bisection, secant, and inverse-quadratic interpolation. Superlinear convergence in practice; falls back to bisection when necessary. Raises EInvalidArgument if f(A) · f(B) > 0.

Best for: general use — the recommended default solver.

Secant

class function Secant(F: TScalarFunc; X0, X1: Double;

  Tol: Double = 1E-10; MaxIter: Integer = 100): Double;

class function SecantResult(F: TScalarFunc; X0, X1: Double;

  Tol: Double = 1E-10; MaxIter: Integer = 100): TRootResult;

Derivative-free quasi-Newton using two initial guesses. Superlinear convergence but may diverge without a good bracket.

Best for: when a derivative is unavailable and you have a good initial estimate.

---

Numerical Integration

TrapezoidalRule

class function TrapezoidalRule(F: TScalarFunc; A, B: Double;

  N: Integer = 1000): Double;

Composite trapezoidal rule with N sub-intervals. Error O(h²). Exact for linear functions.

SimpsonRule

class function SimpsonRule(F: TScalarFunc; A, B: Double;

  N: Integer = 1000): Double;

Composite Simpson's 1/3 rule with N sub-intervals. Values below 2 are rejected, and an odd N is incremented to the next even value. Error O(h⁴). Exact for polynomials of degree ≤ 3.

GaussLegendre5

class function GaussLegendre5(F: TScalarFunc; A, B: Double): Double;

5-point Gauss-Legendre quadrature on [A, B]. Exact for polynomials of degree ≤ 9. Very accurate for smooth functions with only 5 function evaluations.

Comparison:

---

ODE Solvers

Solve dy/dt = f(t, y), y(t₀) = y₀ over the interval [T0, T1].

Single Steps

class function EulerStep(F: TODEFunc; T0, Y0, H: Double): Double;

class function RK4Step(F: TODEFunc; T0, Y0, H: Double): Double;

Return the solution at T0 + H given the current state (T0, Y0).

Full Solvers

class function EulerSolve(F: TODEFunc; T0, Y0, T1: Double; N: Integer): TODESolution;

class function RK4Solve(F: TODEFunc; T0, Y0, T1: Double; N: Integer): TODESolution;

Integrate from T0 to T1 using N uniform steps of size h = (T1−T0)/N. Return a TODESolution with T and Y arrays of length N+1.

Accuracy comparison for dy/dt = y, y(0) = 1 over [0, 1] (exact: e ≈ 2.71828):

---

Interpolation

LinearInterp

class function LinearInterp(const XKnots, YKnots: TDoubleArray; X: Double): Double;

Piecewise linear interpolation between sorted knots using binary search. Clamps to endpoint values outside the knot range.

XKnots must be strictly increasing and the two arrays must have equal, non-zero length. Array lengths, finite values, ordering, and duplicate knots are validated.

LagrangeInterp

class function LagrangeInterp(const XKnots, YKnots: TDoubleArray; X: Double): Double;

Global Lagrange polynomial interpolation through all N knots. Exact at every knot.

The arrays must have equal, non-zero length and finite values, and X knots must be distinct. These conditions are validated before interpolation.

> Warning: ill-conditioned for N > ~10 (Runge phenomenon). Prefer CubicSplineBuild for larger datasets.

CubicSplineBuild / CubicSplineEval

class function CubicSplineBuild(const XKnots, YKnots: TDoubleArray): TCubicSpline;

class function CubicSplineEval(const S: TCubicSpline; X: Double): Double;

Natural cubic spline (zero second-derivative boundary conditions) solved via the Thomas tridiagonal algorithm. Exact at every knot; smooth C² between knots. Clamps to endpoint values outside the knot range.

Spline construction requires equal-length arrays with at least two strictly increasing finite knots. Dimensions, ordering, duplicates, and finite values are validated during construction; evaluation also validates spline dimensions and its finite query value.

---

Quick Start

uses MathBase.SharedTypes, NumericsLib.Numerics;



{ --- functions for root finding and integration --- }

function F(X: Double): Double; begin Result := X*X - 2; end;

function DF(X: Double): Double; begin Result := 2*X; end;

function G(X: Double): Double; begin Result := X*X; end;



{ --- ODE: dy/dt = y  →  exact solution y(t) = e^t --- }

function DYDT(T, Y: Double): Double; begin Result := Y; end;



var

  Root, Integral: Double;

  Sol: TODESolution;

  XK, YK: TDoubleArray;

  Spline: TCubicSpline;

begin

  { Root: Brent's method — recommended default }

  Root := TNumericsKit.Brent(@F, 1.0, 2.0);

  Writeln('sqrt(2) = ', Root:0:10);      // 1.4142135624



  { Root: Newton-Raphson — fast with derivative }

  Root := TNumericsKit.NewtonRaphson(@F, @DF, 1.5);

  Writeln('sqrt(2) = ', Root:0:10);



  { Integration: ∫₀¹ x² dx = 0.3333... }

  Integral := TNumericsKit.SimpsonRule(@G, 0, 1, 1000);

  Writeln('Integral = ', Integral:0:6);  // 0.333333



  { ODE: RK4 is accurate with far fewer steps than Euler }

  Sol := TNumericsKit.RK4Solve(@DYDT, 0, 1.0, 1.0, 100);

  Writeln('y(1) = e = ', Sol.Y[100]:0:8);  // 2.71828183



  { Spline interpolation through y = x² at integer knots }

  XK := TDoubleArray.Create(0, 1, 2, 3, 4);

  YK := TDoubleArray.Create(0, 1, 4, 9, 16);

  Spline := TNumericsKit.CubicSplineBuild(XK, YK);

  Writeln('Spline(1.5) = ', TNumericsKit.CubicSplineEval(Spline, 1.5):0:4);  // ≈ 2.2321

end.

Expected output:

sqrt(2) = 1.4142135624

sqrt(2) = 1.4142135624

Integral = 0.333333

y(1) = e = 2.71828183

Spline(1.5) = 2.2321

---

Error Handling

---

Design Notes

• All functions are class-static — pass function pointers, not method pointers. In FPC, unit-level function declarations are compatible with TScalarFunc / TODEFunc.
• The cubic spline uses natural boundary conditions (S''(x₀) = S''(xₙ) = 0). This is the standard choice when no derivative information is available at the boundaries.
• GaussLegendre5 performs only 5 function evaluations regardless of the smoothness of f. For oscillatory functions, use SimpsonRule with a large N instead.
• TODESolution arrays are zero-indexed: Sol.T[0] = T0, Sol.T[N] = T1.
• Root ...Result methods expose root, residual, iteration count, and a convergence flag. Scalar wrappers raise ENumericsConvergenceError on exhaustion.
• Root, integration, ODE, and interpolation entry points validate callbacks, finite inputs and callback outputs, positive controls, array dimensions, and strictly increasing/distinct interpolation knots as applicable.
• Bracketed solvers return endpoint roots immediately and reject brackets that contain no sign change.

Common mistakes

• Bracketing requires a sign change. Bisection and Brent reject brackets where f(A) and f(B) share a sign with EInvalidArgument.
• Interpolation is not fitting. Interpolators pass exactly through every knot; a fitted curve does not. LagrangeInterp is ill-conditioned beyond about ten knots (Runge phenomenon) — prefer CubicSplineBuild.
• The ODE solvers are explicit and non-stiff. EulerSolve and RK4Solve suit ordinary initial-value problems; stiff systems are a post-2.0 roadmap gap and are not supported here.
• Fixed-evaluation quadrature. GaussLegendre5 uses only five function evaluations; for oscillatory integrands use SimpsonRule with a large N.