🥮 mooncakes.io

    numsolve

    Pure MoonBit ODE suite with dense output, events, direct/iterative/banded linear systems, root finding, and interpolation

    root-finding
    numerical-analysis
    differential-equations
    ode
    runge-kutta
    solve-ivp
    linear-algebra
    linear-solvers
    iterative-solvers
    banded-solvers
    dense-output
    event-detection
    interpolation
    spline
    scientific-computing
    Download zip
    Version
    0.6.0
    License
    Apache-2.0
    Repository
    Last updated
    15 days ago
    Downloads
    43

    No dependencies

    #numsolve

    Numerical solvers in pure MoonBit: root finding, ordinary differential equations, linear systems, and interpolation. The algorithms are standard published methods, implemented as a small MoonBit-first toolkit; selected results are checked against analytic solutions or SciPy where that is a meaningful reference.

    This is not presented as a new algorithm or a drop-in replacement for mature numerical ecosystems. Its value is a coherent, tested MoonBit API that covers common scalar and small-system problems without a foreign-language runtime.

    #Installation

    moon add noweizongqing/numsolve

    The package is published on mooncakes.io. The ODE layer supports both the original scalar convenience functions and vector systems through solve_ivp, with cubic Hermite dense output between accepted steps and first-event detection on that interpolation. Implicit or stiff solvers are not implemented yet and are called out below.

    ///|
    let root = @numsolve.brent(fn(x) { x * x - 2.0 }, 1.0, 2.0) // 1.4142135623731364

    ///|
    let faster = @numsolve.newton(fn(x) { x * x - 2.0 }, fn(x) { 2.0 * x }, 1.5) // 1.4142135623730951, in four steps

    #Scope

    • The package favours correctness, explicit error values, and small dependency surfaces over matching the full feature set of SciPy, Sundials, or DifferentialEquations.jl.
    • ODE support now includes vector-valued states, relative and absolute tolerances, a common solve_ivp dispatch, cubic Hermite dense output, and first-event detection with direction filtering. It does not yet provide sensitivity analysis or implicit stiff solvers.
    • Reference comparisons are part of the test suite where they add information; they are not a claim that every routine is a line-for-line port.

    #Root finding

    roots.mbt finds where a function of one variable is zero. Three of the five searches keep a bracket — two points whose function values have opposite signs, which promises a root between them — and two are unbracketed, starting from a point instead.

    FunctionMethodSteps to 1e-10 on x² − 2 from [1, 2]
    bisectionhalve the interval, keep the half the sign change is in32
    false_positionthe line through the two ends of the bracket, with the Illinois halving8
    secantthe line through the last two points, unbracketed6
    newtonthe tangent at the current point, derivative needed4
    brentbisection, the secant and inverse quadratic interpolation, chosen between6

    Bisection is the slow one and the only one whose accuracy can be written down before it runs. False position is bisection's bracket with a much better point than the midpoint, and needs the Illinois modification to be usable at all: on x¹⁰ − 1, which curves the same way across the whole bracket, the plain method is two thousandths out after a hundred steps while the halving version is exact in seventeen. The secant gets speed without a derivative, at the cost of the bracket. Newton is the fastest and the most demanding. Brent is the one to reach for by default: it keeps bisection's bracket, converges as fast as the secant on smooth functions, and falls back to a bisection step whenever a fitted one looks like a bad idea.

    ///|
    let cubic = @numsolve.brent(fn(x) { x * x * x - x - 2.0 }, 1.0, 2.0)

    ///|
    let sine = @numsolve.false_position(fn(x) { @math.sin(x) - 0.5 }, 0.0, 1.0)

    The tolerance is measured in the variable being solved for, and each method reads it in the way its own search can honour:

    MethodStops when
    bisection, false_position, brentthe interval left is no wider than tolerance · (1 + \|x\|)
    secant, newtonthe step is that small, or the function value is no larger than tolerance

    The (1 + |x|) factor is what makes one tolerance mean the same thing for a root near zero and a root a long way from it. A tolerance of exactly zero turns the tolerance test off, and the search then runs the full max_iterations steps and returns what it has — that is how a fixed number of steps is asked for, and what the convergence tests are built on. A function value that comes out exactly zero ends a search at once, at any tolerance, and an endpoint that is already a root is not searched for at all.

    What a search cannot do, or cannot finish, is raised rather than answered around:

    VariantPayloadRaised when
    NotBracketed(Double, Double)f(a) and f(b)the two points given as a bracket have the same sign
    ZeroDerivative(Double)the point where it happeneda step would divide by a slope of exactly zero
    NonPositiveTolerance(Double)the tolerancethe tolerance is negative
    MaxIterations(Int)the limitthe search took max_iterations steps without meeting its tolerance

    ///|
    let root = @numsolve.bisection(fn(x) { x * x - 2.0 }, 1.0, 2.0) catch {
      _ => 0.0
    }

    A bracket may be given either way round, the two points of a bracket that hold opposite signs promise a root and anything else is refused with both values, and an unbracketed search — which promises nothing — reports a stall instead: the secant on x¹⁰ − 1 from 0 and 1.5 walks off to 1.9e11 and back, lands on a point it repeats exactly, and reports the zero slope it is left with. Newton on x³ − 2x + 2 from zero steps to one and back to zero for ever, and reports running out of iterations.

    The tests measure the convergence each algorithm claims, not only its final answer:

    ClaimMeasured
    Bisection after k steps is within (b − a) / 2^k of the rootheld for k = 1, 2, 4, 8, 12, 16, 20
    Newton's error is about the square of the error before it8.6e-2, 2.5e-3, 2.1e-6, 1.6e-12, 0
    The secant's error falls as the golden ratio's power, 1.618log-ratios of 1.83, 1.68 and 1.67 for the last steps
    The Illinois halving is what keeps false position off an endexact in 17 steps where the plain method is 2.2e-3 out at 100

    #Ordinary differential equations

    ode.mbt keeps the original scalar API for y'(t) = f(t, y): euler, midpoint, rk4, and adaptive Dormand-Prince rk45. ode_system.mbt adds the form needed by most real models, where the state is an Array[Double]. Both layers return a trajectory containing the times the solver stopped at and the value or state at each one. Integration backwards is allowed, and the steps then have the sign of t1 - t0.

    ///|
    let doubling = @numsolve.rk4(fn(_t, y) { y }, 1.0, 0.0, 1.0, 10) // 2.718279744135166

    ///|
    // The harmonic oscillator y'' = -y, written as [position, velocity]' = [v, -y].
    let oscillator = @numsolve.solve_ivp(
      fn(_t, y) { [y[1], -y[0]] },
      [1.0, 0.0],
      0.0,
      6.283185307179586,
      solver=@numsolve.OdeMethod::RK4,
      steps=2000,
    )

    APIStateStep selection
    euler, midpoint, rk4scalar Doublefixed number of equal steps
    rk45scalar Doubleadaptive Dormand-Prince 5(4)
    *_system counterpartsArray[Double]fixed or adaptive, as above
    solve_ivpArray[Double]dispatches to Euler, Midpoint, RK4, or RK45
    first_eventVectorTrajectorylocates the first sign change and its time/state
    truncate_atVectorTrajectoryreturns the prefix ending exactly at a requested time

    The scalar fixed-step methods remain the simplest entry point. solve_ivp is the recommended entry point for new code because it gives systems and scalar solvers a common shape without hiding the method choice. It is intentionally a small, familiar API, not a claim to cover SciPy's full solve_ivp feature set.

    A VectorTrajectory keeps the state and its derivative at every accepted sample. interpolate_at(t) evaluates a cubic Hermite segment between two samples, and sample_at(times) does the same for a list of requested times. The sample times may be unsorted, and a time outside the trajectory raises OutOfRange instead of silently extrapolating.

    first_event(trajectory, event, direction=...) scans the accepted samples for a sign change in the requested direction and refines the crossing with Brent's method on the same Hermite interpolation. The result carries the event time, interpolated state, and interval index; no extra integration step is taken. truncate_at(event.time()) returns a new trajectory whose last sample is the interpolated event state, which gives a terminal-event result without mutating the original trajectory.

    StepperOrderEvaluations per stepError when the step halves
    euler11÷ 2
    midpoint22÷ 4
    rk444÷ 16
    rk455, adaptive7set by the tolerances

    On y' = y from 1 over [0, 1], where the answer is e, ten steps of each fixed-step method give:

    StepperTen stepsError
    euler2.59374246011.2e-1
    midpoint2.7140808466082244.2e-3
    rk42.7182797441351662.1e-6
    rk45 at 1e-82.71828183620885347.8e-9, in 11 steps and 77 evaluations

    rk45 uses seven evaluations to produce both a fifth-order and a fourth-order answer; their difference estimates the error of the step. A step is accepted only when that estimate meets the tolerance, and the next step is sized from the same estimate. The scalar controller uses |y5 - y4| / (1 + |y5|). The system controller uses the root-mean-square of the component errors after scaling each component by both tolerances: atol + rtol * max(|y_old|, |y_new|). This is why a vector problem with one large and one near-zero component does not let either component dominate the step decision.

    tolerance, rtol, and atol must be positive. The trajectory that comes back is sampled where the adaptive solver chose to stop, not at an evenly spaced grid. max_evaluations is a budget on calls to the right-hand side rather than on steps; every adaptive step costs seven evaluations. Running out of that budget raises instead of returning a trajectory that stops short of t1.

    VariantPayloadRaised when
    EmptyStatenonea vector system was given an empty initial state
    DimensionMismatch(Int, Int)state dimension, derivative dimensionthe right-hand side returned an array whose length does not match the state
    NonPositiveSteps(Int)the step counta fixed-step stepper was given zero steps or fewer
    NonPositiveTolerance(Double)the offending tolerancean adaptive solver was given a tolerance that is not positive
    MaxEvaluations(Int)the evaluation budgetan adaptive solver spent its budget without reaching t1
    OutOfRange(Double, Double, Double)start, end, requested timea dense-output sample lies outside the trajectory

    The tests assert the order each fixed-step method claims, the adaptive solver's tolerance behaviour, a two-dimensional oscillator with a known period, backward integration, dimension validation, and empty-state validation. The ODE tests also cover cubic Hermite sampling on fixed and adaptive trajectories, out-of-range errors, event detection on scalar and adaptive vector systems, and terminal trajectory truncation. Stiff/implicit methods are not claimed yet.

    A longer design note, including the ecosystem comparison and roadmap, is in docs/ode.md.

    #Linear systems

    linear.mbt solves A · x = b for a dense square matrix, given as an Array[Array[Double]] of rows.

    ///|
    let a = [[2.0, 1.0, 1.0], [4.0, 3.0, 3.0], [8.0, 7.0, 9.0]]

    ///|
    let x = @numsolve.solve_linear_system(a, [1.0, 2.0, 4.0]) // [0.5, 0.0, 0.0]

    FunctionMethodUse it when
    LU::decompose and solve / solve_many / solve_transposedGaussian elimination with partial pivoting, P · A = L · Uthe matrix is general, and one factorization is reused for several right-hand sides or the transpose
    Cholesky::decompose and solve / solve_manyA = L · Lᵀ, no pivotingthe matrix is symmetric positive definite: half the work of an LU
    solve_tridiagonalthe Thomas recurrence over three diagonalsthe matrix is tridiagonal — O(n) instead of O(n³)
    solve_bandedbanded elimination with pivoting inside the lower bandwidththe matrix has a declared p/q bandwidth — O(n·(p+q)²) instead of O(n³)
    residual_normcompensated residual ‖b − A·x‖∞a computed solution has to be judged, especially when A is ill-conditioned
    LU::solve_refined / solve_linear_system_refinedresidual correction passes over one factorizationa direct solve needs another pass without refactoring
    jacobi / gauss_seidelstationary iterations using only matrix-vector productsthe system is large or diagonally dominant and a factorization is not wanted
    conjugate_gradient / preconditioned_conjugate_gradientKrylov iterations for symmetric positive definite systemsthe matrix is SPD and a factorization is not wanted

    A decomposition is worth keeping when there is more than one right-hand side: LU::solve_many answers with a list of them after one factorization, and LU::inverse is that call with the identity. LU::solve_transposed answers Aᵀ · x = b from the factors of A instead of factoring the transpose, and Cholesky::solve_many is the plural form of the symmetric solve. LU::determinant and Cholesky::determinant come from the diagonal of the factorization at no extra cost.

    solve_banded generalizes the structured solve to p sub-diagonals and q super-diagonals. It validates that the dense input has no non-zero entry outside the declared band, then limits pivot search and elimination to the band and its fill-in. A matrix that violates its declared structure raises NotBanded instead of being solved with the wrong assumption.

    jacobi and gauss_seidel are the stationary iterative alternatives: they use only matrix-vector products, need no factorization, and can stop as soon as the residual meets a tolerance. They are not unconditional solvers: the diagonal must be non-zero, the iteration has to converge, and a run that exhausts its budget raises NoConvergence. On the same diagonally dominant system, Gauss-Seidel usually takes fewer iterations because it uses each component as soon as it is updated.

    conjugate_gradient adds a Krylov method for symmetric positive definite systems. It also uses matrix-vector products, but chooses each search direction so that the residual decreases in the A-norm. The Jacobi-preconditioned variant divides by the diagonal before each inner product, which often reduces the iteration count when the diagonal varies widely. A non-positive curvature means the matrix is not SPD or the iteration has broken down, and is raised as Breakdown.

    The residual is computed with a compensated dot product rather than a plain matrix-vector multiply. That matters for the last few digits: on Hilbert 8 the plain residual is already close to the rounding floor, and a naive residual would be rounded to zero before it could be used as a correction. solve_refined uses that residual to solve A·δ = r with the existing factors and returns x + δ. It can reduce the residual further, but it does not promise a fixed number of correct digits in x: the forward error of a very ill-conditioned system remains bounded by its condition number and by the precision of the arithmetic.

    Pivoting is not an optimization here but the difference between an answer and a wrong one. On [[0, 2], [3, 4]] there is no answer at all without a row exchange — the pivot is zero — and on [[1e-18, 1], [1, 1]], whose exact solution is [1, 1], elimination in the order the rows come in divides by 1e-18 and answers [0, 1]. Both cases are in the tests, the second one against an elimination written without the pivot search.

    What the arithmetic costs shows up in the condition number. The Hilbert matrix of order 8 — H[i][j] = 1 / (i + j + 1), condition number 1.5e10 — solved for the right-hand side made of its row sums has the exact solution of all ones, and comes back with a largest error of 1.4e-7: about ten of the sixteen digits a double has are left, and that is the matrix rather than the method. The test asserts 1e-5 rather than that figure, so that a different order of summation inside the same algorithm cannot fail it.

    VariantPayloadsRaised when
    NotSquare(Int, Int)rows, columnsthe matrix is not square, or its rows are not all the same length
    SizeMismatch(Int, Int)needed, givena right-hand side or a diagonal is the wrong length
    Singular(Int)the columna pivot is exactly zero and no row below it can be exchanged for it
    NotPositiveDefinite(Int)the indexCholesky was given a matrix whose pivot is not positive
    NonPositiveRefinementSteps(Int)the step countiterative refinement was asked for zero passes or fewer
    NonPositiveTolerance(Double)the tolerancea tolerance for refinement or an iterative solve is negative
    NonPositiveMaxIterations(Int)the limitan iterative solve was given zero iterations or fewer
    ZeroDiagonal(Int)the rowJacobi or Gauss-Seidel met a zero diagonal entry
    NoConvergence(Int, Double)iterations, residual norman iterative solve exhausted its budget without meeting the tolerance
    Breakdown(Int)the iterationconjugate gradient met a non-positive pᵀ·A·p, so the matrix is not positive definite or the iteration broke down
    InvalidBandwidth(Int)the bandwidththe banded solver was given a negative bandwidth
    NotBanded(Int, Int)row, columnthe dense matrix has a non-zero entry outside its declared band

    #Relationship to mizchi/numbt

    mizchi/numbt already provides native BLAS/LAPACK-backed dense solves through fmat_solve and related functions. numsolve overlaps with that basic operation, and says so explicitly rather than claiming a new dense solver. The extension boundary is different: numsolve is pure MoonBit, uses Double, is designed for portable/Wasm use, returns reusable LU and Cholesky factors, adds multiple-right-hand-side and transpose solves, provides compensated residual diagnostics and iterative refinement, adds Jacobi/Gauss-Seidel iterations, and keeps the O(n) tridiagonal path that the spline module depends on. The full comparison, interoperation guidance, and non-goals are in docs/linear.md.

    #Interpolation

    interp.mbt draws a curve through a table of points — two arrays, xs and ys, of the same length, with the xs in increasing order — and reads the value it takes between them.

    ///|
    let table_x = [0.0, 1.0, 2.0, 3.0]

    ///|
    let table_y = [1.0, 2.0, 0.0, 4.0]

    ///|
    let line = @numsolve.LinearInterp::new(table_x, table_y)

    ///|
    let spline = @numsolve.CubicSpline::new(
      table_x,
      table_y,
      @numsolve.SplineBoundary::Natural,
    )

    TypeThe curveUse it when
    LinearInterpthe straight line between the two knots around xthe table is dense, or a kink at every knot does not matter
    CubicSplineone cubic per interval, joined so that value, slope and second derivative agree at every knotthe curve has to be smooth: it is C², and its slope is continuous across the knots

    A spline is the piecewise cubic whose second derivatives M[i] at the knots solve h[i-1]·M[i-1] + 2(h[i-1] + h[i])·M[i] + h[i]·M[i+1] = 6·(s[i] - s[i-1]), where h is an interval's width and s its slope — n unknowns and n - 2 equations, so two more conditions are needed and SplineBoundary is which ones:

    BoundaryThe two extra conditionsEnd behaviour
    NaturalM[0] = M[n-1] = 0free to bend at both ends: the choice when nothing is known about them
    Clamped(d0, dn)S'(x0) = d0 and S'(xn) = dnprescribed slopes at the two ends
    NotAKnotthe third derivative is continuous at the second and second-to-last knotno assumption about the ends: the two outermost cubics are one polynomial

    The system is solved by solve_tridiagonal from linear.mbt, so a spline through n points costs O(n) — one pass down and one back — rather than the O(n³) of a dense elimination. The not-a-knot condition is the one that does not fit that system as a row: it relates three unknowns where row 0 has room for two, so the end curvature is eliminated with it rather than solved for, and read back once the interior ones are known.

    On the table above, the three boundaries give:

    BoundaryS(0.5)S(1.5)S(2.5)
    Natural1.950.7751.325
    Clamped(0, 0)1.650.6252.225
    NotAKnot2.43750.81250.6875

    moments() is the second derivatives M the cubics are built from — [0,-7.2, 10.8, 0] for the natural spline above — and coefficients() is the cubics themselves, one row [a, b, c, d] per interval for the polynomial a·x³ + b·x² + c·x + d in the global x rather than in the distance from the knot at the left of the interval. Two neighbouring rows of it describe the same value, the same slope and the same second derivative at the knot they share, which is what the tests assert to check the C² claim at every knot under all three boundaries.

    A table of two points is one interval, and under Natural or Clamped that is a straight line; a not-a-knot spline needs a knot that is not one, so it needs three.

    A point outside the range of the table has no answer — a spline describes the interval it was built from, not the line beyond it — so evaluate returns None there, and what cannot be built at all is raised:

    VariantPayloadRaised when
    TooFewKnots(Int)the number of pointsfewer than two points were given, or fewer than three for a not-a-knot spline
    UnevenLengths(Int, Int)the two lengthsxs and ys are not the same length
    NotSorted(Double)the first knot out of orderthe knots do not increase
    DuplicateKnots(Double)the repeated knottwo knots are the same point

    ///|
    let value = @numsolve.CubicSpline::new(table_x, table_y, @numsolve.SplineBoundary::Natural) catch {
      _ => abort("the table is a valid one")
    }.evaluate(1.5) // Some(0.775)

    The tests hold the fit to what SciPy's CubicSpline returns for the same table under the same boundary, and to the properties a spline is chosen for: every knot is reproduced exactly, the three boundaries agree on a table that lies on a straight line, and the cubics of neighbouring intervals meet in value, slope and second derivative. Sampling e^(x) at 0, 0.25, 0.5 and 1 and fitting with its true end slopes 1 and e puts the curve within 2.2e-5 of it at 0.125 and within 4.4e-4 at 0.75 — the price of a cubic per interval.

    #Command line

    The solvers are also wired to a command-line demo under cmd/main, which runs one solver on one problem and prints its answer next to the number the problem is known to have:

    moon run cmd/main -- root sqrt2
moon run cmd/main -- ode growth --steps 100
moon run cmd/main -- spline natural
moon run cmd/main -- solve hilbert

    CommandSubjectsSolverPrints
    rootsqrt2, cubic, trigbrent, false_positionthe root, against √2 or π/6
    odegrowth, decay, optionally --steps Nrk45, or rk4 with N stepsthe end value, e or 1/e, and the error
    splinenatural, clamped, notaknotCubicSpline through (0,1), (1,2), (2,0), (3,4)every piece's coefficients and three samples
    solvedemo, hilbertsolve_linear_systemthe 3×3 solution, or how close Hilbert 8's is to all ones

    --help prints the same list, and anything that is not a command is refused with a message rather than answered with a guess — the parser is tested on both. The demo lives in a package of its own, so its tests also hold the package boundary to its word: they call brent and the other solvers from outside, with and without the optional arguments, and assert the same values the library's own tests do.

    #Contributing

    Development setup, formatting, tests, benchmarks, and release checks are in CONTRIBUTING.md. Release history is in CHANGELOG.md.

    #License

    numsolve is licensed under the Apache License 2.0. See LICENSE.

    #Testing

    moon test

    104 tests. Reference values were computed first with SciPy 1.18 and NumPy 2.5 and then written down in full: SciPy's brentq on x² − 2 returns 1.4142135623731364, which is what the test holds Brent to, and the convergence rates above are read off the iterates rather than assumed. The ODE goldens come from a model of the same arithmetic — the same step sizes, the same order of operations — so the adaptive stepper's eleven steps and 77 evaluations are what the test asserts, not a range it is allowed to fall in. The vector ODE tests add analytic solutions and structural checks for the new system API. The linear tests add compensated residuals, iterative refinement, transpose and multiple-right-hand-side solves, Jacobi/Gauss-Seidel and conjugate-gradient convergence checks, general banded elimination, dense-output sampling, and event detection.

    The package also carries benchmarks, run against the same functions:

    moon bench --release

    InterpError

    source
    pub suberror InterpError {
      TooFewKnots(Int)
      UnevenLengths(Int, Int)
      NotSorted(Double)
      DuplicateKnots(Double)
    }

    A table that cannot be interpolated.

    LinError

    source
    pub suberror LinError {
      NotSquare(Int, Int)
      SizeMismatch(Int, Int)
      Singular(Int)
      NotPositiveDefinite(Int)
      NonPositiveRefinementSteps(Int)
      NonPositiveTolerance(Double)
      NonPositiveMaxIterations(Int)
      ZeroDiagonal(Int)
      NoConvergence(Int, Double)
      Breakdown(Int)
      InvalidBandwidth(Int)
      NotBanded(Int, Int)
    }

    A system that cannot be set up, or that has no single solution.

    OdeError

    source
    pub suberror OdeError {
      EmptyState
      DimensionMismatch(Int, Int)
      NonPositiveSteps(Int)
      NonPositiveTolerance(Double)
      MaxEvaluations(Int)
      OutOfRange(Double, Double, Double)
    }

    A trajectory that could not be run, or that ran out of its budget.

    RootError

    source
    pub suberror RootError {
      NotBracketed(Double, Double)
      ZeroDerivative(Double)
      NonPositiveTolerance(Double)
      MaxIterations(Int)
    }

    A root search that could not be run, or that ran out of iterations.

    Cholesky

    source
    pub struct Cholesky {
      // private fields
    }

    The Cholesky factorization A = L · Lᵀ of a symmetric positive definite matrix, as the lower factor.

    The diagonal entries are built by taking a square root, and the argument of that square root is positive exactly when the matrix is positive definite: a matrix that is symmetric and not positive definite — or not symmetric, which this does not check — has a pivot that is zero or negative, and is refused with the index of that pivot.

    The work is n³/6 against the LU's n³/3, and there is no pivoting to pay for: the property that makes the pivots positive is what makes the method stable without row exchanges.

    Cholesky::decompose

    source
    fn Cholesky::decompose(a : Array[Array[Double]]) -> Cholesky raise LinError

    Factor a symmetric positive definite matrix.

    Only the lower triangle of the argument is read: an upper triangle that disagrees with it — a matrix that is not symmetric — is ignored rather than averaged in, and the factors are those of the symmetric matrix the lower triangle describes.

    Cholesky::determinant

    source
    fn Cholesky::determinant(self : Cholesky) -> Double

    The determinant, as the square of the product of the diagonal of L: the two triangular factors have the same diagonal, so det(L · Lᵀ) is det(L)², and it is positive the way a positive definite matrix's determinant has to be.

    Cholesky::size

    source
    fn Cholesky::size(self : Cholesky) -> Int

    How wide the factored system is.

    Cholesky::solve

    source
    fn Cholesky::solve(self : Cholesky, b : Array[Double]) -> Array[Double] raise LinError

    Solve A · x = b with the factors: substitute forward through L, then back through Lᵀ, which is the same array read the other way.

    Cholesky::solve_many

    source
    fn Cholesky::solve_many(self : Cholesky, bs : Array[Array[Double]]) -> Array[Array[Double]] raise LinError

    Solve A · X = B once for each right-hand side, reusing the factor. The Cholesky factor is read by solve; this is the plural form of the same operation, useful when a model is evaluated against several loads or initial conditions.

    CubicSpline

    source
    pub struct CubicSpline {
      // private fields
    }

    A cubic spline: one cubic per interval, joined so that the value, the slope and the second derivative agree at every knot.

    CubicSpline::coefficients

    source
    fn CubicSpline::coefficients(self : CubicSpline) -> Array[Array[Double]]

    The coefficients of the cubic on each interval as a polynomial in the global x, one row [a, b, c, d] per interval of a · x³ + b · x² + c · x + d.

    Two neighbouring rows describe the same value, the same slope and the same second derivative at the knot they share — which is what makes the spline C², and what a reader checking the fit rather than sampling it wants.

    CubicSpline::evaluate

    source
    fn CubicSpline::evaluate(self : CubicSpline, x : Double) -> Double?

    The value of the spline at x, or None if x is outside the table.

    CubicSpline::moments

    source
    fn CubicSpline::moments(self : CubicSpline) -> Array[Double]

    The second derivatives at the knots, in the order the knots are given: the M of the system above, and the same numbers the cubics are built from.

    They are zero at both ends exactly when the spline is natural.

    CubicSpline::new

    source
    fn CubicSpline::new(xs : Array[Double], ys : Array[Double], boundary : SplineBoundary) -> CubicSpline raise InterpError

    Fit a cubic spline through xs and ys under the boundary condition.

    A NotAKnot spline needs three knots to have a knot that is not one: with two points it is refused, and with the other two boundaries a table of two points is a straight line, whose second derivatives are zero.

    Event

    source
    pub struct Event {
      // private fields
    }

    The first event found in a trajectory.

    Event::index

    source
    fn Event::index(self : Event) -> Int

    The index of the trajectory interval in which the event was found.

    Event::state

    source
    fn Event::state(self : Event) -> Array[Double]

    The interpolated state at the event time.

    Event::time

    source
    fn Event::time(self : Event) -> Double

    The time at which the event occurred.

    EventDirection

    source
    pub(all) enum EventDirection {
      Any
      Increasing
      Decreasing
    }

    Which sign changes count as an event.
    pub struct LU {
      // private fields
    }

    The LU factorization P · A = L · U of a square matrix, with partial pivoting.

    The two factors are held in one array: below the diagonal is L without its unit diagonal, which is not stored because it is ones, and on and above the diagonal is U. pivots holds the row permutation — row i of the permuted matrix is row pivots[i] of A — and swaps counts the exchanges, which is what the sign of the determinant turns on.

    Nothing here is exposed: a decomposition is used through solve, determinant and inverse, which are the questions it answers.

    LU::decompose

    source
    fn LU::decompose(a : Array[Array[Double]]) -> LU raise LinError

    Factor a in place — in a copy of it, that is — by Gaussian elimination with partial pivoting.

    At each column the largest remaining element of that column is exchanged into the diagonal before the column is eliminated with it. Without that, a pivot that is small makes every multiplier below it large, and the rounding of the elimination is multiplied by as much: [[1e-18, 1], [1, 1]] has the exact solution [1, 1], and elimination in the order the rows come in answers [0, 1] — the pivot is not zero, so nothing complains.

    LU::determinant

    source
    fn LU::determinant(self : LU) -> Double

    The determinant, as the product of the diagonal of U with the sign the row exchanges gave it. A singular matrix never gets this far: its decomposition raises.

    LU::inverse

    source
    fn LU::inverse(self : LU) -> Array[Array[Double]] raise LinError

    The inverse, by solving A · X = I one unit vector at a time. It is the same work as n solves, and is worth it only when the inverse is what is wanted: solving for a right-hand side is cheaper and better conditioned than multiplying by an inverse that had to be formed first.

    LU::size

    source
    fn LU::size(self : LU) -> Int

    How wide the factored system is.

    LU::solve

    source
    fn LU::solve(self : LU, b : Array[Double]) -> Array[Double] raise LinError

    Solve A · x = b with the factors: permute b the way the rows were permuted, substitute forward through L, and back through U.

    LU::solve_many

    source
    fn LU::solve_many(self : LU, bs : Array[Array[Double]]) -> Array[Array[Double]] raise LinError

    Solve A · x = b once for each right-hand side, in one pass over the factors per column: the cost of factoring is paid once and the cost of the substitutions is paid per right-hand side.

    LU::solve_refined

    source
    fn LU::solve_refined(self : LU, a : Array[Array[Double]], b : Array[Double], steps? : Int, tolerance? : Double) -> Array[Double] raise LinError

    Solve A · x = b and then improve the answer with residual-based iterative refinement.

    The factorization is computed once; each pass forms r = b - A · x, solves A · δ = r with the existing factors, and replaces x by x + δ. This is useful when an ill-conditioned direct solve has a small residual but a larger forward error, as the Hilbert systems in the tests show. The matrix passed here must be the matrix that was factored.

    steps is the number of refinement passes. A positive tolerance stops early once the infinity norm of the residual is at or below it; zero runs every requested pass.

    LU::solve_transposed

    source
    fn LU::solve_transposed(self : LU, b : Array[Double]) -> Array[Double] raise LinError

    Solve Aᵀ · x = b with the factors of A, which is what a caller wants when the transpose is the system of interest and refactoring Aᵀ would be wasted work.

    From P · A = L · U the transpose is Aᵀ = Uᵀ · Lᵀ · P. The solve substitutes through Uᵀ forward, through Lᵀ backward, and then undoes the row permutation. The factorization is not modified.

    LinearInterp

    source
    pub struct LinearInterp {
      // private fields
    }

    The straight line through a table of points.

    LinearInterp::evaluate

    source
    fn LinearInterp::evaluate(self : LinearInterp, x : Double) -> Double?

    The value of the line at x, or None if x is outside the table.

    LinearInterp::new

    source
    fn LinearInterp::new(xs : Array[Double], ys : Array[Double]) -> LinearInterp raise InterpError

    Build the line through xs and ys.

    OdeMethod

    source
    pub(all) enum OdeMethod {
      Euler
      Midpoint
      RK4
      RK45
    }

    The method used by solve_ivp.

    SplineBoundary

    source
    pub(all) enum SplineBoundary {
      Natural
      Clamped(Double, Double)
      NotAKnot
    }

    How the two ends of a cubic spline are pinned down: the two conditions that turn n - 2 equations in n unknowns into a system.

    Trajectory

    source
    pub struct Trajectory {
      // private fields
    }

    Where a stepper went: the times it stopped at and the value it found at each, with times[0] == t0 and the two arrays the same length.

    The times are the points the method actually used. The fixed-step steppers space them evenly; rk45 does not, and the values of its times are the ones its step-size control chose — that is the part of its trajectory worth reading.

    Trajectory::end_value

    source
    fn Trajectory::end_value(self : Trajectory) -> Double

    The last value, which is the answer the run was for. A trajectory always has at least one sample — the starting point — so there is one.

    Trajectory::length

    source
    fn Trajectory::length(self : Trajectory) -> Int

    How many samples the trajectory holds: one more than the number of steps taken.

    Trajectory::times

    source
    fn Trajectory::times(self : Trajectory) -> Array[Double]

    The times the trajectory was sampled at, oldest first.

    Trajectory::values

    source
    fn Trajectory::values(self : Trajectory) -> Array[Double]

    The value at each of times, in the same order.

    VectorTrajectory

    source
    pub struct VectorTrajectory {
      // private fields
    }

    A vector ODE trajectory: the times the solver stopped at and the state at each one, in matching order. The first sample is always the initial state.

    VectorTrajectory::derivatives

    source
    fn VectorTrajectory::derivatives(self : VectorTrajectory) -> Array[Array[Double]]

    The derivatives at the trajectory samples. These are stored so the trajectory can interpolate between accepted steps without asking the caller for the right-hand side again.

    VectorTrajectory::dimension

    source
    fn VectorTrajectory::dimension(self : VectorTrajectory) -> Int

    The number of components in each state vector.

    VectorTrajectory::end_state

    source
    fn VectorTrajectory::end_state(self : VectorTrajectory) -> Array[Double]

    The final state vector, which is the result of the run.

    VectorTrajectory::interpolate_at

    source
    fn VectorTrajectory::interpolate_at(self : VectorTrajectory, t : Double) -> Array[Double] raise OdeError

    Interpolate the trajectory at one time with a cubic Hermite segment.

    The value and derivative at both ends of the segment are known, so the cubic matches y and y' at the knots. Requesting a time outside the trajectory raises OutOfRange; requesting a knot returns its stored state without interpolation.

    VectorTrajectory::length

    source
    fn VectorTrajectory::length(self : VectorTrajectory) -> Int

    The number of samples in the trajectory.

    VectorTrajectory::sample_at

    source
    fn VectorTrajectory::sample_at(self : VectorTrajectory, sample_times : Array[Double]) -> Array[Array[Double]] raise OdeError

    Sample the trajectory at each requested time, in the order requested. The times need not be sorted; each one is located independently.

    VectorTrajectory::states

    source
    fn VectorTrajectory::states(self : VectorTrajectory) -> Array[Array[Double]]

    The state vectors at times, in the same order.

    VectorTrajectory::times

    source
    fn VectorTrajectory::times(self : VectorTrajectory) -> Array[Double]

    The solver times, oldest first.

    VectorTrajectory::truncate_at

    source
    fn VectorTrajectory::truncate_at(self : VectorTrajectory, t : Double) -> VectorTrajectory raise OdeError

    Truncate a trajectory at a time inside its range.

    The samples after t are removed and a final interpolated sample is added at t, so the result can represent a terminal event returned by first_event. The endpoint keeps the derivative of the Hermite segment, which means it can still be sampled or truncated again.

    bisection

    source
    fn bisection(f : (Double) -> Double, a : Double, b : Double, tolerance? : Double, max_iterations? : Int) -> Double raise RootError

    The root of f between a and b, by bisection.

    The two points have to bracket a root — their function values of opposite signs — and either of them may be the root itself. The answer is the midpoint of the last interval, so it is within tolerance · (1 + |x|) of the root, and within (b - a) / 2^k after k steps.

    It is the method to reach for when nothing is known about the function beyond its sign at two points, and it never fails on a promise that holds.

    brent

    source
    fn brent(f : (Double) -> Double, a : Double, b : Double, tolerance? : Double, max_iterations? : Int) -> Double raise RootError

    The root of f between a and b, by Brent's method.

    The two points have to bracket a root. The answer is within tolerance · (1 + |x|) of the root, and usually at the cost of far fewer steps than bisection needs for the same tolerance — six against thirty-two on the tests' functions — with no derivative to supply.

    conjugate_gradient

    source
    fn conjugate_gradient(a : Array[Array[Double]], b : Array[Double], tolerance? : Double, max_iterations? : Int) -> Array[Double] raise LinError

    Solve a symmetric positive definite system with the conjugate-gradient method.

    Each pass uses one matrix-vector product and two inner products, and the residual decreases in the A-norm. A non-positive pᵀ·A·p means the matrix is not positive definite or the iteration has broken down, and is raised as Breakdown rather than silently continued.

    euler

    source
    fn euler(f : (Double, Double) -> Double, y0 : Double, t0 : Double, t1 : Double, steps : Int) -> Trajectory raise OdeError

    Euler's method: follow the slope at the point you are standing on for the whole step.

    One evaluation per step, first-order accurate, and the method the other three are improvements on.

    euler_system

    source
    fn euler_system(f : (Double, Array[Double]) -> Array[Double], y0 : Array[Double], t0 : Double, t1 : Double, steps : Int) -> VectorTrajectory raise OdeError

    Euler's method for a vector state: one derivative evaluation per step.

    false_position

    source
    fn false_position(f : (Double) -> Double, a : Double, b : Double, tolerance? : Double, max_iterations? : Int) -> Double raise RootError

    The root of f between a and b, by false position with the Illinois modification.

    The bracket is read the same way as bisection's, and on most functions the answer arrives in a fraction of the steps: six to nine rather than thirty-odd for a tolerance of 1e-10. On a function that curves one way across the whole bracket the modification is what keeps it from stalling — the tests hold x^10 - 1 as the case in point, where the same search without the halving is still two thousandths out after a hundred steps.

    first_event

    source
    fn first_event(trajectory : VectorTrajectory, event : (Double, Array[Double]) -> Double, direction? : EventDirection, tolerance? : Double) -> Event? raise OdeError

    Find the first event in a trajectory, or None if there is no crossing.

    event(t, y) is evaluated at the accepted samples. A zero at a sample is returned directly. A sign change is refined with Brent's method on the cubic Hermite interpolation, so the event time lies inside the accepted step and no extra integration step is taken.

    gauss_seidel

    source
    fn gauss_seidel(a : Array[Array[Double]], b : Array[Double], tolerance? : Double, max_iterations? : Int) -> Array[Double] raise LinError

    Solve A · x = b by the Gauss-Seidel iteration: each component uses the values already updated in the current pass, so the method usually needs fewer iterations than Jacobi on the same matrix.

    It has the same diagonal, tolerance and iteration-budget contract as jacobi, and it raises NoConvergence when the budget runs out.

    jacobi

    source
    fn jacobi(a : Array[Array[Double]], b : Array[Double], tolerance? : Double, max_iterations? : Int) -> Array[Double] raise LinError

    Solve A · x = b by the Jacobi iteration: at each pass every component is updated from the previous iterate, so the updates can be computed in any order.

    The method converges when the diagonal is non-zero and the spectral radius of the iteration matrix is below one; a diagonally dominant matrix is the common sufficient case. A run that exhausts max_iterations raises NoConvergence rather than returning the last iterate as if it were a solution.

    midpoint

    source
    fn midpoint(f : (Double, Double) -> Double, y0 : Double, t0 : Double, t1 : Double, steps : Int) -> Trajectory raise OdeError

    The midpoint method: take a Euler step, read the slope there, and use that slope for the whole step instead.

    Two evaluations per step for second-order accuracy: the slope is sampled in the middle of the step rather than at its start, which cancels the term Euler's error comes from.

    midpoint_system

    source
    fn midpoint_system(f : (Double, Array[Double]) -> Array[Double], y0 : Array[Double], t0 : Double, t1 : Double, steps : Int) -> VectorTrajectory raise OdeError

    The explicit midpoint method for a vector state: two derivative evaluations per step and second-order accuracy.

    newton

    source
    fn newton(f : (Double) -> Double, df : (Double) -> Double, x0 : Double, tolerance? : Double, max_iterations? : Int) -> Double raise RootError

    A root of f, by Newton's method, from x0, with the derivative df.

    The starting point may be a root itself, in which case it is returned unchanged. There is no bracket, so the root it finds is the one its starting point leads to; on a function with several roots, start where the root is wanted.

    preconditioned_conjugate_gradient

    source
    fn preconditioned_conjugate_gradient(a : Array[Array[Double]], b : Array[Double], tolerance? : Double, max_iterations? : Int) -> Array[Double] raise LinError

    Solve a symmetric positive definite system with Jacobi-preconditioned conjugate gradient. The preconditioner is the diagonal of A, which costs one division per component and often reduces the iteration count of systems whose diagonal varies by orders of magnitude.

    residual_norm

    source
    fn residual_norm(a : Array[Array[Double]], b : Array[Double], x : Array[Double]) -> Double raise LinError

    The infinity norm of the residual b - A · x.

    This is a diagnostic for a computed solution, not a solver: a small residual says the equation is nearly satisfied, while the error in x also depends on the condition number of A. Iterative refinement uses the residual as a correction equation when a factored solve is not accurate enough on an ill-conditioned system.

    rk4

    source
    fn rk4(f : (Double, Double) -> Double, y0 : Double, t0 : Double, t1 : Double, steps : Int) -> Trajectory raise OdeError

    The classic fourth-order Runge-Kutta: four slopes per step — at the start, twice in the middle, at the end — averaged with weights 1/6, 2/6, 2/6, 1/6.

    The weights are Simpson's rule, which is why a right-hand side that does not depend on y is integrated exactly: a polynomial of degree three or less is what this method reproduces without error.

    rk45

    source
    fn rk45(f : (Double, Double) -> Double, y0 : Double, t0 : Double, t1 : Double, tolerance? : Double, max_evaluations? : Int) -> Trajectory raise OdeError

    The adaptive stepper: Dormand-Prince 5(4)7M, with the step size chosen by the tolerance rather than given.

    Seven evaluations of f per step produce both a fifth-order answer and a fourth-order one; their difference is the error estimate, and the step is taken only if that estimate, divided by 1 + |y|, is at or below tolerance. An accepted step is followed by one sized from the same estimate, so a right-hand side that is easy in one stretch and hard in another is spent steps where it is hard.

    The trajectory it returns is sampled where it chose to stop, which is not evenly spaced and is further apart than a fixed-step method's: for y' = y on [0, 1] a tolerance of 1e-8 reaches an error of 8e-9 in eleven steps and 77 evaluations, where rk4 needs forty steps and 160 evaluations for an error of the same size.

    max_evaluations is a budget on calls to f rather than on steps, and each step costs seven of them. Running out of it raises rather than returning a trajectory that stops short of t1.

    rk45_system

    source
    fn rk45_system(f : (Double, Array[Double]) -> Array[Double], y0 : Array[Double], t0 : Double, t1 : Double, rtol? : Double, atol? : Double, max_evaluations? : Int) -> VectorTrajectory raise OdeError

    Dormand-Prince 5(4) for a vector state, with adaptive step-size control.

    rtol and atol are both positive. Each component is scaled by atol + rtol * max(|y_old|, |y_new|), and the accepted step must have a root-mean-square scaled error no greater than one.

    rk4_system

    source
    fn rk4_system(f : (Double, Array[Double]) -> Array[Double], y0 : Array[Double], t0 : Double, t1 : Double, steps : Int) -> VectorTrajectory raise OdeError

    The classic fourth-order Runge-Kutta method for a vector state.

    secant

    source
    fn secant(f : (Double) -> Double, x0 : Double, x1 : Double, tolerance? : Double, max_iterations? : Int) -> Double raise RootError

    A root of f, by the secant method, from the two starting points x0 and x1.

    The two points say which part of the function to work in; a root is not promised between them, and a start far from a root can return a point that is not one. Either starting point may be a root itself.

    solve_banded

    source
    fn solve_banded(lower_bandwidth : Int, upper_bandwidth : Int, a : Array[Array[Double]], b : Array[Double]) -> Array[Double] raise LinError

    Solve a banded system A · x = b by elimination with partial pivoting restricted to the lower bandwidth.

    lower_bandwidth is the number of diagonals below the main diagonal that may be non-zero, and upper_bandwidth is the number above it. Entries outside that band are refused with NotBanded, because a dense matrix that violates the declared structure would otherwise be solved with the wrong arithmetic assumptions.

    solve_ivp

    source
    fn solve_ivp(f : (Double, Array[Double]) -> Array[Double], y0 : Array[Double], t0 : Double, t1 : Double, solver? : OdeMethod, steps? : Int, rtol? : Double, atol? : Double, max_evaluations? : Int) -> VectorTrajectory raise OdeError

    Solve a vector initial-value problem with a common, solve_ivp-style entry point.

    method selects a fixed-step stepper (Euler, Midpoint, RK4) or the adaptive RK45 method. Fixed-step methods read steps; RK45 reads rtol, atol, and max_evaluations. Integration backwards is allowed.

    solve_linear_system

    source
    fn solve_linear_system(a : Array[Array[Double]], b : Array[Double]) -> Array[Double] raise LinError

    Solve A · x = b by factoring A and substituting. The one-call form of LU::decompose followed by LU::solve.

    solve_linear_system_refined

    source
    fn solve_linear_system_refined(a : Array[Array[Double]], b : Array[Double], steps? : Int, tolerance? : Double) -> Array[Double] raise LinError

    Factor A, solve A · x = b, and apply iterative refinement. The one-call form of LU::decompose followed by LU::solve_refined.

    solve_tridiagonal

    source
    fn solve_tridiagonal(lower : Array[Double], diagonal : Array[Double], upper : Array[Double], rhs : Array[Double]) -> Array[Double] raise LinError

    Solve a tridiagonal system by the Thomas algorithm: eliminate the sub-diagonal forward, then substitute back.

    The matrix is given as its three diagonals rather than as a matrix — the sub-diagonal first, then the diagonal, then the super-diagonal, each with one fewer entry than the diagonal — because a tridiagonal matrix stored as rows is mostly zeros, and O(n) work on three arrays is what the shape is worth.

     | d₀ u₀  0 |      lower = [l₀, l₁]
 | l₀ d₁ u₁ |      diagonal = [d₀, d₁, d₂]
 |  0 l₁ d₂ |      upper = [u₀, u₁]

    A zero on the diagonal after elimination is a zero pivot and nothing to exchange it with within the band: the matrix is singular, and the index of that entry is what the error carries.

    Source Files

    Powered by MoonBit

    Site sourceReport issueModulesSkillsCLIBuild queueStatistics

    © 2026 mooncakes.io