diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,3 +5,8 @@ 1.0.0.1 ------- * removed the upper bounds of the dependencies + +1.1.0.0 +------- +* new function `integratePolynomialOnSimplex`, returning the exact value of + the integral of a polynomial over a simplex diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -2,6 +2,18 @@ Pure Haskell implementation of simplicial cubature (integration over a simplex). +This library is a port of a part of the R package **SimplicalCubature**, +written by John P. Nolan, and which contains R translations of +some Matlab and Fortran code written by Alan Genz. +It is also a port of a part of the R package **SphericalCubature**, also +written by John P. Nolan. +In addition it provides a function for the exact computation of the integral +of a polynomial over a simplex. + +___ + +## Integral of a function on a simplex + ```haskell integrateOnSimplex :: (VectorD -> VectorD) -- integrand @@ -14,7 +26,7 @@ -> IO Results -- values, error estimates, evaluations, success ``` -## Example +### Example ![equation](http://latex.codecogs.com/gif.latex?%5Cint_0%5E1%5Cint_0%5Ex%5Cint_0%5Ey%5Cexp%28x+y+z%29%5C,%5Cmathrm%7Bd%7Dz%5C,%5Cmathrm%7Bd%7Dy%5C,%5Cmathrm%7Bd%7Dx=%5Cfrac%7B1%7D%7B6%7D%28e-1%29%5E3%5Capprox%20.8455356853) @@ -28,7 +40,7 @@ :} ``` -Define the simplex: +Define the simplex (tetrahedron in dimension 3) by the list of its vertices: ```haskell simplex = [[0, 0, 0], [1, 1, 1], [0, 1, 1], [0, 0, 1]] @@ -65,11 +77,61 @@ -- , success = True } ``` + +## Exact integral of a polynomial on a simplex + +```haskell +integratePolynomialOnSimplex + :: (C a, Fractional a, Ord a) -- `C a` means that `a` must be a ring + => Spray a -- ^ polynomial to be integrated + -> [[a]] -- ^ simplex to integrate over + -> a +``` + +### Example + +We take as an example the rational numbers for `a`. Thus we must take a +polynomial with rational coefficients and a simplex whose vertices +have rational coordinates. Then the integral will be a rational number. + +Our polynomial is + +![equation](https://latex.codecogs.com/gif.image?\dpi{110}P(x,&space;y,&space;z)&space;=&space;x^4&space;+&space;y&space;+&space;2xy^2&space;-&space;3z.) + +It must be defined in Haskell with the +[**hspray**](https://github.com/stla/hspray) library. + +```haskell +import Numeric.Integration.IntegratePolynomialOnSimplex +import Data.Ratio +import Math.Algebra.Hspray + +:{ +simplex :: [[Rational]] +simplex = [[1, 1, 1], [2, 2, 3], [3, 4, 5], [3, 2, 1]] +:} + +x = lone 1 :: Spray Rational +y = lone 2 :: Spray Rational +z = lone 3 :: Spray Rational + +:{ +poly :: Spray Rational +poly = x^**^4 ^+^ y ^+^ 2.^(x ^*^ y^**^2) ^-^ 3.^z +:} + +integratePolynomialOnSimplex poly simplex +-- 1387 % 42 +``` + + ## 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). +### Example + For example take the first orthant in dimension 3: ```haskell @@ -99,3 +161,14 @@ -- , evaluations = 17065 -- , success = True } ``` + + +## References + +- A. Genz and R. Cools. + *An adaptive numerical cubature algorithm for simplices.* + ACM Trans. Math. Software 29, 297-308 (2003). + +- Jean B. Lasserre. + *Simple formula for the integration of polynomials on a simplex.* + BIT Numerical Mathematics 61, 523-533 (2021). diff --git a/scubature.cabal b/scubature.cabal --- a/scubature.cabal +++ b/scubature.cabal @@ -1,6 +1,6 @@ cabal-version: >=1.10 name: scubature -version: 1.0.0.1 +version: 1.1.0.0 license: GPL-3 license-file: LICENSE copyright: 2022 Stéphane Laurent @@ -25,6 +25,7 @@ exposed-modules: Numeric.Integration.SimplexCubature Numeric.Integration.SphericalSimplexCubature + Numeric.Integration.IntegratePolynomialOnSimplex hs-source-dirs: src other-modules: @@ -41,4 +42,6 @@ vector >=0.12.3.1, matrix >=0.3.6.1, containers >=0.6.4.1, - ilist >=0.4.0.1 + ilist >=0.4.0.1, + hspray >=0.1.1.0, + numeric-prelude >=0.4.4 diff --git a/src/Numeric/Integration/IntegratePolynomialOnSimplex.hs b/src/Numeric/Integration/IntegratePolynomialOnSimplex.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Integration/IntegratePolynomialOnSimplex.hs @@ -0,0 +1,43 @@ +module Numeric.Integration.IntegratePolynomialOnSimplex where +import Algebra.Ring ( C ) +import Data.List ( transpose ) +import Data.Matrix ( detLU + , fromLists + ) +import Math.Algebra.Hspray ( (*^) + , Spray + , (^+^) + , bombieriSpray + , composeSpray + , constantSpray + , lone + , toList + ) + +-- | Exact integral of a polynomial over a simplex +integratePolynomialOnSimplex + :: (C a, Fractional a, Ord a) + => Spray a -- ^ polynomial to be integrated + -> [[a]] -- ^ simplex to integrate over + -> a +integratePolynomialOnSimplex p simplex = + s * abs (detLU $ fromLists b) / (fromIntegral $ product [2 .. n]) + where + v = last simplex + n = length v + b = map (\column -> zipWith (-) column v) (take n simplex) + vb = zip v (transpose b) + variables = map lone [1 .. n] + newvariables = map + (\(vi, bi) -> + (constantSpray vi) ^+^ foldl1 (^+^) (zipWith (*^) bi variables) + ) + vb + q = composeSpray p newvariables + qterms = toList $ bombieriSpray q + s = sum $ map f qterms + where + f (exponents, coef) = if d == 0 + then coef + else coef / (fromIntegral $ product [n + 1 .. n + d]) + where d = sum exponents diff --git a/src/Numeric/Integration/SimplexCubature.hs b/src/Numeric/Integration/SimplexCubature.hs --- a/src/Numeric/Integration/SimplexCubature.hs +++ b/src/Numeric/Integration/SimplexCubature.hs @@ -6,7 +6,7 @@ import Data.Array.Unsafe (unsafeThaw) import qualified Data.Vector.Unboxed as UV import Numeric.Integration.Simplex.Simplex ( isValidSimplices, Simplices ) -import Numeric.Integration.SimplexCubature.Internal ( VectorD, IO3dArray, adsimp ) +import Numeric.Integration.SimplexCubature.Internal ( VectorD, IO3dArray, adsimp ) data Results = Results { values :: [Double]