Data.Massiv.Array.Numeric.Integral

Inspiration for the code in this module was taken from Paul Dawkins Online Notes. In particular the page on Integral Approximation, so if you need to brush up on some theory it is a great place to start.

Implementation-wise, integral approximation here relies heavily on stencils with stride, because such computation is fast and is automatically parallelizable.

Here are some examples of where this can be useful:

Integral of a function on a region

Say we have a gaussian f(x) = e^(x^2) on interval [0, 2] (as in Paul's tutorial above). For this we define a function f, an array with equally spaced (dx) sample input values and apply the function to that array, which will give us an array of n + 1 sample points, or looking from a different angle n intervals.

>>> import Data.Massiv.Array>>> f x = exp ( x ^ (2 :: Int) ) :: Float>>> fromFunction Seq (\ scale x -> f (scale x)) 0 2 (Sz1 1) 4Array D Seq (Sz1 5)  [ 1.0, 1.2840254, 2.7182817, 9.487736, 54.59815 ]

Once we have that array of sample points ready, we could use integralApprox and one of the stencils to compute an integral, but there are already functions that will do both steps for you:

>>> simpsonsRule Seq U (\ scale x -> f (scale x)) 0 2 (Sz1 1) 4Array D Seq (Sz1 1)  [ 17.353626 ]

scale is the function that will change an array index into equally spaced and appropriately shifted values of x, y, ... before they can get applied to f(x, y, ...)

Accurate values of a function

Another very useful place where integral approximation can be used is when a more accurate representation of a non-linear function is desired. Consider the same gaussian function applied to equally spaced values, with zero being in the middle of the vector:

>>> xArr = makeArrayR D Seq (Sz1 4) $ \ i -> fromIntegral i - 1.5 :: Float>>> xArrArray D Seq (Sz1 4)  [ -1.5, -0.5, 0.5, 1.5 ]>>> fmap f xArrArray D Seq (Sz1 4)  [ 9.487736, 1.2840254, 1.2840254, 9.487736 ]

The problem with above example is that computed values do not accurately represent the total value contained within each vector cell. For that reason if your were to later use it for example as convolution stencil, approximation would be very poor. The way to solve it is to approximate an integral across each cell of vector by drastically blowing up the xArr and then reducing it to a smaller array by using one of the approximation rules:

>>> startValue = -2 :: Float>>> distPerCell = 1 :: Float>>> desiredSize = Sz1 4 :: Sz1>>> numSamples = 4 :: Int>>> xArrX4 = fromFunction Seq ($) startValue distPerCell desiredSize numSamples>>> xArrX4Array D Seq (Sz1 17)  [ -2.0, -1.75, -1.5, -1.25, -1.0, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0 ]>>> yArrX4 = computeAs U $ fmap f xArrX4>>> integralApprox trapezoidStencil distPerCell desiredSize numSamples yArrX4Array D Seq (Sz1 4)  [ 16.074406, 1.4906789, 1.4906789, 16.074408 ]

We can clearly see the difference is huge, but it doesn't mean it is much better than our previous estimate. In order to get more accurate results we can use a better Simpson's rule for approximation and many more sample points. There is no need to create individual arrays xArrX4 and yArrX4, there are functions like simpsonsRule that will take care of it for us:

>>> simpsonsRule Seq U (\ scale i -> f (scale i)) startValue distPerCell desiredSize 128Array D Seq (Sz1 4)  [ 14.989977, 1.4626511, 1.4626517, 14.989977 ]