Calculus

(ns calculus
  (:require [fastmath.calculus :as calc]
            [fastmath.solver :as solver]
            [fastmath.dev.codox :as codox]))

Integration

Differentiation

Solvers

Reference

fastmath.calculus

Integration and derivatives

Integrate univariate and multivariate functions.

  • VEGAS / VEGAS+ - Monte Carlo integration of multivariate function
  • h-Cubature - h-adaptive integration of multivariate function
  • Guass-Kronrod and Gauss-Legendre - quadrature integration of univariate functions
  • Romberg, Simpson, MidPoint and Trapezoid

Integrant is substituted in case of improper integration bounds.

Derivatives (finite differences method):

  • derivatives of any degree and any order of accuracy
  • gradient and hessian for multivariate functions

cubature

  • (cubature f lower upper)
  • (cubature f lower upper options)

Cubature - h-adaptive integration of multivariate function, n>1 dimensions.

Algorithm uses Genz Malik method.

In each iteration a box with biggest error is subdivided and reevaluated.

Improper integrals with infinite bounds are handled by a substitution.

Arguments:

  • f - integrant
  • lower - seq of lower bounds
  • upper - seq of upper bounds

Options:

  • :initvid - initial subdivision per dimension, default: 2.
  • :max-evals - maximum number of evaluations, default: max integer value.
  • :max-iters - maximum number of iterations, default: 64.
  • :rel - relative error, 1.0e-7
  • :abs - absolute error, 1.0e-7
  • :info? - return full information about integration, default: false

Function returns a map containing (if info? is true, returns result otherwise):

  • :result - integration value
  • :error - integration error
  • :iterations - number of iterations
  • :evaluations - number of evaluations
  • :subdivisions - final number of boxes
  • :fail? - set to :max-evals or :max-iters when one of the limits has been reached without the convergence.

derivative

  • (derivative f)
  • (derivative f n)
  • (derivative f n options)

Create nth derivative of f using finite difference method for given accuracy :acc and step :h.

Returns function.

Arguments:

  • n - derivative
  • :acc - order of accuracy (default: 2)
  • :h - step, (default: 0.0, automatic)
  • :method - :central (default), :forward or :backward
  • :extrapolate? - creates extrapolated derivative if set to true or a map with extrapolate function options

extrapolate

  • (extrapolate g)
  • (extrapolate g options)

Richardson extrapolation for given function g=g(x,h). Returns extrapolated function f(x).

Options:

  • :contract - shrinkage factor, default=1/2
  • :power - set to 2.0 for even functions around x0, default 1.0
  • :init-h - initial step h, default=1/2
  • :abs - absolute error, default: machine epsilon
  • :rel - relative error, default: ulp for init-h
  • :tol - tolerance for error, default: 2.0
  • :max-evals - maximum evaluations, default: maximum integer

f’

  • (f' f)

First central derivative with order of accuracy 2.

f’’

  • (f'' f)

Second central derivative with order of accuracy 2.

f’’’

  • (f''' f)

Third central derivative with order of accuracy 2.

fx->gx+h

  • (fx->gx+h f)

Convert f(x) to g(x,h)=f(x+h)

fx->gx-h

  • (fx->gx-h f)

Convert f(x) to g(x,h)=f(x-h)

gradient

  • (gradient f)
  • (gradient f options)

Create first partial derivatives of multivariate function for given accuracy :acc and step :h.

Returns function.

Options:

  • :acc - order of accuracy, 2 (default) or 4.
  • :h - step, default 1.0e-6

hessian

  • (hessian f)
  • (hessian f options)

Creates function returning Hessian matrix for mulitvariate function f and given :h step (default: 5.0e-3).

integrate

  • (integrate f)
  • (integrate f lower upper)
  • (integrate f lower upper {:keys [rel abs max-iters min-iters max-evals info? integrator integration-points], :or {rel BaseAbstractUnivariateIntegrator/DEFAULT_RELATIVE_ACCURACY, abs BaseAbstractUnivariateIntegrator/DEFAULT_ABSOLUTE_ACCURACY, min-iters BaseAbstractUnivariateIntegrator/DEFAULT_MIN_ITERATIONS_COUNT, max-evals Integer/MAX_VALUE, integration-points 7, integrator :gauss-kronrod, info? false}, :as options})

Univariate integration.

Improper integrals with infinite bounds are handled by a substitution.

Arguments:

  • f - integrant
  • lower - lower bound
  • upper - upper bound

Options:

  • :integrator - integration algorithm, one of: :romberg, :trapezoid, :midpoint, :simpson, :gauss-legendre and :gauss-kronrod (default).
  • :min-iters - minimum number of iterations (default: 3), not used in :gauss-kronrod
  • :max-iters - maximum number of iterations (default: 32 or 64)
  • :max-evals - maximum number of evaluations, (default: maximum integer)
  • :rel - relative error
  • :abs - absolute error
  • :integration-points - number of integration (quadrature) points for :gauss-legendre and :gauss-kronrod, default 7
  • :initdiv - initial number of subdivisions for :gauss-kronrod, default: 1
  • :info? - return full information about integration, default: false

:gauss-kronrod is h-adaptive implementation

Function returns a map containing (if info? is true, returns result otherwise):

  • :result - integration value
  • :error - integration error (:gauss-kronrod only)
  • :iterations - number of iterations
  • :evaluations - number of evaluations
  • :subdivisions - final number of boxes (:gauss-kronrod only)
  • :fail? - set to :max-evals or :max-iters when one of the limits has been reached without the convergence.

vegas

  • (vegas f lower upper)
  • (vegas f lower upper options)

VEGAS+ - Monte Carlo integration of multivariate function, n>1 dimensions.

Improper integrals with infinite bounds are handled by a substitution.

Arguments:

  • f - integrant
  • lower - seq of lower bounds
  • upper - seq of upper bounds

Additional options:

  • :max-iters - maximum number of iterations, default: 10
  • :nevals - number of evaluations per iteration, default: 10000
  • :nintervals - number of grid intervals per dimension (default: 1000)
  • :nstrats - number of stratifications per dimension (calculated)
  • :warmup - number of warmup iterations (results are used to train stratification and grid spacings, default: 0
  • :alpha - grid refinement parameter, 0.5 slow (default for vegas+), 1.5 moderate/fast (defatult for vegas)
  • :beta - stratification damping parameter for startification adaptation, default: 0.75
  • :rel - relative accuracy, default: 5.0e-4
  • :abs - absolute accuracy, default: 5.0e-4
  • :random-sequence - random sequence used for generating samples: :uniform (default), low-discrepancy sequences: :r2, :sobol and :halton.
  • :jitter - jittering factor for low-discrepancy random sequence, default: 0.75
  • :info? - return full information about integration, default: false
  • :record-data? - stores samples, number of strata, x and dx, default: false (requires, :info? to be set to true)

For original VEGAS algorithm set :nstrats to 1.

:nstrats can be also a list, then each dimension is divided independently according to a given number. If list is lower then number of dimensions, then it’s cycled.

Function returns a map with following keys (if info? is true, returns result otherwise):

  • :result - value of integral
  • :iterations - number of iterations (excluding warmup)
  • :sd - standard deviation of results
  • :nintervals - actual grid size
  • :nstrats - number of stratitfications per dimension
  • :nhcubes - number of hypercubes
  • :evaluations - number of function calls
  • :chi2-avg - average of chi2
  • :dof - degrees of freedom
  • :Q - goodness of fit indicator, 1 - very good, <0.25 very poor
  • :data - recorded data (if available)

fastmath.solver

cubic

  • (cubic a b c d)

Real solution of cubic formula ax3+bx2+cx+d=0. Returns NaNs where no real solutions are found.

find-root

  • (find-root f lower-bound upper-bound)
  • (find-root f lower-bound upper-bound {:keys [absolute-accuracy relative-accuracy max-iters initial-value solver], :or {max-iters 100, solver :brent}})

Find zero (root) of a function f in given range [lower-bound, upper-bound].

Optional parameters:

  • :absolute-accuracy - default 1.0e-8
  • :relative-accuracy
  • :max-iters - maximum iterations (default: 100)
  • :initial-value - algorithm starting value
  • :solver - one of: :brent (default), :bisection, :illinois, :muller, :muller2, :pegasus, :regula-falsi, :ridders and :secant.

quadratic

  • (quadratic a b c)

Real solutions of quadratic formula ax^2+bx+c=0. Returns NaNs where no real solutions are found.