scubature: Multidimensional integration over simplices

[ gpl, integration, library, numeric ] [ Propose Tags ]

This library allows to evaluate integrals over Euclidean and spherical simplices.


[Skip to Readme]

Downloads

Maintainer's Corner

For package maintainers and hackage trustees

Candidates

Versions [RSS] 1.0.0.0
Change log CHANGELOG.md
Dependencies array (>=0.5.4.0 && <0.6), base (>=4.7 && <5), containers (>=0.6.4.1 && <0.7), ilist (>=0.4.0.1 && <0.5), matrix (>=0.3.6.1 && <0.4), vector (>=0.12.3.1 && <0.13) [details]
License GPL-3.0-only
Copyright 2022 Stéphane Laurent
Author Stéphane Laurent
Maintainer laurent_step@outlook.fr
Category Numeric, Integration
Home page https://github.com/stla/scubature#readme
Source repo head: git clone https://github.com/stla/scubature
Uploaded by stla at 2022-08-06T15:00:39Z
Distributions NixOS:1.0.0.0
Downloads 18 total (2 in the last 30 days)
Rating (no votes yet) [estimated by Bayesian average]
Your Rating
  • λ
  • λ
  • λ
Status Docs available [build log]
Last success reported on 2022-08-06 [all 1 reports]

Readme for scubature-1.0.0.0

[back to package description]

scubature

Pure Haskell implementation of simplicial cubature (integration over a simplex).

integrateOnSimplex
    :: (VectorD -> VectorD)   -- integrand
    -> Simplices              -- domain of integration (union of the simplices)
    -> Int                    -- number of components of the integrand
    -> Int                    -- maximum number of evaluations
    -> Double                 -- desired absolute error
    -> Double                 -- desired relative error
    -> Int                    -- integration rule: 1, 2, 3 or 4
    -> IO Results             -- values, error estimates, evaluations, success

Example

equation

Define the integrand:

import Data.Vector.Unboxed as V
:{
f :: Vector Double -> Vector Double
f v = singleton $ exp (V.sum v)
:}

Define the simplex:

simplex = [[0, 0, 0], [1, 1, 1], [0, 1, 1], [0, 0, 1]]

Integrate:

import Numeric.Integration.SimplexCubature
integrateOnSimplex f [simplex] 1 100000 0 1e-10 3
-- Results { values = [0.8455356852954488]
--         , errorEstimates = [8.082378899762402e-11]
--         , evaluations = 8700
--         , success = True }

For a scalar-valued integrand, it's more convenient to define... a scalar-valued integrand! That is:

:{
f :: Vector Double -> Double
f v = exp (V.sum v)
:}

and then to use integrateOnSimplex':

integrateOnSimplex' f [simplex] 100000 0 1e-10 3
-- Result { value         = 0.8455356852954488
--        , errorEstimate = 8.082378899762402e-11
--        , evaluations   = 8700
--        , success       = True }

Integration on a spherical triangle

The library also allows to evaluate an integral on a spherical simplex on the unit sphere (in dimension 3, a spherical triangle).

For example take the first orthant in dimension 3:

import Numeric.Integration.SphericalSimplexCubature
o1 = orthants 3 !! 0
o1
-- [ [1.0, 0.0, 0.0]
-- , [0.0, 1.0, 0.0]
-- , [0.0, 0.0, 1.0] ]

And this integrand:

:{
integrand :: [Double] -> Double
integrand x = (x!!0 * x!!0 * x!!2) + (x!!1 * x!!1 * x!!2) + (x!!2 * x!!2 * x!!2)
:}

Compute the integral (the exact result is pi/4 ≈ 0.7853981634):

integrateOnSphericalSimplex integrand o1 20000 0 1e-7 3
-- Result { value         = 0.7853981641913279
--        , errorEstimate = 7.71579524444753e-8
--        , evaluations   = 17065
--        , success       = True }