github ↗
docs menu
on this page

Module Jera.Quad

Integrals.

An integrand is an elementwise function: it receives points with the rule's node axes in front and the range's shape behind, and returns one value per point. Over a range every element is its own integral, so the integrand must not reduce or mix along any axis; over a box every lane is its own integral, and the integrand reduces only the last, coordinate axis. The axes in front and their sizes are the method's, and may change between calls.

ocaml
(* ∫₀¹ x^a dx for each a, by 10-point Gauss–Legendre *)
let moments a =
  Quad.fixed (Quad.Rule.gauss 10)
    (fun x -> Nx.pow x a)
    (Quad.Range.v (Nx.zeros_like a) (Nx.ones_like a))

An integral of samples is the integral of their interpolant (Piecewise.integral).

Rules

#module Rule : sig ... end

Quadrature rules on [−1, 1].

Ranges

#module Range : sig ... end

Ranges of integration, one per element.

#module Box : sig ... end

Boxes of integration, one per lane.

#type ('b, 'd) integrand = (float, 'b, 'd) Nx.t -> (float, 'b, 'd) Nx.t

The type for integrands: f x is f at each point of x: of x's shape over a range, and over a box of x's shape without its last, coordinate axis.

Formulas

#val fixed : [> `Formula ] Rule.t -> ('b, 'd) integrand -> ('b, 'd) Range.t -> (float, 'b, 'd) Nx.t

fixed r f range is r's sum for the integral of f over each element of range, of range's shape. f receives the points of shape [m] @ shape, m the rule's nodes. A finite range maps the rule linearly; from a by x = a + (1 + u) / (1 − u) and line c by x = c + u / (1 − u²), so an infinite range needs an integrand that decays fast enough for the rule's degree to show.

Error. An n-point Gauss sum on a finite range is exact for polynomials of degree 2n − 1, and for a smooth f its error falls geometrically in n. Cost. One call of f on m points per element. Derivative. The sum's: in the integrand's parameters, in the ends and in the points f reads.

Raises Invalid_argument if f's result has another shape than its points.

#val cumulative : [> `Formula ] Rule.t -> ('b, 'd) integrand -> (float, 'b, 'd) Nx.t -> (float, 'b, 'd) Nx.t

cumulative r f knots is the integral of f from the first knot to each, by r on each interval between knots: for knots of shape [n] @ shape, a result of the same shape whose row 0 is zero and row i the sum of the first i intervals' integrals. f receives points of shape [m; n − 1] @ shape.

Raises Invalid_argument if knots has no axis or no knot, or as fixed does.

Solves

Each solve is elementwise: every element of the range is its own integral, with its own status. Its search runs on detached values; the answer of a converged element is its rule over its final decisions, tracked, so its derivative is that rule's, and an element that did not converge returns its detached best estimate.

#val adaptive : [> `Embedded ] Rule.t -> tol:Tol.t -> budget:int -> ('b, 'd) integrand -> ('b, 'd) Range.t -> ((float, 'b, 'd) Nx.t, 'd) Solution.t

adaptive r ~tol ~budget f range is the integral of f over each element of range. f receives points of shape [m; k] @ shape, m the rule's points and k pieces.

Method. The rule r on a partition it refines: it bisects the piece of largest error until the error meets tol. An infinite range is mapped as fixed maps it. Error. e is the sum over the pieces of |K − G|, the Kronrod sum's difference from its embedded Gauss sum, and y the integral. An element whose worst piece is at level 62, or holds no float strictly inside, ends Stalled; a non-finite sum ends it Not_finite; budget pieces end it Budget_spent. A feature narrower than the first rule's nodes can be invisible to every estimate, and an element can converge without it. Cost. 2n + 1 points per piece, and the answer evaluates the final partition again, in chunks of 32 pieces under a Rune.scan, so reverse mode keeps one chunk's values. Derivative. The final partition's rule's: the pieces are integers (level, index) whose ends are a + (b − a) index / 2^level, so it reaches the ends.

Raises Invalid_argument if budget < 1, or as fixed does.

#val tanh_sinh : tol:Tol.t -> ('b, 'd) integrand -> ('b, 'd) Range.t -> ((float, 'b, 'd) Nx.t, 'd) Solution.t

tanh_sinh ~tol f range is the integral of f over each element of range by a double-exponential rule: tanh-sinh on a finite range, exp-sinh on Range.from and sinh-sinh on Range.line. Its nodes crowd toward the ends double-exponentially, so it converges for an integrable singularity at an end, and on a half-line or the line for an integrand that decays. Scale the variable so the integrand's width is near 1. f receives points of shape [64] @ shape, a chunk of nodes.

Method. The trapezoidal rule in t after the change of variable, at steps 1, 1/2, 1/4, ...: each level adds the odd multiples of its step, in fixed-size chunks, to the finest level the dtype calls for (2^-7 in float64, 2^-6 in float32). The nodes stop where their numbers leave the dtype's normal floats; a node whose point rounds to an end is unused. Near an end a = 0 a point is its own distance to the end, so the nodes reach the singularity in full precision; another end loses the digits a's magnitude rounds away, unless the integrand computes the distance to the end from its argument. Error. e is the difference of the last two levels and y the integral. A lane whose terms at the truncation do not fall below the tolerance, or whose finest level does not meet it, ends Stalled; a non-finite sum ends it Not_finite. Cost. At most the nodes of the finest level, in chunks of 64, per element, and the answer evaluates the final level's nodes again. Derivative. That of the sum over the final level, through the ends and the integrand's parameters.

#val cubature : tol:Tol.t -> budget:int -> ('b, 'd) integrand -> ('b, 'd) Box.t -> ((float, 'b, 'd) Nx.t, 'd) Solution.t

cubature ~tol ~budget f box is the integral of f over each lane's box, of d dimensions with 2 ≤ d ≤ 10. f receives points of shape [m; k] @ lanes @ [d], m the rule's points and k boxes, and reduces only their last, coordinate axis.

Method. Genz and Malik's (1980) adaptive rule of degree 7 with an embedded degree 5: it bisects the box of largest error across the axis of largest fourth difference. Error. e is the sum over the boxes of the two rules' difference, and y the integral. A lane whose worst box is at level 62 along its axis, or holds no float strictly inside along it, ends Stalled; a non-finite sum ends it Not_finite; budget boxes end it Budget_spent. Cost. 2^d + 2d² + 2d + 1 points per box, budget boxes at most per lane, and the answer evaluates the final partition again. Derivative. The final partition's rule's, through the corners and the integrand's parameters.

Raises Invalid_argument if d is not in [2, 10], if budget < 1, or if f's result is not the points' shape without its last axis.

#val qmc : 'd Nx.Rng.key -> tol:Tol.t -> budget:int -> ('b, 'd) integrand -> ('b, 'd) Box.t -> ((float, 'b, 'd) Nx.t, 'd) Solution.t

qmc key ~tol ~budget f box is the integral of f over each lane's box, of any dimension d up to 1111, by randomised quasi-Monte Carlo. f receives points of shape [64; 16] @ lanes @ [d], a chunk of 64 points under each of 16 shifts, and reduces only their last, coordinate axis.

Method. The mean of f over a Sobol sequence (Joe and Kuo's direction numbers) under 16 independent random digital shifts drawn from key. It adds the sequence in chunks of 64 points and tests at each power of two, where a Sobol prefix is balanced. Points are (i + ½) / 2^k after the shift, k the bits the dtype holds below 1 (32 in float64), so none lies on the box's boundary. Error. e is the standard error of the mean over the shifts, an estimate of a standard deviation: the test is statistical. y is the integral. Each estimate at a fixed point count is unbiased, and the stopped one to within its standard error. budget chunks end a lane Budget_spent. Cost. 1024 points per chunk, 64 under each of 16 shifts, budget chunks at most, and the answer evaluates the final points again. Derivative. The mean's over the final points: an estimate of the integral's derivative where the integrand is Lipschitz in the parameter.

Raises Invalid_argument if d is above 1111, if budget < 1, or if f's result is not the points' shape without its last axis.