Haskell numerical integration via Trapezoidal rule results in wrong sign

Question

I've written some code that's meant to integrate a function numerically using the trapezoidal rule. It works, but the answer it produces has a wrong sign. Why might that be?

The code is:

integration :: (Double -> Double) -> Double -> Double -> Double
integration f a b = h * (f a + f b + partial_sum)
    where 
        h = (b - a) / 1000 
        most_parts  = map f (points (1000-1) h) 
        partial_sum = sum most_parts

points  :: Double -> Double -> [Double]
points x1 x2 
    | x1 <= 0 = []
    | otherwise = (x1*x2) : points (x1-1) x2

Trapezoidal rule

The code is probably inelegant, but I'm only a student of Haskell and would like to deal with the current problem first and coding style matters after that.

Answer 1

_{Note: This answer is written in literate Haskell. Save it with .lhs as extension and load it in GHCi to test the solution.}

Finding the culprit

First of all, let's take a look at integration. In its current form, it contains only summation of function values f x. Even though the factors aren't correct at the moment, the overall approach is fine: you evaluate f at the grid points. However, we can use the following function to verify that there's something wrong:

ghci> integration (\x -> if x >= 10 then 1 else (-1)) 10 15
-4.985

Wait a second. x isn't even negative in [10,15]. This suggests that you use the wrong grid points.

Grid points revisited

Even though you've linked the article, let's have a look at an exemplary use of the trapezoidal rule (public domain, original file by Oleg Alexandrov):

Although this doesn't use a uniform grid, let's suppose that the 6 grid points are equidistant with grid distance h = (b - a) / 5. What are the x coordinates of those points?

x_0 = a + 0 * h (== a)
x_1 = a + 1 * h
x_2 = a + 2 * h
x_3 = a + 3 * h
x_4 = a + 4 * h
x_5 = a + 5 * h (== b)

If we use set a = 10 and b = 15 (and therefore h = 1), we should end up with [10, 11, 12, 13, 14, 15]. Let's check your points. In this case, you would use points 5 1 and end up with [5,4,3,2,1].

And there's the error. points doesn't respect the boundary. We can easily fix this by using pointsWithOffset:

> points  :: Double -> Double -> [Double]
> points x1 x2 
>     | x1 <= 0 = []
>     | otherwise = (x1*x2) : points (x1-1) x2
>
> pointsWithOffset :: Double -> Double -> Double -> [Double]
> pointsWithOffset x1 x2 offset = map (+offset) (points x1 x2)

That way, we can still use your current points definition to generate grid points from x1 to 0 (almost). If we use integration with pointsWithOffset, we end up with

integration :: (Double -> Double) -> Double -> Double -> Double
integration f a b = h * (f a + f b + partial_sum)
    where 
        h = (b - a) / 1000 
        most_parts  = map f (pointsWithOffset (1000-1) h a)  
        partial_sum = sum most_parts

Tying up loose ends

However, this doesn't take into account that you use all inner points twice in the trapezoid rule. If we add the factors, we end up with

> integration :: (Double -> Double) -> Double -> Double -> Double
> integration f a b = 
>      h / 2 * (f a + f b + 2 * partial_sum)
>    --    ^^^              ^^^
>    where 
>        h = (b - a) / 1000 
>        most_parts  = map f (pointsWithOffset (1000-1) h a)  
>        partial_sum = sum most_parts

Which yields the correct value for our test function above.

Exercise

Your current version only supports 1000 grid points. Add an Int argument so that one can change the number of grid points:

integration :: Int -> (Double -> Double) -> Double -> Double -> Double
integration n f a b = -- ...

Furthermore, try to write points in different ways, for example go from a to b, use takeWhile and iterate, or even a list comprehension.

Answer 2

Yes it indeed was the points plus you had some factors wrong (the inner points are multiplied by 2) - this is the fixed version of your code:

integration :: (Double -> Double) -> Double -> Double -> Double
integration f a b = h * (f a + f b + innerSum) / 2
    where 
        h = (b - a) / 1000 
        innerPts  = map ((2*) . f . (a+)) (points (1000-1) h) 
        innerSum = sum innerPts

points  :: Double -> Double -> [Double]
points i x 
    | i <= 0 = []
    | otherwise = (i*x) : points (i-1) x

which gives sensible approximations (to 1000 points):

λ> integration (const 2) 1 2
2.0
λ> integration id 1 2
1.5

λ> integration (\x -> x*x) 1 2
2.3333334999999975

λ> 7/3
2.3333333333333335