viva_math/calculus

Numerical calculus: differentiation and integration.

Closed-form symbolic calculus is out of scope. This module provides finite-difference approximations for derivatives and quadrature rules for integrals.

Order of accuracy

Function	Order	Notes
`forward_diff` / `backward`	O(h)	One-sided
`central_diff`	O(h²)	Recommended default
`five_point_diff`	O(h⁴)	Five evaluations, very accurate
`second_derivative`	O(h²)	Central stencil
`trapezoid` / `simpson`	O(h²/h⁴)	Composite rules
`romberg`	adaptive	Richardson extrapolation

Values

backward_diff

</>

pub fn backward_diff(
  f: fn(Float) -> Float,
  x: Float,
  h: Float,
) -> Float

Backward difference (f(x) - f(x-h)) / h. O(h) accurate.

central_diff

</>

pub fn central_diff(
  f: fn(Float) -> Float,
  x: Float,
  h: Float,
) -> Float

Central difference (f(x+h) - f(x-h)) / (2h). O(h²) accurate.

five_point_diff

</>

pub fn five_point_diff(
  f: fn(Float) -> Float,
  x: Float,
  h: Float,
) -> Float

Five-point stencil O(h⁴): the gold standard for smooth functions.

f’(x) ≈ (-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)) / (12h)

forward_diff

</>

pub fn forward_diff(
  f: fn(Float) -> Float,
  x: Float,
  h: Float,
) -> Float

Forward difference (f(x+h) - f(x)) / h. O(h) accurate.

gradient

</>

pub fn gradient(
  f: fn(List(Float)) -> Float,
  point: List(Float),
  h: Float,
) -> List(Float)

Gradient of a multivariate function f: List(Float) -> Float.

Uses central differences along each axis with step h.

norm_l2

</>

pub fn norm_l2(values: List(Float)) -> Float

romberg

</>

pub fn romberg(
  f: fn(Float) -> Float,
  a: Float,
  b: Float,
  levels: Int,
) -> Float

Romberg integration via Richardson extrapolation.

Repeatedly halves h and combines trapezoid estimates to cancel error terms. levels controls the depth (typically 5-8 suffices).

second_derivative

</>

pub fn second_derivative(
  f: fn(Float) -> Float,
  x: Float,
  h: Float,
) -> Float

Second derivative via central stencil. O(h²).

f’’(x) ≈ (f(x+h) - 2f(x) + f(x-h)) / h²

simpson

</>

pub fn simpson(
  f: fn(Float) -> Float,
  a: Float,
  b: Float,
  n: Int,
) -> Result(Float, Nil)

Composite Simpson’s rule. n must be even and ≥ 2.

∫f ≈ (h/3) · (f₀ + 4f₁ + 2f₂ + 4f₃ + … + 4fₙ₋₁ + fₙ)

trapezoid

</>

pub fn trapezoid(
  f: fn(Float) -> Float,
  a: Float,
  b: Float,
  n: Int,
) -> Float

Composite trapezoid rule over n equal subintervals.

∫f ≈ h · ( ½f(a) + Σ f(a + i·h) + ½f(b) ), h = (b - a) / n