CODEWITHSUNDEEP
haskellintegration
bubu
bubu

Reputation: 179

Trapezoid Rule in Haskell

I'm trying to define the trapezoid rule in Haskell. I created a helper function innerSum, which is sum portion of the trapezoid rule. Then in the integral definition I multiple by the distance and take in the lower bound, upper bound, a function, and some n number of trapezoids.As n increases the answer should become more accurate. My function seems to work for most cases

except this case (and likely others): definiteIntegral (-1) 1 (\x->x^100) 20.

As I change the value for 20 instead of my answers diverging to a certain and getting more accurate, the numbers just jump randomly. I cannot seem to find the mistake

definiteIntegral :: Double -> Double -> (Double -> Double) -> Integer -> Double
definiteIntegral a b g n | a <= b = (dist a b n)*(innerSum a b g (dist a b n))
                         | otherwise = (dist a b n)*(innerSum b a g (dist b a n))
                         where dist a b n = (b-a)/(fromIntegral n::Double)

innerSum :: Double -> Double -> (Double -> Double) -> Double -> Double
innerSum a b g d | a >= b = 0
                 | otherwise = (((g a) + (g (a + d)))/2)+(innerSum (a + d) b g d)

Upvotes: 1

Views: 632

Answers (1)

leftaroundabout
leftaroundabout

Reputation: 120731

The problem is your end condition. You calculate the step size so in principle, you should end up exactly at the right boundary after n steps. And indeed that will work if you make sure the step size can be exactly represented in floating-point, i.e. a multiple of a power of two:

*Main> definiteIntegral (-1) 1 (^100) <$> [64, 128, 256, 512, 1024]
[3.396429282002939e-2,2.3718601030590182e-2,2.08093271780123e-2,2.0055667986086628e-2,1.9865519301509465e-2]
     -- convergence towards 1.98×10⁻²

But most numbers can not be exactly represented in FP, therefore you'll regularly not hit the right boundary exactly. If the last point falls slightly short of that bound, the algorithm will do an entire extra step, i.e. you blow up a small float inaccuracy to an entire-step inaccuracy. Numerically speaking, it's “just” an order-1 error, because that single step size gets ever smaller as you increase the resolution. (It's still bad, because the trapezoidal rule should actually be second-order accurate!)

Problem with (^100) in particular is: that function is very close to zero on most of the interval [-1,1], but grows very quickly in the vicinity of ±1. Therefore the result of the interval is dominated by the trapezes right next to the boundary, and even more by a trapeze outside the boundary!

The easiest fix is to get rid of that innerSum function and use built-in tools instead:

definiteIntegral :: Double -> Double -> (Double -> Double) -> Integer -> Double
definiteIntegral a b g n
    | a <= b = - definiteIntegral b a g n
    | otherwise = δ * sum [ (g x + g (x+δ))/2
                          | x <- [a, a+δ .. b] ]
 where δ = (b-a) / fromIntegral n

It is explained in this question why that works reliable even when the step size is not exactly represented: float ranges like [a, a+δ .. b] refuse to take an extra step if the previous one is already close to the boundary and the next one would take you far beyond the boundary.

Upvotes: 5

Related Questions